Natural organic matter (NOM) and humic substances (HS) as its part are key players in diverse biogeochemical processes occurring at interfaces at different scales, contributing to the carbon cycling, sorption reactions, transport and environmental fate of chemicals and biological activity. One important vision considers HS composed of small molecules interacting through non-covalent forces rather than representing high-molecular-weight substances. There is a continued interest in understanding the structure of NOM/HS and modeling their physico-chemical characteristics using molecular simulations. This work aimed to examine the stability of HS models, represented as aggregates formed by small molecules, upon dilution in water, by performing molecular dynamics (MD) simulations. The HS models were built to reproduce a 13C NMR-derived composition of standard Leonardite humic acid (HA) by using the Vienna Soil Organic Matter Modeler. The models varied in molecular size, amount of water and the choice of counter cation (Na+ vs Ca2+), while keeping the chemical composition identical, to link low-hydrated HS systems with aqueous NOM solutions in the concentration range of environmental significance. MD simulations demonstrated that HS aggregates disintegrate to some degree by diluting the system. The process occurs more easily in the presence of monovalent Na+ compared to divalent Ca2+. Strong interactions within the HS aggregates may remain at the level of dimers of small molecules revealing a higher stability against dilution. The simulation results support a step-wise mode of HS aggregation. It is suggested also that the existence of large-size HS molecules possibly undergoing step-wise aggregation should be taken into account while examining HS properties in aqueous solutions.