MASKE: A kinetic simulator of coupled chemical and mechanical processes driving microstructural evolution

A typical simulation features spherical particles in a prismatic simulation box with volume $V_{box}$ and periodic or fixed boundary conditions: see Fig. 1. Box and particles are generated and tracked using the LAMMPS C++ library interface (Thompson et al., 2022). The particles may represent any explicit phase, hereafter assumed to be solid ones. The particles interact mechanically with each other through user-defined interaction potentials, which are also defined and managed in LAMMPS: e.g. the $U_{ij}(r_{ij})$ interaction in Fig. 1. The current version of MASKE can only manage pairwise interactions, but extensions to any of the potentials in LAMMPS would be quite straightforward: see A for more discussion on this. Particle displacements caused by the mechanical interactions are also computed in LAMMPS, via usual methods such as molecular dynamics or energy minimization (respectively using the $run$ and $minimize$ commands in LAMMPS). This interface with LAMMPS gives the user flexibility to study complex systems with various mechanical interaction potentials and loading conditions.

The portion of a simulation box that is not occupied by solid particles contains a volume $voidV$ of empty space, and a volume $V_{sol}$ of implicit solution made up of various, user-defined molecular species. The current version of MASKE assumes that the molar concentrations of the solvated species are uniform in $V_{sol}$, with no diffusion taking place. There is also an additional volume $\Delta V_{box}$ attached to the simulation box, which can receive some of the molecules released during dissolution reactions, or provide some of the molecules consumed during precipitation reactions. This $\Delta V_{box}$ can be used to keep $V_{box}$ limited, which is where time-consuming operations take place; for example, Fig. 2 depicts a scenario where this approach allows focussing on dissolution and precipitation near the surface of a large grain, while preserving realistic concentrations in solution and surface-to-volume ratio of suspended grains at a larger scale. The initial molar concentrations of all the solvated molecules in $V_{box}$ and in $\Delta V_{box}$ are provided by the user in the input script, along with fixing the temperature of the solution and the $A$ and $B$ solvent parameters to be used in Debye-Hückel calculations of activity coefficients.

MASKE simulates discrete events, which happen once every $\Delta t_{k}$ units of time (where $\Delta t_{k}$ varies during a simulation depending on mechanical interactions and solution composition), but also continuous processes that are numerically integrated over a user-specified time interval $\Delta t_{i}$ with time step $dt_{i}$ (for the $i^{th}$ continuous process). In the current version of MASKE, the possible discrete events are full dissolution (called $deletion$) of any existing particle, or full precipitation (called $nucleation$) of a new particle at any possible location in $V_{box}$. The only continuous process implemented so far is called $CmstoreLMP$: it repeatedly executes a set of stored LAMMPS commands, specified by the user in the MASKE input script (see C for more details on the input/output structure of the code). A $CmstoreLMP$ process will be used in an example here, to impose a constant strain rate. The occurrence of discrete events and continuous processes is staggered onto a unique time line, $t^{*}$, which is advanced by whichever event or process has the shortest time of future occurrence, or $FOT$. The details of how this unique time line is managed can be found in B. The occurrence of discrete events is sampled using Kinetic Monte Carlo (KMC), thus the list of all possible discrete events is associated with a unique $FOT_{k}$, which is the future occurrence time of any one out of all the possible discrete events. If $FOT_{k}$ is smaller than $FOT_{i}$ of all the other continuous processes, i.e. if a discrete event is next in line to be executed, then one discrete event is picked out from the list of all the possible ones, with a probability that is proportional to its individual rate. This is the usual approach in KMC simulations; specifically, MASKE uses a rejection-free KMC algorithm that is time-dependent (Prados et al., 1997), to account for changes in solution composition or stress state that continuous processes may induce between two successive discrete events. The equations of time-dependent KMC are also presented and discussed in B.

In the KMC method, $FOT_{k}$ is obtained from the cumulative rate:

which is the sum of the individual rates $r_{k,i}$ of all $N$ possible particle deletions plus all $M$ particle nucleations at all possible locations in the simulation box. The expressions for $r_{k,j}$ are presented in the next section. $N$ is finite, whereas the number of possible locations for particle precipitation is infinite and, in principle, one should sample them all to use KMC. To treat this problem, MASKE asks the user to select a region within $V_{box}$ and to create a 3D lattice of $M$ trial particles of the desired type, viz. with desired chemical composition and interaction potentials. Each of these $M$ trial particles is then driven towards its closest local minimum of interaction energy with the other already-existing, so-called $real$ particles (no interactions between trial particles are considered, and the positions of the real particles are kept fixed while moving the trial ones). This energy minimization is performed using a bespoke algorithm, called $quickmaske$, which the $quickmin$ algorithm in LAMMPS except that each trial particles is independently relaxed. The fineness of the user-specified lattice determines the volume $\Delta V$ of the generic lattice cell. The energy minimum obtained for a trial particle is assumed to characterize all the particles that may nucleate anywhere else in $\Delta V$. In this way, the cumulative rate of all the $M$ sampled nucleation events is effectively a numerical integral in the sampled region, with resolution $\Delta V$. The most accurate result is obtained when $\Delta V\rightarrow 0$, but in practice a small but finite $\Delta V$ is already sufficient to estimate the cumulative nucleation rate accurately. How small $\Delta V$ should be, it depends on the employed interaction potentials and on the morphology of the simulated microstructure. An assumption underling this approach is that the energy minima alone fully control the precipitation rate cumulated over the whole $V_{box}$. In reality, thermal fluctuations would add contributions also from the locations surrounding the minima; future versions of MASKE may sample such fluctuations by relaxing the trial particles using molecular dynamics instead of energy minimization. F offers more discussion on interaction energy landscapes and their sampling for dissolution and precipitation events.

The flowchart in Fig. 3 depicts the algorithm for the kinetic simulations in MASKE. A brief description is given here, while more details can be found in E. The left side of the flowchart summarises the sampling of $FOT$ for all the discrete and continuous processes in a simulation. The shortest $FOT$ determines which type of process or event is to be executed next. If a continuous process is next, it is directly integrated numerically over its time period $\Delta t_{i}$ (user-specified in the input script). After that, the individual rates $r_{k,j}$ of the sampled discrete events are multiplied times $\Delta t_{i}$ and these products are stored in a per-particle array called $ratesT$, of size $N+M$. The $ratesT$ array is used in the next kinetic step for computing $FOT_{k}$. If the next event to execute is instead a discrete one, the $ratesT$ values are immediately updated, adding to each entry the product $r_{k,j}(FOT_{k}-t^{*})$. The $ratesT$ values are then used to select which discrete event to carry out specifically, out of all the possible ones. After executing a discrete event, the concentrations of the solvated molecules are updated accordingly, as discussed in the next section, and the $ratesT$ values are reset to zero. Subsequently, irrespective of whether the executed process is continuous or discrete, MASKE restores the initial position of all trial particles and executes a third type of possible processes, called $Every$. Unlike continuous or discrete processes, whose occurrence is associated to a time line, processes of type $Every$ are simply executed once every $n$ kinetic steps, where a step is one full loop in Fig. 3. $n$ is user-specified as an input, and presently MASKE implements only one type of $Every$ process, called $EmstoreLMP$, which executes a set of LAMMPS command that the user can specify and store in the input script. A typical use of such an $Every$ process is to execute a full energy minimization of the system with $n=1$, i.e. after each execution of a continuous or discrete process.

The algorithm in Fig. 3 is well suited to parallelisation, since each discrete event (i.e. each of the $N$ particle deletions and $M$ particle nucleations) can be sampled independently. Furthermore, spatial domain decomposition can be exploited during energy minimization. To exploit this, MASKE features a two-level parallelisation. On the first level, the user can specify multiple groups of processors, called $subcommunicators$, and attribute the sampling of one or more processes to each subcommunicator (e.g. one subcommunicator samples particle deletion, another one samples nucleation). On the second level, each subcommunicator uses all its processors to run any invoked LAMMPS command, e.g. energy calculations and minimisation using domain decomposition. This article does not report on the scalability and efficiency of the parallelisation, which is amenable to optimisation in the current version of MASKE. However, previously published works have already exploited both levels of parallelisation to simulate hundreds of thousands particles (Coopamootoo and Masoero, 2020) and multiple types of solid particles, all sampled for deletion and nucleation (Alex et al., 2023).

MASKE computes the rate of any discrete event, such as particle deletion or nucleation, based on two user-defined entities: a chemical reaction and a reaction mechanism, both specified in a chemical database text file, called $ChemDB$, whose syntax is presented in D. The chemical reaction is assumed to be one-step and is defined primarily through its formula and stoichiometric coefficients, e.g. :

for the dissolution of a calcium hydroxide molecule. The stoichiometric coefficients are specified in the $ChemDB$ file. The molecules appearing in the reaction formulas are also user-defined in the $ChemDB$ file, along with a set of useful parameters. For the $i^{th}$ solvated molecule, these parameters are its apparent volume $v_{app,i}$, hydrated radius $a^{o}_{i}$, ionic charge $z_{i}$, and the $b_{i}$ parameter for Debye-Hückel calculations of activity coefficients, discussed later. For solid molecular species, only linear sizes are required, which are used to compute their molecular volumes. $ChemDB$ gathers also the thermodynamic and kinetic parameters associated to each chemical reaction: standard free energy barriers $\Delta G^{o{\ddagger}}$, standard concentrations of activated complexes $c^{o{\ddagger}}$ (in the sense of Transition State Theory, which will be used below to express rate equations), equilibrium constants $K_{eq}$, and interfacial energies $\gamma$. See D for more details on this quantities.

The $ChemDB$ file also contains the definition of reaction mechanisms, which are the sets of assumptions that MASKE makes when assembling single chemical reactions to delete or nucleate a full particle. Indeed, a particle may be a coarse-grained entity whose deletion or creation may entail more than just a single reaction, e.g. $n_{r,j,diss}$ or $n_{r,j,nuc}$ reactions respectively to dissolve or nucleate a particle. Classical Nucleation Theory is an example of multiple individual reactions being assembled to express a reaction rate at the particle scale. In the current version of MASKE, the only implemented mechanism is called $allpar$ and it assumes that, irrespective of how big a particle is, its full deletion or nucleation only takes the time for one single-step chemical reaction to occur. This is exact if the volume $V_{p,j}$ of the generic particle $j$ equals the sum of the volumes of solid molecules dissolving or precipitating in one reaction. For bigger particles instead, the $allpar$ mechanism effectively assumes that the same chemical reaction takes places everywhere in $V_{p,j}$ at the same time. This is clearly unrealistic for particles beyond the nanometre or so, therefore in this article all the examples will feature particles whose $V_{p,j}$ entail $n_{r,j,diss}=n_{r,j,nuc}=1$. Hereafter this resolution is called ’unimolecular’, although in principle a single chemical reaction might also involve multiple solid molecules at once.

Whenever a particle is deleted or nucleated, a chemical reaction is carried out as part of an $allpar$ mechanisms. MASKE ensures mass conservation by consistently managing molecule exchanges between solid phases and solution, following the user-specified reaction formulas in $ChemDB$. Chemical reactions may entail differences in volume between reactants and products; these are accommodated in the $voidV$ volume introduced in the previous section, which may thus be positive or negative. The current version of MASKE does not simulate chemical reactions taking place in the implicit solution, hence it does not minimize the free energy of the solution itself. An interface between MASKE and the thermodynamic simulator PHREEQC (Parkhurst, 1995) is currently being developed, to address this limitation.

The rate equations for the $allpar$ mechanis in MASKE are rooted in Transition State Theory (TST), but with an extension to account for excess enthalpy coming from the mechanical interactions between particles. This excess enthalpy term, originally proposed in Shvab et al. (2017), is essential for capturing the effect of local morphology and stress state on reaction rates. The rate equations for one-step dissolution and precipitation reactions are:

Eqs. 3 and 4 will be hereafter called ’straight’ rate equations, as opposed to ’net’ rates defined later in this section. $\kappa$ is the probability of an activated complex evolving into reaction products; MASKE makes the common assumption that $\kappa=1$ (Lasaga, 2014). $k_{B}$ and $h$ are the Boltzmann and Planck constants, and $T$ is the temperature in degrees Kelvin. $c^{o{\ddagger}}_{d}$ and $c^{o{\ddagger}}_{pr}$ are the standard concentrations of activated complexes for the dissolution and precipitation reactions, $\gamma^{{\ddagger}}_{d}$ and $\gamma^{o{\ddagger}}_{pr}$ are their activity coefficients, and $\Delta G^{o{\ddagger}}_{d}$ and $\Delta G^{o{\ddagger}}_{pr}$ are their standard free energies of activation. $c^{o{\ddagger}}$ and $\Delta G^{o{\ddagger}}$ can be estimated experimentally, noting that the whole prefactors in Eqs. 3 and 4 defines the rate constant $k$:

Eq. 5 implies that the arbitrary choice of $c^{o{\ddagger}}$ must be compensated by the value of $\Delta G^{o{\ddagger}}$, so that a physical, univocally defined rate constant is preserved. The details of this compensation are discussed in F.

The $V_{t}$ term in Eqs. 3 and 4 is the tributary volume of the particle being sampled for deletion or nucleation, viz. the volume of the interaction energy basin, whose local energy minimum is where the sampled particle sits. The user of MASKE must specify $V_{t}$ in the $ChemDB$ file; for solid phases governed by short-range interactions, $V_{t}$ is typically in the range of one particle volume (see F for more discussion and examples of $V_{t}$). $\alpha$ is the spatial dimensionality of $c^{o{\ddagger}}_{d}$ or $c^{o{\ddagger}}_{pr}$; namely, $\alpha=3$ if the reactions are per unit volume or, as typical for dissolution/precipitation, $\alpha=2$ indicating reactions per unit surface. $\Delta V$ is the cell volume of the user-defined lattice of trial particles sampling candidate nucleation sites. When the KMC algorithm sums all the nucleartion rates for trial particles pertaining to the same interaction energy basin, the sum of all $\Delta V$’s equals $V_{t}$ (within the limits posed by the finite, and not infinitesimal, value of $\Delta V$). As a result, $V_{t}^{\alpha/3}$ in Eqs. 3 equals the sum of all $\Delta V\cdot V_{t}^{\alpha/3-1}$ pertaining to the same interaction energy basin. This aspect is discussed in more detail in G.

The square brackets in Eqs. 3 and 4 quantify the excess enthalpy of a particle, resulting from two contributions: $\Delta U_{d}$, which is the change in total interaction energy in the system caused by the particle deletion ($\Delta U_{pr}$ is the analogous term for a nucleation even), and$\gamma$, which is the solid-solution interfacial energy attributed to the sampled particle, and user-defined in the $ChemDB$ file. The $\Omega$ term, which multiplies $\gamma$ in the excess enthalpy, is the surface area of the sampled particle. $n_{r,j,diss}$ and $n_{r,j,nuc}$ appear in the excess enthalpy because the $allpar$ mechanism may be used for particles whose deletion or nucleation entails multiple chemical reactions; for the unimolecular resolution in this article, $n_{r,j,diss}=n_{r,j,nuc}=1$. The user-defined parameter $\chi$, with value between 0 and 0.5, specifies which fraction of excess enthalpy $U+\gamma\Omega$ is still present in the system when the reaction is in its activated complex state. $\chi=0$ is the usual assumption in TST, but results in this article will show that $\chi>0$ may sometimes improve the efficiency of a simulation.

The last terms to be defined in the rate equations are $Q_{r,d}$ and $Q_{r,pr}$, which are the activity products of the reactants in the dissolution and precipitation reactions. If $i$ is the generic solvated reactant:

where $c_{i}$ is the concentration in solution of the molecular species $i$, $\nu_{i}$ is its stoichiometric coefficient in the reaction, and $\gamma_{i}$ is its activity coefficient. MASKE computes $\gamma_{i}$ using the WATEQ Debye-Hückel formula (Truesdell and Jones, 1974):

$z_{i}$ is the charge of the $i^{th}$ species. $A$ and $B$ are solvent-specific constants (0.51 and 3.29 nm^-1 for water). $a_{i}^{o}$ is the hydrated radius of the molecule in solution, and $b_{i}$ is a molecule-specific dimensionless constant; both these parameters are provided by the user in the $ChemDB$ file. $I$ is the ionic strength of the fluid, featuring $n_{s}$ solvated species:

When $a_{i}^{o}$ or $b_{i}$ are unknown, MASKE computes the activity coefficients using established simplified methods (Langmuir, 1997), detailed in F.

The straight rate expressions in Eqs. 3 and 4 sample all the fluctuations between particle deletion and nucleation that, on average, guide the overall morphology evolution of the system. However, many such fluctuations may cancel out, especially near equilibrium, leading to inefficient simulations where the overall morphology makes very little progress on average. In such scenarios it may be convenient to remove the fluctuations by considering net rates instead, defined as the difference between straight rates:

These net rate equations are also implemented in MASKE, as a possible option for the $allpar$ mechanism. Their derivation is detailed in F. $Q_{p,d}$ and $Q_{p,pr}$ are the activity products of the products of the dissolution and precipitation reactions, with $K_{eq,d}$ and $K_{eq,pr}$ their equilibrium constants. The reactants of a dissolution reaction and the products of a precipitation reaction are solid phases, conventionally in standard state when stress-free, thus $Q_{r,d}=Q_{p,pr}=1$. G shows how the net rates in Eqs. 9 and 10 effectively scale linearly with the saturation index $\beta$ of the solution for the associated chemical reaction; namely, $r_{1,d}^{net}\sim(1-\beta)$ and $r_{1,pr}^{net}\sim(\beta-1)$, as usual in Transition State Theory. The definition of $\beta$ is:

noting that both the reactants of precipitation and the products of dissolution are simply the solvated species participating in the reaction. In a scenario where there is no excess enthalpy from changes in mechanical stress and interfacial energies, $\beta=1$ would imply equilibrium, i.e. no dissolution nor precipitation taking place, whereas $\beta>1$ would favour precipitation and $\beta<1$ dissolution.

The energy $U$ from mechanical interactions plays an important role in the reaction rates in MASKE, via the excess enthalpy terms $\Delta U_{d}-\gamma\Omega$ and $\Delta U_{pr}+\gamma\Omega$ in Eqs. 3 and 4, as implemented in the $allpar$ mechanism. All the examples in this article will employ a truncated harmonic potential of interaction between pairs of identical, spherical particles, $i$ and $j$:

plotted in Fig. 4. $r_{ij}$ is the interparticle distance. When the particles are at their equilibrium distance $r_{0}$, the interaction energy $U_{ij}$ is minimum and equal to the bond energy $-\varepsilon_{0}$. The stiffness of the interaction is dictated by the $k_{0}$ parameter. The cutoff distance for the interaction, $r_{c}$, is set here as the distance at which $U_{ij}$ returns to zero when $r_{ij}>r_{0}$; because of this definition, $r_{c}$ is not an independent parameter, but it is fixed once $r_{0}$, $k_{0}$, and $\varepsilon_{0}$ are decided:

$k_{0}$ is obtained from the Young modulus $E$ of the material of the particles:

The particle diameter $D_{0}$ is close but not equal to $r_{0}$, and their relationship depends on the modelling choice of the ‘limit crystal lattice’ of the solid phase. Depending on the type and range of the employed interaction potential, there are only several ordered arrangements, viz. crystal lattices, that are mechanically stable and where all the interacting particles sit at equilibrium distances. For instance, limit crystal lattices for the short-ranged, first-neighbour-only, radial interactions in Fig. 4, would be face centred cubic (FCC) or hexagonal close packing (HCP), both with $n_{b}=12$ first neighbours in the bulk. The mismatch between $r_{0}$ and $D_{0}$ emerges because any limit crystal lattice of spherical particles features a porosity and fills the space only up to a volume fraction $\eta<1$. For example, both FCC and HPC lattices have packing fraction $\eta=0.74$, hence a 26% porosity. However, if these limit lattices are used to discretise a non-porous continuum in a simulation, the maximum packing $\eta$ should represent a solid with same density as a non-porous phase, viz. with the same number of molecules, or mass, per unit volume. To exactly compensate for the artificial porosity $1-\eta$, Coopamootoo and Masoero (2020) have proposed to use an $r_{0}$ smaller than $D_{0}$ by:

The expressions of excess enthalpy in Eqs. 3 and 4 ($\Delta U\pm\gamma\Omega$) impose a dependence between interaction potential $\varepsilon_{0}$ and solid-solution interfacial energy $\gamma$, which must be respected in order to produce physically meaningful simulations. Consider a particle in kink position in the crystal lattice, which by definition has $n_{b}/2$ first neighbours. Precipitating or dissolving a stress-free particle in a kink position should not change the excess enthalpy of the system, because it creates the same amount of solid-solution interfacial area as it removes. Such a dissolution (and similarly, in reverse, for precipitation) can be modelled, first, as a detachment of the kink particle from the solid, which causes a loss of $n_{b}/2$ inter-particle bonds each with energy $\varepsilon_{0}$ (accounted for in the $\Delta U$ term), and second, as a dissolution of the detached particle which removes $\gamma\Omega$ interfacial energy. This leads to:

One example in this article will feature mechanical interactions between particles of two different types, with different sizes and representing different materials. For such scenarios, Alex et al. (2023) have proposed an extension of the truncated harmonic potential in Fig. 4:

The equilibrium distance $r_{12}=(r_{0,1}+r_{0,2})/2$ is the average of the equilibrium distances between identical particles for the two different types, 1 and 2. The stiffness term is expressed as $k_{12}=2\frac{k_{1}k_{2}}{k_{1}+k_{2}}$, i.e. idealizing the interacting particles as two springs in series, with individual stiffness as per Eq. 14. The bond energy is $\varepsilon_{12}=\frac{\gamma_{12}(\Omega_{1}+\Omega_{2})}{n_{b}}$, with $\gamma_{12}$ the solid-solid interfacial energy between the particles and $\Omega_{1}$ and $\Omega_{2}$ their individual surface areas. Eq. 13 is still valid for the cutoff distance, using $r_{12}$, $\varepsilon_{12}$, and $k_{12}$ in it.

A nanocrystal of calcium hydroxide, Ca(OH)₂, is considered in this example. The crystal is immersed in an aqueous solution of Ca²⁺ and OH^- ions at room temperature $T=298$ K. The crystal is made of 4,055 spherical particles, each with the volume of one Ca(OH)₂ molecule, $V_{p,CH}=\frac{V_{M,CH}}{N_{a}}=0.0557$ nm³, where $N_{a}$ is the Avogadro number and $V_{M,CH}=33.53$ cm³/mol is the molar volume of calcium hydroxide (indicated as $CH$ in the subscript). This brings to a particle diameter $D_{CH}=0.4737$ nm. The particles are initially arranged on an FCC lattice with cell size $a=\sqrt{2}\;r_{0,CH}$, which implies that $r_{0,CH}$ is the distance between first-neighbour particles in the lattice. Given the packing density $\eta=0.74$ of the FCC lattice, Eq. 15 provides $r_{0,CH}=\sqrt[3]{\eta}\;D_{CH}=0.4285$ nm.

The interaction potential between particles is the truncated harmonic one from Eq. 12. From the Young modulus $E_{CH}=48$ GPa (Wittmann, 1986) of calcium hydroxide, Eq. 14 provides $k_{0,CH}=19,745$ MPa$\cdot$nm. The bond energy $\varepsilon_{0,CH}$ comes from Eq. 16, having considered an interfacial energy between calcium hydroxide and surrounding solution $\gamma=68.4$ mJ/m² (Estrela et al., 2005) and $n_{b}=12$ first neighbours for a bulk particle in the FCC lattice. The interaction cutoff distance thus becomes $r_{c,CH}=0.457$ nm, from Eq. 13; this entails a failure strain $\frac{r_{c,CH}-r_{0,CH}}{r_{0,CH}}\approx 6.5\%$, which is reasonable at the molecular scale. With these interactions, an initially cubic grain of FCC-arranged calcium hydroxide particles is pre-dissolved as in Masoero (2023b), obtaining the rounded grain with diameter of ca. 7.3 nm in Fig. 5, which is the configuration for this example.

A very large volume of solution is placed in contact with the solid grain, by setting $\Delta V_{box}=10^{20}$ nm³. In this way, particle dissolution and precipitation negligibly change the concentration of solvated ions, which retain their initially assigned values ¹¹1The apparent volumes of solvated Ca²⁺ and OH^- ions, set to -0.025129603 nm³ and 0.001030404 nm³ respectively (Lothenbach et al., 2019) in the $ChemDB$ database, have a negligible effect in this example due to the large $\Delta V$ attached to the simulation box. The chemical reactions underlying dissolution and precipitation here are:

The value of $K_{eq,d}$ is taken from Lothenbach et al. (2019). The solvated ions are assigned the following parameters to calculate their activities as well as the ionic strength of the solution: $z_{Ca}=+2$, $a_{Ca}^{o}=0.486$ nm, $b_{Ca}=0.15$, and $z_{OH}=-1$, $a_{OH}^{o}=1.065$ nm, $b_{OH}=0.064$ (Lothenbach et al., 2019). In this way, the activity product of the solvated ions for Eq. 18 can be computed for any initially imposed set of concentrations. The explored concentrations range from 0 for both Ca²⁺ and $OH^{-}$ (saturation index $\beta=0$) to $c_{Ca}=39.1$ mmol/L, $c_{OH}=78.2$ mmol/L (i.e. pH $=12.89$ and $\beta=10$). Chemical equilibrium in stress-free conditions is expected at $\beta=1$, which implies $c_{Ca}=16.43$ mmol/L and $c_{OH}=32.86$ mmol/L (pH $=12.51$) for the fixed $c_{OH}/c_{Ca}=2$ ratio assumed here to ensure charge neutrality of the solution.

For the dissolution reaction in Eq. 18, the activation energy $\Delta G^{o{\ddagger}}_{d}$ and standard state concentration of the activated complex $c^{{\ddagger}}_{d}$ are both computed from the rate constant of Ca(OH)₂ dissolution in Bullard (2007): $k=0.66\;\mu$mol m^-2 s^-1. The conversion from $1\;\mu$mol m^-2 to number of reactions per $nm^{2}$, which are the simulation units, yields $c^{{\ddagger}}_{d}=10^{-6}N_{a}/10^{18}=0.6022$ nm^-2. Eq. 5 then provides $\Delta G^{o{\ddagger}}_{d}=74$ kJ/mol $=122.9$ ag nm² ns^-2 in simulation units. For the precipitation reaction in Eq. 18, $c^{{\ddagger}}_{pr}=c^{{\ddagger}}_{d}$ is assumed, thus $\Delta G^{o{\ddagger}}_{pr}=\Delta G^{o{\ddagger}}_{d}+k_{B}T\ln(K_{eq,d})=45$ kJ/mol $=74$ ag nm² ns^-2. The activity coefficients of the activated complex are taken as $\gamma^{\ddagger}_{d}=\gamma^{\ddagger}_{pr}=1$, as usual for reactions involving solids (Domínguez et al., 1998).

Deletion of existing particles and nucleation of trial particles are sampled using two different subcommunicators: $subA$ with 1 processor, and $subB$ with 2 processors. For particle nucleation, a cubic lattice of 50,656 trial particles is employed, with lattice spacing $r_{0,CH}/2$ in all three directions. The trial lattice is only created in a spherical shell region with inner and outer radii of 1.8 and 5 nm respectively, as shown in Fig. 5; this avoids the computational cost of sampling unlikely precipitation events far from the surface of the initial crystal, while leaving room for enough growth or dissolution to reach a steady state regime with constant average reaction rate. For the trial particles to move into a local minimum of interaction energy, 600 steps of energy minimisation are performed using the aforementioned $quickmaske$ algorithm, with time step of 0.45 ps and a maximum displacement per step of $r_{0,CH}/200$. Such a short minimization is appropriate because, for the interactions used here, a particle only needs to move by $r_{0,CH}$ or less, to find its closest local energy minimum. The ions added/removed to/from the solution after particle dissolution/precipitation are distributed between the volume of solution $V_{sol}$ in the simulation box and the additional $\Delta V_{box}$ in proportion to their volumes. However, since $\Delta V_{box}\gg V_{sol}$, changes in ion concentrations and saturation index saturation $\beta$ within the simulation box are negligible here. After each accepted dissolution or precipitation event, a process of type $Every$ is invoked, which minimizes the total interaction energy by displacing all the real particles in the system; this is done via 10,000 steps of the LAMMPS $quickmin$ algorithm, with timestep of 1 ps and a maximum displacement per step of 0.2 nm.

The nanocystal of Ca(OH)₂ from the previous example is now brought into initial contact with two platens, via a procedure detailed in Masoero (2023b): see Fig. 6.a. The platens are arbitrarily discretized with particles of diameter $D_{S}=0.6$ in an FCC lattice, hence with equilibrium distance $r_{0,S}=\sqrt[3]{\eta}D_{S}=0.5427$. The mechanical interactions between particles in the platens are of the form in Eq. 12, with elastic modulus typical of steel, $E_{S}=200$ GPa, which yields a high interaction stiffness, $k_{0,S}=208,397$ MPa$\cdot$nm A high interfacial energy is then assumed between platen particles and the solution, $\gamma_{S}=200$ mJ m${}^{-}2$, which leads to $\varepsilon_{0,S}=37.70$ MPa$\cdot$nm³. The resulting cutoff for the platen particles is $r_{c,S}=0.5617$ nm, viz. a strain at failure of ca. 3.5%. The mechanical interactions between platen and Ca(OH)₂ particles are of the form in Eq. 17, with equilibrium distance $r_{0,CH-S}=0.4856$ nm and interaction stiffness $k_{CH-S}=33,200$ MPa$\cdot$nm, obtained from the already computed, individual $r_{0}$ and stiffness of the two particle types. The platen-Ca(OH)₂ bond energy is set to $\varepsilon_{CH-S}=0$, making the interaction purely repulsive, i.e. $r_{c,CH-S}=r_{0,CH-S}$ and no bond under tension.

The volume of the simulation box is set to $V_{box}=8,000$ nm³, with an additional volume in contact with it $\Delta V_{box}=10,000$ nm³. The initial volume of the crystal is given by the initial number of Ca(OH)₂ particles times their individual volume, $V_{0,CH}=N_{0,CH}\cdot V_{p,CH}=4055\cdot 0.0557=226$ nm³, hence the initial volume of solution (in $V_{box}$ plus in $\Delta V_{box}$) is $V_{sol,tot}\approx 18,000$ nm³. The change in concentration of Ca²⁺ and OH^- ions in solution caused by the dissolution of one Ca(OH)₂ can be approximately estimated by neglecting the changes in solution volume caused by the reaction (MASKE accounts for them, through the volume of the Ca(OH)₂ particle and the apparent volumes of solvated Ca²⁺ and OH^- ions given in Section 2.5). For Ca²⁺, a dissolution event releases one molecule, hence changing $c_{Ca}$ approximately by $1/V_{sol,tot}\approx 92$ $\mu$mol/L; for OH^- ions the change is $2/V_{sol,tot}\approx 184$ $\mu$mol/L. Therefore, in this example the dissolution of Ca(OH)₂ particles during a simulation may significantly change the concentrations of ions in solution, and thus the saturation index $\beta$.

Only the Ca(OH)₂ particles are allowed to dissolve and to nucleate, via the same chemical reactions, thermodynamic and kinetic parameters, and trial precipitation lattice as in Example 1. The initial ion concentrations are set to $c_{Ca}=0.01643$ mol/L and $c_{OH}=0.03286$ mol/L, implying $\beta=1$. This means that no sustained dissolution nor precipitation is expected in stress-free conditions. All the particles in both platens are constrained not to move under the action of forces; this is done using the $setforce$ fix in LAMMPS. However, the particles in the top platen are rigidly displaced towards the bottom platen by 0.01 nm (via the $displace\_atoms$ command in LAMMPS) every $dt$ nanoseconds of simulated time in MASKE. This is done by using a continuous process of $mstoreLMP$ type in the MASKE input script. The chosen value of $dt$ controls the downward velocity of the top platen: $v_{z}=0.01/dt$ nm/ns. The initial distance between the centres of mass of the two platens is 9.12 nm, hence the imposed strain rate is estimated in $\dot{\epsilon_{z}}\approx 0.01/dt/9.12=1.1\cdot 10^{-3}/dt$ ns^-1.

The platens put the Ca(OH)₂ crystal under compression: see Fig. 6.b. The average compressive stress in the crystal, $\sigma_{z}$, is estimated as the per-particle virial stress in $z$ direction for each Ca(OH)₂ particles, $\sigma_{z,i}$ (computed in LAMMPS via a $compute$ of type $stress/atom$), averaged over all the $N_{CH}$ particles of Ca(OH)₂ existing at any time in the simulation:

$V_{p,i}$ is the volume of the i^th particle, j is the index of a neighbouring particle (of any type), $N_{n,i}$ is the number of neighbours mechanically interacting with particle $i$, $z$ is the Z coordinate of the two interacting particles, and $F_{z}$ is the Z component of the force on each particle stemming from the pairwise interaction. The presence of per-particle stresses entails a nonzero excess enthalpy in the rate equations for dissolution and precipitation. As a result, the compression is expected to trigger dissolution of Ca(OH)₂ even if the solution is initially at $\beta=1$. As the compressive strain increases, dissolution will help relaxing the increasing $\sigma_{z}$, but the concurring increase in ion concentrations in solution, and thus of $\beta$, will try to slow down dissolution and even cause precipitation of new Ca(OH)₂ particles in stress-free regions of the crystal surface. All this will result in a complex kinetic balance between solution chemistry and mechanical stress, which will depend on the imposed strain rate and generate a non-trivial evolution of the crystal morphology through dissolution/precipitation events and plastic deformations.

This section considers the Ca(OH)₂ nanocrystal described in Section 2.5, with the aim of assessing: (1) whether the simulations correctly predict the transition from crystal dissolution to growth that, for stress-free crystals, is expected as the solution goes from undersaturated ($\beta<1$) to supersaturated ($\beta>1$); (2) the variability in the results stemming from the random numbers involved in the Kinetic Monte Carlo algorithm (see B to appreciate the role of random numbers); (3) how the fineness of the lattice of trial particles may impact the results; (4) the pros and cons of using straight or net rates (viz. Eqs. 3 and 4 rather than Eqs. 9 and 10).

Fig. 7 shows precipitation and dissolution rates obtained from simulations using net rates. The rates are expressed per unit area of crystal surface $S$ and per unit time. To estimate $S$, knowing the volume $V_{p}$ of a Ca(OH)₂ particle and the number $N_{CH}$ of particles at any give time in a simulation, one can compute the total volume of the crystal $V_{p}N_{CH}$, and set $S=\pi(6V_{p}N_{CH}/\pi)^{2/3}$, having approximated the whole crystal as a sphere. Fig. 7.a shows a typical curve of precipitation rate (averaged over the 100 KMC events before the current one) vs. simulated time. The time scale extends to various seconds, despite the nanometre scale of the system: this ability to sample long time scales irrespective of the length scale being considered, is a strength of the KMC approach. Fig. 7.a also shows a steady state of constant rate in the initial part of the simulation; this happens for all $\beta\leq 1$ and $\beta\geq 8$ values, but not when $1<\beta<8$. The steady state values are plotted against $\beta$ in Fig. 7.b; for the $1<\beta<8$ range, the plotted rate is the one averaged over the first 100 KMC steps of a simulation. Fig. 7.b shows that: (i) dissolution (viz. a negative precipitation rate) is correctly predicted when $\beta<1$, as well as precipitation when $\beta<1$ and and chemical equilibrium (viz. zero rate) when $\beta=1$. The figure also shows that the sequence of pseudo-random numbers used in the KMC algorithm has a negligible effect on the average rates (‘seed 1’ and ‘seed 2’ indicate two difference sequences).

Transition State Theory (TST) predicts a linear relationship between rates per unit surface and $\beta$; by contrast the simulations, despite using rates that are based on TST, produce a nonlinearity in the $1<\beta<8$ range: see Fig. 7.b. This effect is due to the crystal morphology. An overall linear relationship, as per TST, requires that kink sites are not exhausted during precipitation; this is possible if there is a very large number of such sites in the initial crystal, which is not the case for the nanocrystal simulated here. Alternatively, new kink sites must be created rapidly, which is possible if $\beta$ is large compared to the energy barriers of precipitating particles with fewer neighbours $n$ than a kink site, $\Delta G_{n}=\varepsilon_{0}(n_{kink}-n)$, or if rare, energetically unfavourable fluctuations favour precipitation in under-coordinated sites, such as adatoms, which provide starting points for the growth of new layers. The latter scenario is impossible when net rates are used, as they average out all the energetically unfavourable fluctuations. As a result in Fig. 7.b, when $1<\beta<8$, precipitation is short-lived as it only consumes the initial kink sites. The continuous exhaustion of kink sites causes a continuous decrease in rate, which explains why a steady state is not established in said range. When $\beta\geq 8$, instead, new kink sites can be formed on certain crystal facets and precipitation can proceed further and faster. However, when the layers growing from facets on which new kinks can form at $\beta=10$ are exhausted, precipitation stops and the crystal ends up with the clear-cut morphology in Fig. 7.c. Even higher $\beta$ values would enable formation of new kinks on other facets too, thus sustaining further precipitation; a systematic exploration of these $\beta$ thresholds is not of primary interest here, but an indication of this process at play will be given later in Example 2. A close inspection of Fig. 7.b reveals a nonlinear increase in dissolution rate also as $\beta$ approaches 0. This happens because a very low $\beta$ enables the dissolution of sites that are more coordinated than kinks, creating vacancies (or ‘pits’) on the crystal surfaces where new dissolution fronts can emanate, accelerating the overall process.

The effect of energetically unfavourable fluctuations is shown in Fig. 8.a,b, with simulation results from using straight rates. The $\chi$ parameter in the rate equations controls how prevalent the effects of fluctuations are. The explored values are $\chi=0$ (as usual in Transition State Theory), which attributes all the excess enthalpy from mechanical interactions to the dissolution reaction only and maximizes the fluctuations, to $\chi=0.5$, which equally distributes the excess enthalpy between dissolution and precipitation, minimizing fluctuations but not removing them (as opposed to net rates). Fig. 8.a shows that all the results from straight rates predict a linear relationship between overall rate and $\beta$, also in the $1\leq\beta\leq 8$ range where net rates were unable to nucleate new crystal layers. Straight rates instead always enable the formation and stabilization of adatoms (e.g. those circled in Fig. 8.b) which allow for indefinite growth. A close inspection of Fig. 8.a reveals a slight shift of equilibrium (i.e. rate = 0) towards $\beta$ values that are slightly greater than 1; this is the effect of relying on an estimated value of $V_{t}$ when explicitly sampling both dissolution and precipitation events (here $V_{t}$ was set to $V_{p,CH}\frac{r_{c,CH}}{r_{0,CH}}$).

All the simulations with straight rates, $\beta>1$, and $\chi=0.5$ produce the same morphology as in Fig. 8.b, which is reached when the entire nucleation lattice is covered by real particles; further growth is expected if bigger lattices were used. The same morphology is expected to emerge when $\chi=0$ and $\chi=0.05$, but in this case the simulations were stopped after a few thousands atoms precipitated, because of the high computational cost of sampling many unlikely fluctuations. Indeed, the precipitation of ca. 500 particles at $\beta=10$ required ca. 10,000 KMC steps in simulations with $\chi=0$ and 8,000 steps when $\chi=0.05$, as opposed to 800 steps with $\chi=0.5$, and 500 steps for net rates. On the other hand, Fig. 8.a shows that the precipitation rates compute with $\chi=0$ are significantly underestimated in all the other simulations with $\chi>0$, as well as when net rates are used. This confirms that rare fluctuations are important for quantitatively capturing the nucleation of new layers. Therefore, when using MASKE one should judge on a case-by-case basis whether to use straight or net rates and, in case, which value of $\chi$, in order to balance the requirements of running rapid calculations, predicting qualitatively correct mechanisms, and obtaining quantitatively correct rates.

A last parametric study here concerns the impact of the fineness of the nucleation lattice on the resulting average rates. For simulations using net rates, Fig. 8.c shows the temporal increase in $N_{CH}$ at $\beta=10$ when cubic nucleation lattices of different spacing $\zeta\cdot r_{0,CH}/2$ are employed. All simulations with lattice spacing less than $r_{0,CH}/2$, viz. with $\zeta\leq 1$, return very similar results ²²2The simulation with $\zeta=0.33$ was stopped earlier because it entailed the sampling of ca. 1.4 million trial particles, as opposed to ca. 50,000 when $\zeta=1$ and ca. 6,000 when $\zeta=2$. Coarser lattices however, with $\zeta>1$, underestimate precipitation because they miss energetically convenient minima in the interaction energy landscape. This indicates that a nucleation lattice spacing of $r_{0,CH}/2$ or less is appropriate for the simulations and the mechanical interaction potentials considered here.

In this example, the Ca(OH)₂ crystal from the previous section is surrounded by a relatively small volume of solution, which is initially at chemical equilibrium with the crystal itself, viz. $\beta=1$. The crystal is then subjected to a compressive strain with fixed rate, which should favour the dissolution of particles under high mechanical stress. The dissolving particles increase the concentration of Ca²⁺ and OH^- in solution, hence its saturation index $\beta$ too; this may cause precipitation of new Ca(OH)₂ particles in regions of the crystal surface that are under low mechanical stress. At high imposed strain rates, plastic deformations may occur too.

Fig. 9 shows the deformation and stress field of the Ca(OH)₂ crystal when compressed at different strain rates $\dot{\epsilon}_{z}$. The KMC algorithm in MASKE allows exploring a broad range of strain rates, here from $10^{8}$ s^-1 to $10^{-4}$ s^-1; this is valuable because the experiments typically impose rates in the order of the lowest values here, which are difficult to reach with other simulation techniques (e.g. molecular dynamics; Cao et al. (2019)). All the simulations in Fig. 9 employed net rates; some considerations on the effect of straight rates will be made later in this section.

At low strain rates, $\dot{\epsilon}_{z}\sim 10^{-4}$ s^-1 in Fig. 9, the dissolution triggered by the per-particle stress (via the excess enthalpy in the rate equations) is fast compared to the deformation process. As a result, the dissipation of local stress $\sigma_{z,i}$ through dissolution is efficient, therefore the local stress remains relatively low even at high $\epsilon_{z}$, as shown by the dark violet color in the $\epsilon_{z}=40$% snapshot in Fig. 9. The same snapshot shows significant recrystallization, driven by the increase of $\beta$ in solution eventually triggering precipitation of new Ca(OH)₂ in low-stress regions.

Higher strain rates induce higher per-particle stresses $\sigma_{z,i}$: see $\dot{\epsilon}_{z}=1.1$ and $\sim 10^{8}$ s^-1 in Fig. 9. For compatibility of displacements, the rate at which a layer of particles in contact with the platens dissolves must equal the velocity at which the top platen moves. A higher strain rate entails a higher platen velocity, which requires a higher dissolution rate, achieved by allowing a slightly higher strain before dissolution, which in turn increases the local stress and thus the excess enthalpy term in the rate. Fig. 9 also shows that little recrystallization takes place at high strain rates. This is because the precipitation rate scales linearly with $\beta$, and therefore with the number of dissolved particles (assuming that no reprecipitation takes places at all, which is a good approximation at high strain rates). At a given strain level $\epsilon_{z}$, the number of dissolved particles is basically the same irrespective of the imposed strain rate (for sufficiently high $\dot{\epsilon}_{z}$). By contrast, the strain rates have been increased by orders of magnitude here, as well as the dissolution rates required to accommodate the strain; in the resulting short timescales of the deformation/dissolution processes, precipitation events become rare. In the simulation with highest strain rate, $\sim 10^{8}$ s^-1, the animation of the full chemo-mechanical process in Fig. 9 shows that plastic shear slips along crystalline planes occur at high strain levels, $\epsilon_{z}\gtrsim 30\%$. Such slips are not recorded at lower strain rates, where the dissolution-induced relaxation of local stresses prevents the activation of plastic deformations.

Fig. 9 also shows results for a purely mechanical response, i.e. with both dissolution and precipitation turned off. The simulation was conducted at a low strain rate $\dot{\epsilon}_{z}\sim 10^{-4}$ s^-1, but the same result would be obtained at higher strain rates because the particles are brought to static equilibrium between two successive execution of the platen-moving continuous process in MASKE (static equilibrium is imposed via a MASKE process of type $Every$ performing energy minimization in MASKE). Therefore inertial effects do not contribute to rate-dependence here. The purely mechanical response in Fig. 9 is somewhat similar to the coupled chemo-mechanical one at high $\dot{\epsilon}_{z}\sim 10^{8}$ s^-1, but with some important differences too. At low strain levels (see $\epsilon_{z}=10\%$ in Fig. 9), even the limited dissolution taking place at $\dot{\epsilon}_{z}\sim 10^{8}$ s^-1 is sufficient to reduce the local stresses compared to the purely mechanical test. At $\epsilon_{z}\sim 20\%$ and $40\%$ the purely mechanical test displays a wider low-stress region near the surface of the crystal; the particles in this region do not contribute significantly to the overall mechanical response to the compression. Partial dissolution at $\dot{\epsilon}_{z}\sim 10^{8}$ instead moulds the crystal layers in contact with the platens, removing particles with high $\sigma_{z,i}$ that induce stress concentration, and thus enabling a more uniform mechanical response. Another distinctive aspect of the purely mechanical test is that extensive plastic deformations from shear slips already start at $\epsilon_{z}\approx 17\%$, and a fully developed central shear band controls the entire deformation response at $\epsilon_{z}\gtrsim 30\%$. In the dissolving crystal with $\dot{\epsilon}_{z}\sim 10^{8}$, instead, small and isolated, plastic shear slips only start appearing at $\epsilon_{z}\approx 30\%$.

Fig. 10.a shows the stress strain curves computed at different strain rates $\dot{\epsilon}_{z}$. A higher $\dot{\epsilon}_{z}$ produces a higher average compressive stress $\sigma_{z}$; this reflects the already-discussed need for faster dissolution (whose rate depends on local stresses) to ensure compatibility of displacements at the platen-crystal interface.

A counter-intuitive result in Fig. 10.a is that crystals that are allowed to dissolve, when compressed at high rates, are stronger than crystals that are not allowed dissolve. This is indicated in Fig. 10.a, by the higher peak $\sigma_{z}$ of the curve for $\dot{\epsilon}_{z}\approx 10^{8}$ s^-1, compared to the ‘mech. only’ curve. The ‘mech. only’ curve lacks the mechanism of particle dissolution, which helps relaxing the local stresses on particles while ensuring compatibility of displacement at the crystal-platen interface. Such compatibility, in the ‘mech. only’ simulation, relies entirely on mechanical deformations; these are initially elastic and, when local stresses become too high, plastic shear slips are triggered and relax the local stress. One may thus expect that the ‘mech. only’ simulation, being more constrained, would produce a higher $\sigma_{z}$ than simulations where dissolution is possible too. This is indeed the case at low strain levels: see $\epsilon_{z}<0.2$ in Fig. 10.a. However, the first extensive shear slips in the ‘mech. only’ simulation occur already at $\epsilon_{z}\approx 17\%$. After that point, subsequent shear slips effectively cap the peak stress $\sigma_{z}$ to ca. 2’000 MPa. By contrast, in the $\dot{\epsilon}_{z}\approx 10^{8}$ s^-1 simulation, dissolution of the most stressed particles limits the maximum local stress, which gets redistributed as a lower stress to the particles neighbouring the dissolved one. All this leads to a more uniform stress distribution in the crystal (see Fig. 9) while also delaying the triggering of plastic slips until $\epsilon_{z}\approx 30\%$. As a result, the peak stress $\sigma_{z}$ in the $\dot{\epsilon}_{z}\approx 10^{8}$ s^-1 is ca. 2’500 MPa, greater than in the ‘mech. only’ case.

Another finding from Fig. 10.a is the stepwise stress-strain curve at the lowest strain rate, $\dot{\epsilon}_{z}\approx 10^{-4}$ s^-1. Specifically, $\sigma_{z}$ increases with $\epsilon_{z}$ for a while, then it oscillates around a constant value while $\epsilon_{z}$ keeps increasing, until a new regime of increasing $\sigma_{z}$ is entered, and so forth. Three such steps can be counted in Fig. 10.a. This behaviour is the result of the kinetic chemo-mechanical balance established during the simulations. At low rates $\dot{\epsilon}_{z}$, no plastic deformations are triggered and particle dissolution is the only mechanisms to accommodate the externally imposed $\dot{\epsilon}_{z}$. When compression starts, the first dissolution events release ions in solution, raising the saturation index $\beta$ for Ca(OH)₂ precipitation: see the increasing $\beta$ at $\epsilon_{z}<0.15$ in the $\dot{\epsilon}_{z}\approx 10^{-4}$ s^-1 curve in Fig. 10.b. With the net rate equations used here, Example 1 in the previous section has shown that no significant precipitation is possible when $\beta\lesssim 8$, hence $\beta$ keeps increasing. The increasing $\beta$ slows down dissolution but, since the overall rate $\dot{\epsilon}_{z}$ is imposed and must be matched, a compensating increase in dissolution rate is made possible through an increase in mechanical stress (which increases the excess enthalpy term in the dissolution rate equation). This explains the regimes of increasing $\sigma_{z}$ in Fig. 10.a, which indeed correspond to regimes of increasing $\beta$ in Fig. 10.b. At $\beta\approx 8$, particle precipitation consumes any further ion that dissolution releases in solution, so that the average rate (precipitation minus dissolution over 100 KMC steps) oscillates around zero: see the curve for $\dot{\epsilon}_{z}\approx 10^{-4}$ s^-1 in Fig. 10.c when $0.15<\epsilon_{z}<0.3$. As a result, $\beta$ now remains constant while $\epsilon_{z}$ increases, and the constant overall rate $\dot{\epsilon}_{z}$ implies that the stress $\sigma_{z}$ must remain constant too. This explains the matching plateaus in Figs. 10.a and 10.b. The first plateau continues until exhausting all the surface sites where precipitation at $\beta\approx 8$ is possible: cf. the final morphology in Fig. 7.b. Then a new regime of increasing $\beta$ and $\sigma_{z}$ is entered, until a new plateau occurs at $\beta\approx 50$, when precipitation on other exposed crystal surfaces becomes possible too. The new plateau again corresponds to a net dissolution rate oscillating around zero, in Fig. 10.c. All this describes a complex, kinetic equilibrium between composition of the solution, stress state in the crystal, and microstructural evolution.

The stepwise shape of the stress-strain curve is lost at higher strain rates, as shown by the $\dot{\epsilon}_{z}\approx 1$ s^-1 curves in Fig. 10. Additional results omitted here show a progressive reduction and disappearance of the plateaus for intermediate rate values between $\dot{\epsilon}_{z}\approx 10^{-4}$ s^-1 and 1 s^-1. This means that, despite some particles do still precipitate when $\beta\gtrsim 8$, they cannot compete any more with the faster rate of the dissolving particles. Consistently, when $\dot{\epsilon}_{z}\approx 1$ s^-1 the saturation index $\beta$ increases monotonically in Fig. 10.b and the average dissolution rate never goes to zero in Fig. 10.c. The $\beta(\epsilon_{z})$ curves in Fig. 10.b for $\dot{\epsilon}_{z}\approx 1$ and $10^{+4}$ s^-1 are similar, meaning that dissolution progresses at a similar rate. This is expected because, in a scenario where precipitation is negligible, the total amount of dissolved ions depends only on the accommodated strain $\epsilon_{z}$. If $\beta$ is approximately identical for $\dot{\epsilon}_{z}\approx 1$ and $10^{+4}$ s^-1, then $\sigma_{z}$ is necessarily greater at the higher imposed rate, since only $\sigma_{z}$ now can speed up the dissolution rate to match $\dot{\epsilon}_{z}$. Less intuitively, the $\beta(\epsilon_{z})$ curve for $\dot{\epsilon}_{z}\approx 10^{8}$ s^-1 in Fig. 10.b is initially similar to those for $\dot{\epsilon}_{z}\approx 1$ and $10^{+4}$ s^-1, but it significantly diverges from them at $\epsilon_{z}\approx 0.2-0.3$, leading to less concentrated solutions at high strain levels $\epsilon_{z}$. The reason for this lies in the triggering of plastic deformations at high rate, $\dot{\epsilon}_{z}\approx 10^{8}$ s^-1, which cap the stress level (see $\sigma_{z}$ at $\epsilon_{z}\gtrsim 0.3$ in Fig. 10.a) and provide another mechanism to accommodate the imposed strain. With a capped $\sigma_{z}$, the increasing $\beta$ reduces the dissolution rate, which in turn reduces the slope of the $\beta(\epsilon_{z})$ curve in Fig. 10.b. The resulting interdependence between plastic deformations and composition of the solution is another non-trivial result of the chemo-mechanical coupling in MASKE.

All the results in this section are from simulations where the $allpar$ mechanism in MASKE employed net rates of dissolution and precipitation, as per Eqs. 9 and 10. Additional tests were conducted using straight rates too, as per Eqs. 3 and 4, with $\chi$ values between 0 and 0.5. At high strain strain rates, $\dot{\epsilon}_{z}\gtrsim 1$ s^-1, straight and net rates produce almost identical results; this is expected, since the main effect of using straight straight is to favour precipitation, which plays a negligible role at high $\dot{\epsilon}_{z}$. Instead, at low $\dot{\epsilon}_{z}\lesssim 1$, simulations using straight rates produced lower $\beta(\epsilon_{z})$ than those from net rates at same $\dot{\epsilon}_{z}$. This agrees with the findings from Example 1 in the previous section, where straight rates were shown to enable precipitation at lower $\beta$ values. Therefore, simulations with straight rates are indeed expected to produce plateaus in $\beta$ and $\sigma_{z}$ at lower $\epsilon_{z}$ than simulations with net rates. These aspect, however, should be studied further because the efficiency loss from using straight rates significantly limited the maximum $\epsilon_{z}$ that could be reached during the simulations. For the purpose of this paper, the results presented here for net rates have already demonstrated how MASKE can predict the emergence of a complex kinetic balance between chemical and mechanical processes over long time scales.