Crumbling Crystals: On the Dissolution Mechanism of NaCl in Water

Life on Earth depends upon the dissolution of ionic salts in water, particularly NaCl. However, an atomistic scale understanding of the process remains elusive. Simulations lend themselves conveniently to studying dissolution since they provide the spatio-temporal resolution that can be difficult to obtain experimentally. Nevertheless, the complexity of various inter- and intra-molecular interactions require careful treatment and long time scale simulations, both of which are typically hindered by computational expense. Here, we use advances in machine learning potential methodology to resolve for the first time at an ab initio level of theory the dissolution mechanism of NaCl in water. The picture that emerges is that of a steady ion-wise unwrapping of the crystal preceding its rapid disintegration, reminiscent of crumbling. The onset of crumbling can be explained by a strong increase in the ratio of the surface to volume of the crystal. Overall, dissolution is comprised of a series of highly dynamical microscopic sub-processes, resulting in an inherently stochastic mechanism. These atomistic level insights now pave the way for a general understanding of dissolution mechanisms in other crystals, and the methodology is primed for more complex systems of recent interest such as water/salt interfaces under flow and salt crystals under confinement.


I. INTRODUCTION
Understanding the dissolution of crystals is vital for a myriad of pressing modern day challenges, from technological issues such as battery science [1] and water desalination [2] to drug bioavailability [3] and geochemical weathering [4]. The macroscopic process of dissolution is well described and understood from thermodynamics, where enthalpies of dissolution for example can be readily obtained from experimental techniques such as calorimetry [5]. However despite well understood general theories describing the conceptually reverse process of nucleation [6], at a microscopic level our understanding of dissolution remains largely unresolved.
Dissolution is indeed an old problem [7], however with recent experimental developments, it is an exciting time for the field. Building on early atomic force microscopy studies of the water/NaCl interface by Xu et al. [8], sophisticated imaging techniques such as scanning tunneling microscopy [9] and Bragg coherent diffraction imaging [10] have given insight into low temperature surface ion dissolution events and the role of defects in the dissolution process. In addition, recent impressive highresolution liquid cell transmission electron microscopy work has observed defect-mediated ripening of Cd-CdCl2 core-shell nanoparticles [11]. The solid-liquid interface under flow has also been probed with a combination of surface-specific sum frequency generation spectroscopy * cs2121@cam.ac.uk † am452@cam.ac.uk and microfluidics experiments, revealing drastic changes to the equilibrium of dissolved ions. [12] Recent singlemolecule atomic-resolution real-time electron microscopy experiments have captured in real time a NaCl nucleus emerging during nucleation [13] while atomic electron tomography has been used to monitor the nucleation of FePt nanoparticles over time [14]. However in general, obtaining atomic scale dynamical information from experiment remains a highly challenging endeavour.
In contrast, computer simulations conveniently grant access to theÅngstrom and femtosecond scale resolution required, and have been applied to study the dissolution of a range of materials (see e.g. refs [15][16][17][18][19][20][21][22][23]. Important insights have been obtained from these simulations that have significantly advanced the field, such as the fact that dissolution rates are highly controlled by the crystal structure, and that the existence of site specific barriers to dissolution implies a non-constant dissolution rate throughout the process [24][25][26][27][28]. Nevertheless, simulation studies up to now have had limitations. Notably, both force field (FF) and ab initio studies have been limited to individual trajectories of dissolution events. However, a stochastic dissolution process is expected, as shown for example in larger scale studies in the geosciences [29]. Therefore it is important to establish if the sequence of dissolving atoms can be rationalised at an atomistic level from an ensemble of dissolution trajectories. In addition, nucleation pathways are strongly dependent on solution concentration [30][31][32], raising the question whether dissolution follows the same predictable mechanism or if it is highly sensitive to the conditions. putationally efficient explorations of multiple dissolution trajectories using a methodology that accurately describes the delicate changes in water-water and water-ion interactions that occur during dissolution. Fortunately, with the emergence of machine learning interatomic potentials (MLPs) [33][34][35][36], accurate and efficient potentials can be developed for the treatment of complex processes such as dissolution. Such potentials yield the accuracy of ab initio methods but at a fraction of the computational cost. Of the many flavours of MLP methodology, the approach reported in ref. 37 is particularly suitable as it enables the automated development and validation of MLPs capable of accurately treating complex aqueous systems.
Given the opportunities presented by the recent developments in MLP technology, we have performed a detailed study of NaCl dissolution. NaCl dissolution is a prototypical and widely studied system [38][39][40][41][42][43][44][45][46][47][48], because of its ubiquity and significance to phenomena including biological intracellular reactions [49] and climate chemistry [50,51], not forgetting that two-thirds of the Earth is covered in salty water. We find that the dissolution process is highly dynamic. Many stochastic subprocesses combine to give an overall crumbling mechanism, in which a steady period of ion-wise dissolution precedes the rapid concerted disintegration of the crystal. This disintegration is governed by the steep increase of the surface to volume ratio of the crystal, corresponding to an unfavourable surface to bulk free energy -borrowing from ideas of classical nucleation theory (CNT)that ultimately leads to the rapid collapse of the crystal. Figure 1 provides an overview of a prototypical system considered and the NaCl dissolution process. Throughout this discussion, we show representative trajectories at one temperature but the phenomena observed hold for multiple trajectories at 330 K and 400 K. Full details of the systems and simulation protocols are given in the Methods section. However, in brief, we considered 64 and 216 ion NaCl nanocrystals in simulation boxes containing 1250 to 625 water molecules (giving NaCl concentrations when dissolved of 1.42, 2.84 and 5.61 mol/kg). While we have not computed the solubility of the ML model, the experimental saturation of NaCl in water is 6.15 mol/kg. All systems explored started out with a NaCl nanocrystal in pure water and the final dissolved state is likely undersaturated or at worst at a low supersaturation, thereby facilitating the complete dissolution of the crystal. A typical dissolution trajectory is shown in Figure 1 (a), with the extent of dissolution of the crystal monitored over time. Relevant snapshots from the trajectory are shown in the upper panel. After a steady period corresponding to ions dissolving from typically low coordinated sites of the crystal into solution, there is a rapid increase in the rate of dissolution. This suggests the crystal reaches some critical point of instability after which it rapidly disintegrates or crumbles. Following this, the Cland Na + ions are fully solvated and the system has lost all crystalline order. These general observations for the dissolution mechanisms are consistent over all conditions studied and reproducible over the ensemble of trajecto- ries at a given condition.

II. THE NACL DISSOLUTION MECHANISM
The early stages of dissolution up to the onset of crumbling comprise a set of dynamic processes occurring in equilibrium as shown in Figure 1 (b). The dissolution of any ion from the crystal is first preceded by the ion rolling onto the surface of the crystal (Structures I. & II.). This process involves a simultaneous reduction in coordination of the rolling ion with its neighbours in the lattice until it is just coordinated with one counter ion, and a corresponding increase in coordination of the ion with water. This behaviour is consistent with earlier AIMD studies on NaCl and Li x Mn 2 O 4 (001) surfaces. [40,52] Moreover, there is a difference between the rolling behaviour of Na + and Cl -: Na remains much more stable on the surface, while Clis much quicker to dissolve into solution. This is in agreement with recent work by Silvastri et al. [15], in which they show that Na + ions are stable on terrace sites when they move from kink sites, while Clhas much weaker minima on terrace sites. There are then several possibilities for the fate of this exposed ion: It can become further solvated by water and completely dissolve from the crystal (Structure III.) It can also roll back to its starting position and/or further roll to an orthogonal face (where it again has the possibility to dissolve from the crystal) (Structure IV.). The dissolved ion can also rejoin the crystal (at not necessarily the same location as it originally emerged). In these early stages before the rapid collapse of the crystal, other ions can also simultaneously roll and move about the crystal without necessarily any dissolution. This results in a deformation of the crystal shown in Structure V. The implications of these early dynamic processes with respect to the stochastic nature of dissolution will be discussed in detail in Section III. Figure 2 (a) shows a schematic of the typical 'unwrapping' of the crystal. The ions are coloured according to their time taken to dissolve. Overall there is a preference for ions from low coordinated sites to dissolve first. We observe that dissolution is initiated at corner sites in the lattice, which has been previously noted in several papers [38,40,41,47,53]. The dissolution of ions up to the point of crumbling proceeds in an approximately stepwise manner of ions of opposite charge (Cl -, Na + , Cl -...), thereby minimising charge accumulation on the crystal lattice, confirming earlier force-field simulations on the initial stages of dissolution [47]. At each dissolution step (depending on its location) for an ion to dissolve it must break 3 (corner), 4 (edge), 5 (face) or 6 (centre) ionic bonds. Therefore the edge-wise unwrapping of the crystal shown in Figure 2 (a) follows the hypothesis that at each step, the number of ionic bonds broken is minimised. This holds true up to the point of rapid disintegration, where all atoms irrespective of their current coordination dissolve within a few picoseconds.
These initial observations raise numerous questions re-garding details of the crumbling mechanism proposed. Specifically: (i) Does the crystal splinter into smaller crystalites or completely crumble into ions? (ii) Nucleation studies of NaCl in water, have shown that at very high concentrations there is an amorphous ionic structure before ions from solution form a nucleus [30,31]. Therefore is there also evidence for an amorphous intermediate before the complete dissolution of the crystal? (iii) At what point does the crystal disintegrate, and is this crumbling event determined by the size of the remaining crystalline cluster or something else?
To answer these questions, we now analyse the nature of the crumbling mechanism, and attempt to quantify and describe the driving force for this event. Let us first address the question whether the remaining crystal splinters into smaller subclusters which is suggested in Ref. 41. Another possibility is that it rapidly crumbles into individual ions. Monitoring the sizes of the 4 largest crystalline (ionic) clusters and the crystalline order of the largest cluster over time, as shown in Figure 2 (b) gives insight into the exact mechanism of the crystal collapse. Considering the sizes of the 4 largest ionic clusters, initially there is just one large cluster with 64 ions (i.e. the initial 4 × 4 × 4 crystal), which gradually decreases in size corresponding to the steady period of ions dissolving. At the point of crumbling, the number of ions in the largest cluster rapidly drops and immediately converges to a steady value of just one or two ions rapidly fluctuating over time (corresponding to short-lived ion pairs in solution). Since at no stage is there a second cluster with some intermediate number of ions, this implies that there is a complete disintegration of the crystal, rather than something resembling a shattering process.
To understand if an amorphous transition precedes crumbling, the crystalline order of the ions in the system has been monitored over time.q 8 (defined in the SI) is a variation of the typical Steinhardt bond order parameter q 8 , and averages the bond order vectors over the first shell of neighbouring ions to provide a measure of ordering of the ions in the system.q 8 (Figure 2 (b)) qualitatively tracks the evolution of the largest cluster size, whereby initially there is a steady decrease (but still within the range expected for a crystal) corresponding to the crystal decreasing in size. At a critical value ofq 8 ≈ 0.3, there is a sharp decrease inq 8 , corresponding to a rapid total loss of crystalline order in the system. Therefore for the (low) concentration regimes we explore, our results are in agreement with the nucleation literature, where there is a single step order to disorder transition, with no amorphous intermediate.
We have already suggested that the onset of crumbling occurs when the crystal reaches a critical point of instability. Across the simulated trajectories, there is a wide spread of crystal sizes at which the disintegration occurs (approximately 20-40 ions). So it is not size alone that governs the instability of the nanocrystal. To understand this behaviour we considered the surface area to volume ratio of the convex hull of the largest cluster. Panel (c) in Figure 2 shows this surface area to volume ratio for the crystal over time up to the point of collapse. This gives a measure of the extent of interaction of the ions in the crystal with water over time. There is an initial decrease corresponding to a 'rounding' of the crystal, to a relatively constant value up to the point of crumbling. (This rounding is even more pronounced for the larger 6 × 6 × 6 crystal, Figure S10.) At the onset of crumbling there is then a rapid increase in the surface area accessible to the water. Our system is not in a concentration regime where reversible dissolution and nucleation is to be expected. However, with this ratio we can draw direct parallels with CNT. Within CNT, the free energy of formation of a crystal nucleus comprises a volume term representing the gain in free energy upon formation of a stable crystalline phase, and an unfavourable surface term arising from the cost of forming an interface. The dominance of the surface area over the volume at the point of crumbling of the crystal lattice in this dissolution work is akin to a rapid increase in the (destabilising) surface free energy in CNT. Therefore there is a delicate balance of intra-and inter-molecular interactions that ultimately determine the crystal stability.
We note that the cluster at the point of crumbling here is not equivalent to a critical nucleus in CNT. For the concentration regimes explored here, the critical nucleus would be expected to be much larger [30]. To further test the robustness of our proposed crumbling mechanism, we performed simulations on the same systems using the Joung Cheatham force field for ion-ion interactions with the SPC/E water model [54]. This has a computed solubility of 3.7 mol/kg, and therefore we can be sure that our systems are below saturation when fully dissolved for this model. Full details and results are given in the SI, however we observe qualitatively the same crumbling mechanism, with similar cluster sizes at the point of crumbling. Therefore despite not having the solubility of the ML model, we can be confident that this mechanism is not simply a reverse of CNT.
In summary, this ratio is a simple parameter that explains the onset of disintegration of the crystal. It can also be used directly for other systems, thereby allowing for a simple and intuitive understanding of dissolution processes in general.

III. STOCHASTIC NATURE OF THE NACL DISSOLUTION PROCESS
Thus far, we have described the overall crumbling mechanism of NaCl dissolution in detail, which is completely general for multiple conditions of concentration and temperature. However a pertinent question relevant for many fields is whether the dissolution mechanism is deterministic. That is, given some initial conditions, can we predict at what crystal size/structure and after how long will the crystal dissolve? Previous ab initio studies [40,41] and larger scale FF-based studies on disso- lution have limited their focus to individual trajectories, thus restricting the conclusions that can be drawn regarding the mechanism in this respect. However in this work, access to multiple ab initio quality trajectories via machine learning based simulations has enabled this issue to be addressed through a thorough statistical analysis over multiple trajectories. Figure 3 shows the time evolution of the crystal size distribution for 10 trajectories of a 6 × 6 × 6 nanocrystal. The size distribution is initially very narrow at the beginning of the dissolution. However as dissolution proceeds along the 10 trajectories, the distribution widens significantly. After approximately 15 ns, there is a spread of over 100 ions in the crystal sizes along different trajectories. While all trajectories follow the overall crumbling mechanism described in Section II, this broadening of the distributions illustrates the diverse microscopic paths taken by the crystals en route to dissolution. We previously discussed in Section II the highly dynamic equilibrium of microscopic subprocesses such as ions rolling on the crystal surface and lattice rearrangement occurring during the steady period of ion dissolution from the crystal. These are rooted in inherent stochastic system fluctuations such as lattice and molecular vibrations (on the order of fs) and water dipole and hydrogen bond reorientations (on the order of ps). Therefore the wide variation in the crystal structures at a given time (despite all simulations starting from the same pristine crystal structure) can be directly attributed to the stochastic nature of these microscopic sub-processes. Given the relatively small sizes of the systems studied here, one would expect a much larger spread of times in macroscopic systems. In addition to the range of crystal sizes for different trajectories at a given time, for a given crystal size, there are also numerous possible structures, as shown in the inset of Figure 3, with 3 representative examples of a crystal with 87 atoms.
These diverse structures again arise from a combination of the stochastic microscopic subprocesses. The diversity in microscopic mechanism can also be observed for the other concentration conditions studied here. The largest 6 × 6 × 6 crystal (5.61 mol/kg) takes much longer to dissolve and has a much broader range of dissolution times than the 4 × 4 × 4 crystals (1.42 mol/kg and 2.84 mol/kg), as shown in Figure S10. This broad spread of times can be attributed to an increase in probability of ions rejoining the 6 × 6 × 6 crystal, since it is surrounded by the lowest volume of water, but also there is a cumulative effect of each possible stochastic sub-process for every dissolution step. While a distribution of dissolution times is to be expected, a 2 or 3-fold spread in times is particularly noteworthy. In summary, the dissolution of NaCl therefore proceeds via an overall crumbling mechanism, within which, there is rich structural variety of the dissolving crystal owing to the inherently stochastic microscopic system fluctuations.

IV. CONCLUSION
Despite dissolution being one of the most ubiquitous processes on earth, major gaps in understanding at the atomistic level persist. In this work we provide a rigorous description of the mechanism of NaCl dissolution, through the application of an automated framework to develop accurate MLPs. Surpassing the extent of previous simulation work we have generated multiple trajectories (> 300 ns in total) with ab initio equivalent accuracy, under a range of conditions of NaCl dissolving. We have established a general crumbling mechanism of NaCl dissolution, where following a steady period of ion-by-ion loss, the crystal reaches a critical point of instability and collapses in a concerted fashion. The sequence of dissolving ions in the steady period of dissolution is such that the number of ionic bonds broken and net charge on the crystalline lattice is minimised. The overall stability of a crystalline cluster in water can be straightforwardly reasoned by a delicate balance of interactions between the ions and water, measured via the surface area to volume ratio of the crystal. This simple and physically intuitive concept applies to NaCl and to other ionic crystals in water and can readily be extended to non-ionic systems such as molecular crystals. The crumbling mechanistic insights we have provided are primed for experimental measurement, including confinement controlled monitoring of the dissolution and Bragg coherant diffraction imaging, as have already been successfully applied to crystalline systems [12,13].
We have shown that the overall nature of the dissolution mechanism is highly stochastic, comprising a dynamic equilibrium of sub-processes such as ions rolling, dissolving, rejoining and deformations of the lattice in the crystal. These arise from the stochastic hydrogen bond formation/breaking and thermal vibrations in the crystal and influence the relative stability of the crystal, and therefore the overall time taken for dissolution.
Given the sizes of the nanocrystals studied here, one could envisage a much larger spread of times in macroscopic systems such as those considered in the geosciences. Indeed, studies of dissolution of crystalline and non-crystalline materials in the earth sciences have observed a broad variation in the rate constant of dissolution of up to 2-3 orders of magnitude. [29,55] which they attributed to the inhomogeneity of the initial crystal surface. However our results suggest that this variance in rate is inherent to the system, regardless of initial condition.
Such observations regarding the dynamic nature of the crystal are indeed very amenable to further experimental studies. They suggest a large variation in dissolution mechanism under flow compared to confinement. In the former case the ions are rapidly removed from the surface -similar to the lower ion concentration conditions we have shown here, while confinement essentially traps the ions, creating a high concentration of ions in the vicinity of the surface. We propose that revisiting the microfluidics [12] and electron microscopy [13] experiments described in the Introduction, with a focus on different concentration conditions would be highly insightful.
The accurate treatment of the electronic structure of the system required to capture relevant interactions, yet at a computational cost that facilitates the simulation of many dissolution trajectories is crucial for the observations made here. The high transferability of the methodology and validation process is amenable to the study of other ionic salt systems. Indeed, having now understood NaCl dissolution -in terms of a minimization of charge and of bonds broken, a balance of ion-solvent interactions determining crystal stability, and the dynamic nature of the dissolution process -our work sets the stage for a generalised theory of dissolution of ionic crystals. Computing the solubility of the ML model is also an important next step, given the large differences from experiment in solubilities predicted by different force field models [54]. Finally, this methodology and the insights we have now obtained are primed for application to more challenging and complex systems of recent interest, including highly concentrated electrolyte solutions [56,57] and electrolyte solutions under confinement [58].

V. METHODS
Machine Learning Potential: Machine learning potentials provide a direct functional relationship between atomic positions and forces/potential energy. This bypasses the otherwise computationally expensive requirement to solve Schrödinger-like equations and systems of much larger length and timescales than feasible with AIMD based approaches are now accessible. Building upon the seminal work of Behler and Parrinello in neural network based ML potentials [59], the recent development of an active learning approach using a committee of NNPs has enabled the systematic development and validation of NNP representations of the potential energy surface for complex aqueous systems [37,60] Our model was iteratively trained on forces and energies from carefully chosen ab initio molecular dynamics simulations in the generalised gradient approximation using the rev-PBE functional with Grimme's D3 dispersion correction [61], while using a Coulomb baseline to incorporate the correct long-range electrostatics, evaluated with particle mesh Ewald summation. This baseline uses TIP3P model parameters [62] for water and point charges of +/-1 for Na and Cl. This setup has been shown to be accurate for aqueous systems, and further details on functional selection and long-range electrostatics are described in the SI.
The model was iteratively improved over multiple generations such that the relevant configuration space was accurately covered in the overall training set with an energy and force root mean square error (RMSE) of 1.3 meV per atom and 38.0 meV/Å, respectively. The model was extensively validated in a series of tests incorporating static and dynamical properties of selected test systems of which full details are also given in the SI.
MD simulations: The systems used in all dissolution simulations described comprised a 4×4×4 (6×6×6) NaCl nanocrystal with 32 (108) Na/Cl atoms using a computed lattice constant of 5.72Å [63]. The 4 × 4 × 4 nanocrystal was surrounded by 625 (1250) water molecules in a cubic simulation cell of side length 27.8Å (34.1Å), giving dissolved NaCl concentrations in water of 2.84 mol/kg and 1.43 mol/kg respectively. The 6 × 6 × 6 crystal was surrounded by 1068 waters in a cubic simulation cell of side length 34.1Å giving a dissolved NaCl concentration in water of 5.61 mol/kg. MLP simulations were all carried out using the CP2K/Quickstep code in the NVT ensemble and at a constant temperature of 330 K (chosen to be consistent with previous literature [64]) maintained using the CSVR thermostat [65]. At each concentration, 10 trajectories were simulated using the same initial configuration and initialised with random velocities drawn from the Boltzmann distribution at the target temperature of 330 K. Simulations were run for over 15 ns (5.61 mol/kg), 2 ns (2.84 mol/kg) and 1 ns (1.42 mol/kg) until the crystal was fully dissolved. The same simulation procedure was carried out for 10 trajectories of the 2.42 mol/kg concentration at 400 K for further validation of the generalisability of our conclusions. These results are included in the SI. Overall over 300 ns of ab initio quality machine learning based simulations were performed, far beyond the capabilities of AIMD simulations.
Further details for computing parameters including cluster size,q 8 and surface area to volume ratio of the crystal are given in the SI.

VI. ACKNOWLEDGMENTS
We are grateful to Daan Frenkel for helpful discussions. We thank Julian Gale for his insightful comments on the manuscript. N.O.N acknowledges financial support from the Gates Cambridge Trust. C.S acknowledges financial support from the Alexander von Humboldt Stiftung. S.J.C is a Royal Society University Research Fellow (URF\R 1\211144) at the University of Cambridge. We are grateful for computational support and resources from the UK Materials and Molecular Modeling Hub which is partially funded by EPSRC (Grant Nos. EP/P020194/1 and EP/T022213/1). We are also grateful for computational support and resources from the UK national high-performance computing service, Advanced Research Computing High End Resource (ARCHER2 Validation of Dissolution Process S10 Long-range Effects S12

S1. SELECTING DFT FUNCTIONAL
Two modern DFT functionals that have been shown to perform well for liquid water simulations S1,S2 -the dispersion corrected GGA functional revPBE-D3 S3,S4 and the strongly constrained and appropriately normed (SCAN) meta-GGA functional S5 -were selected for S2 an initial benchmark study to determine the most suitable XC functional for the ab initio simulations upon which the MLP would be trained. It should be noted that dispersion corrected hybrid functionals have been shown to provide slightly better results for liquid water simulations; S6,S7 however initial simulations with revPBE0-D3 showed there is approximately a six-fold increase in computational cost, rendering longer time-scale AIMD simulations unfeasible for the envisaged system sizes.

Na/Cl ion in water
To study the performance of SCAN and revPBE-D3 for ion-water simulations, four separate AIMD simulations were performed for both revPBE-D3 and SCAN for Na + in water and Clin water. The Na(Cl)-O RDF g(r) are compared to experimental values in Figure   S1 and Table I. Both SCAN and revPBE-D3 are in good agreement with the position of the first experimental Cl-O peak (R 1 max ), where both are within 0.07Å of the experimental data. However the height g(R 1 max ) is underestimated by both by approximately 15%. SCAN and revPBE-D3 both qualitatively agree with the position of the second peak maximum, but they both tend to overstructure compared to experiment. The experimental Na-O peak is more accurately reproduced by SCAN than revPBE-D3. The revPBE-D3 peak is shifted 0.17Å from the experimental value, and the peak height is 27% lower. This is also in agreement with previous studies, where SCAN has been shown to outperform GGA functionals in reproducing the sodium-water solvation structure. S8 Despite a lack of experimental data for the second solvation shell, revPBE-D3 is in much closer qualitative agreement with SCAN than for the first peak. The O-O RDF for both solvated sodium and chloride ion systems was also computed and compared to experiment in order to evaluate the performance of the two functionals for water. This analysis is shown below in Figure S2. RDF. As a footnote, SCAN is also two-thirds more expensive than revPBE-D3. Therefore, the AIMD trajectories generated using revPBE-D3 were chosen for the initial development of the C-NNP training set and revPBE-D3 was used for generation of subsequent AIMD data and during the active learning stages of the MLP development.

Automated Work Flow
The procedure for developing the committee neural network potential (C-NNP) was fol-

Details of Model
The chemical environment around each atom was described using a general set of atomcentered symmetry functions. S15 There are 10 radial and 4 angular functions for each pair and triple of atoms, following Ref. S12. All symmetry functions used a cutoff function of angular cosine form with a cutoff radius of 12 Bohr. The committee was comprised of 8 NNP members, of identical architecure with 2 hidden layers and 25 neurons in each layer. In all cases, random sub-sampling was performed to introduce variability between the committee members, where 10% of the total set of structures were discarded. The weights and biases of the NNPs were optimised using the n2p2 code. S16 Individual models during active learning were optimised for 15 epochs, while the final C-NNP model used in simulations was optimised for 50 epochs.
To ensure that long-range effects are accurately captured by the machine learning potential, we explicitly incorporated long-range effects beyond the cutoff of the symmetry functions (12 Bohr). The predicted energy can in general be written as a sum of short range and long range contributions (E sr and E coul respectively): E tot = E sr + E coul . The longrange model was thus trained on the difference between the standard short-ranged model S6 and the Coulomb contribution, calculated using point charges of +/-1 respectively for Na and Cl and using TIP3P model parameters for water. S17 We used this model in all production simulations, where the Coulomb contributions were explicitly included via particle mesh ewald summation. Details on the validation of the final model and the role of long range interactions in these simulations are presented in the next Section. The number of configurations included from each AIMD simulation is given in Generation-01. Thereafter, the number of configurations is the total number of structures in the training set for the given generation.

Validation
The system studied is complex, so care was taken to set up a series of validation tests that incorporates relevant subsystems and properties.
2 × 2 × 2 NaCl nanocrystal The static and dynamical properties of the largest accessible NaCl nanocrystal in water for AIMD simulations (2×2×2 NaCl surrounded by 229 waters -695 atoms) was analysed using a benchmarking scheme proposed in Ref. S12, comparing AIMD to the C-NNP predictions and shown in Figure S4.

Long-range Effects
In order to evaluate the effect of explicitly incorporating long range effects into the model structural properties (RDFs) of the 2.48 mol/kg NaCl solution were computed. These are compared to the short range model in Figure S6. Both predictions are in very good agreement for all cases of water, ions and ion-water RDFs. Nevertheless small discrepancies between long-and short-range model predictions in the cases of the Cl-Cl and Na-Cl RDFs prompt us to use the long-range model for production runs, given the additional physical information included in the model. In general however this analysis suggests that overall for our simulated systems, long-range effects do not have a substantial impact. S12

S13
In summary this set of validation tests confirm the overall quality and accuracy of the C-NNP model. The model has performed comparably to or exceeded literature values on a variety of systems and properties that are very relevant for the dissolution process, from small scale studies of single ions in solution to systems sizes and interactions that will be explored in the detailed analysis of the dissolution process.

System setups
All simulations were carried out in an orthorhombic simulation cell with periodic boundary conditions in x, y, z directions. Details of the simulation setups for all systems are shown in Table ??.

Simulation details
All MD/AIMD simulations were performed using the CP2K/Quickstep code. S19 AIMD simulations were used to generate reference data for development of the model as shown in Table III, with pre-equilibrated structures obtained from force field simulation. C-NNP  production simulations as shown in Table IV were initially equillibrated with fixed ions using the C-NNP model.

Force field pre-equilibration details
A 1 ns force-field based MD simulation was first performed on all Na/Cl containing systems from Table III to obtain a pre-equilibrated initial configuration for subsequent AIMD simulations. The Born-Mayer-Huggins potential for NaCl S20 and TIP3P model for water S17 were used and the water/NaCl interactions were described by the parameters given by

AIMD simulation details
Prior to all AIMD simulations, a short constrained AIMD simulation, keeping ions fixed, was performed on all the systems given in Table III containing ions, to obtain a wellequillibrated water structure around the Na and Cl ions. These were used as starting configurations for simulations in NVT ensemble within the generalised gradient approximation using the revPBE functional S4 with Grimmes' dispersion correction. S23 The electronic density was partitioned into core and valence contributions, with core electrons described using the norm-conserving Goedecker, Teter and Hutter (GTH) pseudopotentials. S24 Na 2s and 2p electrons were also treated explicitly given the well-known issue of non-linear core-valence exchange/correlation. Valence electrons were described using the MOLOPT TZV2P basis set. S25 Stochastic errors for Na forces due to insufficient basis set convergence were observed in additional DFT convergence tests for the cutoff of the plane-wave representation of the charge density. Therefore a plane wave cutoff of 1200 Ry was used to obtain accurate forces.
Deuterium masses were used in lieu of hydrogen, allowing a 1 fs timestep and reducing errors due to neglecting nuclear quantum effects through classical propagation of the nuclei. Simulations were performed at 300 K, maintained using the CSVR thermostat.. S26 All AIMD simulations were at least 40 ps long (except for 2 × 2 × 2 NaCl nanocrystal which was 25 ps due to computational cost.) and can be found in Table III.

Steinhardt order parameters
The Steinhardt order parameter q8 is a particular case of a larger set of bond order parameters based on spherical harmonics. S27 To distinguish between solid and liquid-like ions we first compute the complex vector q 8 (i) for each particle i, where N b (i) is the number of nearest neighbours of particle i, and m runs from -8 to +8.
S17 Y 8,m are the spherical harmonics and r i,j connects particles i and j. Following Lechner et al. S28 we then compute an averaged form of this local bond order parameter by averaging the local q 8,m (i) vectors over the particle i and its surroundings: Computation of these bond order parameters was performed using the Freud S29 python package.

Coordination Numbers and clusters
Computation of all coordination numbers and clusters were performed using Plumed, S30 and analysis carried out using Python. The CN of an atom i with a set of atoms in A is S18 defined by: where s ij is a smooth switching function with a range of [0,1]: r ij is the interatomic distance between atom i and atom j. r 0 for Na-O (2.80Å) and Cl-H using a threshold of a ion-ion CN > 2, and then using these atoms as input for a depth first search graph reduction algorithm to identify the sets of ions clustered together.
To compute the surface area to volume ratio of the dissolving crystal, a DFS algorithm was implemented using the Freud Python package S29 in order to identify the ions in the S19 largest cluster at a given time. The convex hull of this set of points was then obtained (provided there were > 3 ions in the cluster), and the surface area to volume ratio of the hull was computed.

MLP
Results from identical simulations as described above for the C-NNP at 2.48 mol/kg NaCl concentration at 400 K are shown in Figure S9. The additional results for the 1.42 mol/kg and 2.48 mol/kg trajectories at 330 K are also shown in Figure S9.

FF
Results from identical simulations as described above for the C-NNP at 1.42 mol/kg and 2.48 mol/kg at 330 K are shown in Figure S10 for a JC/SCP/E force-field model.