Molecular dynamics (MD) simulations of 13C NMR longitudinal relaxation (T1) distributions were recently established as a powerful tool for characterizing moisture adsorption in natural amorphous polymers. Here, such computational-experimental synergy is demonstrated in a system with intrinsically high structural heterogeneity, namely crystalline cellulose nanofibrils (CNFs) in highly hydrated aggregated state. In such a system, structure-function properties on the nanoscale remain largely uncovered by experimental means alone. In this work, broadly polydispersed experimental 13C NMR T1 distributions could be successfully reproduced in simulations and, for the first time, were decomposed into contributions from distinct molecular sources within the aggregated CNFs, namely, (i) the core and (ii) the less-accessible and accessible surface regions of the CNFs. Furthermore, within the surface groups structurally different sites such as (iii) residues with different hydroxymethyl orientations and (iv) center and origin chains could be discerned based on their distinct molecular dynamics. The MD simulations unravel a direct correlation between dynamical and structural heterogeneity at an atomistic-level resolution that cannot be accessed by NMR experiments. The proposed approach holds the potential to enable quantitative interpretation of NMR data from a range of multicomponent high-performance nanocomposites with significantly heterogeneous macromolecular structure.