Thermodynamically, native cellulose nano fibrils are more stable in an aggregated state. The aggregated state is however not useful from a material development perspective. Therefore much research has been done to stabilize the dispersal of the fibrils. One method to overcome this instability is by surface substitution of the O6 hydroxyl group with carboxylate groups, to make highly charged fibrils in aqueous solutions. It is therefore of much interest to understand the interaction of highly charged fibrils in aqueous solutions. In this study, we aim to model the interaction potential energy between native and surface modified cellulose nanofibrils in order to understand under what conditions the contribution from the dipole interactions can be neglected. To achieve this we propose to use a continuum electrostatic approach, modeling the electrostatic interactions as a function of the fibrils relative dipole orientation, separation, surface charge as well as ionic strength of the solution, by means of using the Poisson-Boltzmann equation.