The hierarchical self-assembly of cellulose nanocrystals (CNCs) is an important phenomenon occurring naturally in plant cell walls. Utilization of this assembly for advanced applications requires a fundamental theoretical understanding of interactions between the CNCs, which is still incomplete. Hence, in this work, we used molecular dynamics simulations to study the effect of surface modification on the interactions between the CNCs and the resulting bundling process. We consider two types of common surface modifications of native CNCs, sulfated CNCs (SCNCs) and TEMPO-oxidized CNCs (TCNCs), in the presence of two types of counterions, Na+ and Ca2+, in solution. We used the umbrella sampling method to calculate the potential of the mean force (PMF), and we found that the strength of interaction between the modified CNCs decreases, compared with the native CNCs. The strength of interaction for TCNCs is almost similar to that for SCNCs at the same level of surface substitution, whereas the type of counterion has a strong effect on the PMF with a higher interaction energy between the CNCs in the presence of a divalent counterion as compared to a monovalent counterion. Finally, we studied the self-assembly of CNCs into a hexagonal bundle for the native CNCs and sulfated CNCs focusing on the twist of the bundle, bound water inside the bundle, inter-CNC gap, and interaction energy between the CNCs in the bundle, and the effect of the counterions on the morphology of the bundle. The equilibrium spacing of the CNCs within the bundle is found to be consistent with the results of PMF calculations for the minimum separation distance between the respective crystal surfaces.