Numerical Phase-Field Model Validation for Dissolution of Minerals

: Modelling of a mineral dissolution front propagation is of interest in a wide range of scientiﬁc and engineering ﬁelds. The dissolution of minerals often involves complex physicochemical processes at the solid–liquid interface (at nano-scale), which at the micro-to-meso-scale can be simpliﬁed to the problem of continuously moving boundaries. In this work, we studied the diffusion-controlled congruent dissolution of minerals from a meso-scale phase transition perspective. The dynamic evolution of the solid–liquid interface, during the dissolution process, is numerically simulated by employing the Finite Element Method (FEM) and using the phase–ﬁeld (PF) approach, the latter implemented in the open-source Multiphysics Object Oriented Simulation Environment (MOOSE). The parameterization of the PF numerical approach is discussed in detail and validated against the experimental results for a congruent dissolution case of NaCl (taken from literature) as well as on analytical models for simple geometries. In addition, the effect of the shape of a dissolving mineral particle was analysed, thus demonstrating that the PF approach is suitable for simulating the mesoscopic morphological evolution of arbitrary geometries. Finally, the comparison of the PF method with experimental results demonstrated the importance of the dissolution rate mechanisms, which can be controlled by the interface reaction rate or by the diffusive transport mechanism. Calculations for this research were conducted on the Lichtenberg high performance computer of the TU Darmstadt.


Introduction
Mathematical modelling of the moving-boundary dissolution fronts of minerals is important in a wide range of engineering technologies. For example, it is of great importance in fields of geochemistry, materials science, hydrometallurgy, etc. Predictions of the moving boundary dissolution phenomena can support in the design of engineering processes where dissolution is desired: e.g., in extraction of elements or reactivity of cementitious minerals, but also when not desired, e.g., in durability (corrosion) issues of building materials (e.g., steel-reinforced concrete frames). In general, minerals dissolve when exposed to aggressive solution environments and form leached layers of varying density and strength [1]. This in turn affects the mechanical and transport properties of the microstructure which further may be relevant at higher scales, for example when the material (rock, concrete or mortars) has structural applications. In addition, the dissolution mechanisms of some special minerals can be of great industrial and environmental interest. For example, the dissolution of scorodite is considered a potentially good carrier for arsenic fixation [2]. Moreover, the application of innovative self-healing concrete in civil engineering has been extensively and intensively studied in recent years. The dissolution of Ca(OH) 2 from the concrete matrix is one of the key processes of the durability and self-healing mechanisms [3][4][5].
The dissolution of minerals often involves complex physico-chemical processes at the solid-liquid interface. However, this can be simplified at the mesoscale to the problem of a continuously moving boundaries. Traditional sharp interface models are thus required temperature gradient at the interface and proved that in dealing with the solidification problem of pure melting, the PF parameters can be accurately determined under the thin interface limit. Based on this model, Xu and Meakin [60] developed a phase-field approach for aqueous dissolution/precipitation reactions assuming first order reaction kinetics. The model was validated by a one-dimensional analytical solution of interface motion due to solute precipitation. Two additional terms were added to the diffusion equation, one corresponding to the discontinuity of the solute concentration gradient, at the interface, while the second one represents the net source (or sink) of the solute, coming from the discontinuity in the solute concentration across the interface. In most of the studies, the values of interface mobility are used as empirical or hypothetical ones [61][62][63][64][65][66][67]. Furthermore, few attempts have been made to explain in detail the calculation of the interface mobility and its relation to other physical parameters [51,[68][69][70][71]. Therefore, tackling of this difficulty will be one of the innovations of this paper.
In contrast to the reaction kinetics controlled case, here we further validate the PF approach on the experimental results (of NaCl dissolutions) and are focusing mainly on the diffusion limited mechanisms. First, the one-dimensional diffusion-controlled dissolution problem will be simulated using an analytical solution and the classical KKS model, separately. The results will then be compared to clarify the estimation and interaction of the interface mobility with other PF parameters. The effect of solid particle shape on the dissolution process is 2D analysed and validated on literature data for NaCl dissolutions. The PF results are then validated against the data obtained from analysis by the video-microscopy images and compared with the analytical model. Finally, a concluding discussion on the whole article is given.

Types of Dissolution
Chemical dissolution of minerals occurs as a congruent or an incongruent reaction, depending upon the type of a mineral [72]. Congruent dissolution of a solid mineral is a chemical reaction which completely dissolves the mineral and all products of this reaction are dissolved species. An obvious example would be calcite CaCO 3 and NaCl [73,74]: when the primary solid phase is altered and at the same time a secondary solid phase is formed, incongruent dissolution occurs, for example the alteration of albite to gibbsite NaAlSi 3 O 8 , or Kaolinite [75,76], which requires a more advanced thermodynamic modelling approaches to be integrated in the PF: Figure 1 shows the dissolved diffusion process of soluble minerals based on the diffusion interface. When dealing with the problem of moving boundaries, the conventional approach separates the different phases by a sharp interface. The interface movement is solved by a partial differential equation describing, for example, mass and thermal diffusion equations. These equations have to be combined with boundary conditions of varying values and positions. When some variables (heat flux or concentration) cross the sharp interface, jump discontinuities can occur, making the calculation very difficult. In the PF model, the interface is described as a diffuse interfacial layer with smooth transitions. Thus, the phase transformation is represented by a change in an order parameter (φ). As shown in Figure 1, the solid phase is represented by "1" while the liquid phase by "0", hence the order parameter varies continuously between 0 and 1 at the solid-liquid interface. Schematic representation of solute concentration of soluble minerals in situ (i.e., green dotted line) and equilibrium states (i.e., red solid line), and the phase transformation within a diffuse and sharp interface, respectively. For the sharp interface, the evolution of the solute concentration is discontinuous at the interface. However, for the diffuse interface, the solute concentration evolves continuously between their equilibrium values at the mineral (c Se ) and solution boundary (c Le ).

Diffusion-Controlled Dissolution Mechanisms
As the congruent dissolution process occurs, the solute is gradually transferred into the solution, the length of the solute base phase decreases, and the solid-liquid diffuse interface gradually moves toward the inside of the solute base phase. The solute concentration in the initial solution is c L0 , while the solute concentration at equilibrium is c Le . The solute concentration in the solid is kept constant when the diffusion phenomenon in the solid is not taken into account. The change of solute concentration in solution with time is related to its position. The solute concentration (c L ) increases with time away from the solute matrix phase; solute concentration (c IS ) decreases with time near the solute matrix phase. The solute concentration in between (c I ) is in a state of dynamic increase or decrease. x Se , x S0 , x Le and x L0 indicates the positions corresponding to the above solute concentrations.

Analytical Solutions
A following 1D planar analytical description of a diffusion-controlled dissolution for phase transformation is considered, where the solid is immersed in a (semi-)infinitive liquid solution. The diffusion equation is a parabolic partial differential equation, which is expressed as follows: where D L is the diffusion coefficient, c(x, t) is the concentration at location x and time t, subject to the conditions: where x = R is the position at the solid-liquid interface, c L0 and c Le represents the initial concentration and the equilibrium concentration of one component in the liquid phase. At the solid-liquid interface, the following independent flux balance condition must be fulfilled: where c S is the concentration in solid phase which is taken as a constant.
The exact analytical solution for the field is [6]: The interface position at current time can be expressed as where R 0 denotes the value of R at the time t = 0, while where λ is given by Different from the planar solid where exact solution is available, only approximately analytical solution model for the diffusion-controlled dissolution of the spherical solid has been found. The stationary-interface approximation is expressed as [6]: where the current interface position is R 2 = R 2 0 − βD L t; β as defined in Equation (10). The implicit expression for the particle radius ratio y (y = R/R 0 ) with respect to time is defined as follows: where,

The Phase-Field (PF) Method
The total free energy of the thermodynamic system drives changes in the microstructure of materials only when the total free energy changes from a high chemical potential (or a higher free energy) state to a low chemical potential (or a lower free energy) state to eventually attain an equilibrium. The total free energy F of the system is a function of the solute concentration c and the phase parameter φ and expressed as: where F loc is the local free energy of the system, F int is the interfacial energy, and k is the gradient energy coefficient. In order to simplify the numerical calculation, the molar concentration of the solute is normalized by the molar concentration of the solid c S , that is, c = c /c S . The molar concentration of the solid is defined as the density of the solid divided by its average molar mass [77]. Each point in the entire domain is a mixture of two phases with different chemical compositions. The Gibbs free energy expression, in the KKS model, has been widely employed in solidification mechanisms of binary alloys [78,79], in addition to the recent extension to the field of electrochemical corrosion [61,68]. The mechanism of the corrosion reaction is similar to that of the dissolution reaction, i.e., both are phase transformations triggered by the diffusion of ions. Hence, the double well potential (i.e., Gibbs free energy density) has two minima at φ = 0, φ = 1 and a maximum at φ = 0.5 (middle of the interface). Based on this theoretical basis, the present model identifies the local free energy f loc (c, φ) as a fractionally weighted average of the solid f S (c S ) and liquid free energies f L (c L ), and imposes a double-well potential ωg(φ) as follows: where the interpolation function h(φ) is built as h(φ) = −2φ 3 + 3φ 2 , and ω is the height of the double-well potential function given by The free energy density of the solid and liquid phase is approximated by a parabolic function with the same curvature A as follows: where c Se = c S /c S = 1 and c Le = c sat /c S are the solute concentrations at the normalized equilibrium of the solid and the liquid phase, respectively. Complementary condition indicates that the phase concentrations are constrained such that the chemical potentials of each phase are equal: The solute composition in the interface area is the fraction-weighted average of the liquid and solid composition, and the same formula is used for the diffusion coefficient, as shown below: The interfacial evolution is controlled through coupled conserved and non-conserved dynamics. Particularly, the Allen-Cahn equation is used to describe the temporal evolution of the non-conserved variable φ. However, the diffusion equation is applied for solving the evolution of the conserved parameter c: where L is the PF mobility: where D is the diffusion coefficient. The energy of the system is minimal when it reaches equilibrium. In order to find the PF profile φ 0 (x) and the composition c 0 (x) at equilibrium state, Kim, et al. [48] deduced a one-dimensional solidification problem with boundary conditions φ 0 | x→−∞ = 1 (solid) and φ 0 | x→+∞ = 0 (liquid). Since the equilibrium state means the vanishing of the driving force, the PF profile φ 0 (x) should satisfy the following Equation: Thus, by combining with the double well equation, the PF profile can be expressed as: Then, the composition is:

Problem Description and Model Tests
In this section a 1D congruent dissolution case study is presented selected as benchmark for verifying the soundness and capability of the proposed PF procedure.

Benchmark with Analytical Model for One-Dimensional (1D) Congruent Dissolution
The model consists of a one-dimensional domain with a size of 20 mm. The solid and liquid domains have a 3:17 ratio (see Figure 2). This is chosen to ensure that the length of the solution must be long enough for diffusion to take place in a system whose domain is considered to be semi-infinite.  Table 1.  Most of the available experimental data on solid dissolution in the literature focus on recording the evolution of solute concentration [11,14,80,81] or dissolved solid mass over time, and do not explicitly observe the movement of the solid-liquid boundary [80,81].
The study of Quilaqueo and Aguilera [95] is one of the few studies that provides detailed experimental dataset to be used for the experimental validation of the PF numerical models. They performed image analysis by coupling a digital camera to a stereo microscope to obtain microscopic images of the dissolution process. Recording started by placing a single NaCl particle in 500 µL of water without stirring at 20 • C. The time profile of dissolution was obtained by calculating the projected area of the single crystal as a function of dissolution time from the video microscope image.
Based on these experimental data, one-dimensional simulations of the dissolution process of a single salt particle performed, in a spherical coordinate system, and twodimensional simulations of the dissolution process of three different shapes of NaCl particles, namely round, ellipsoidal and irregular, have been performed by using the PF model. Neumann no-flux boundary conditions are applied for 1D simulation. Periodic boundary conditions are used for 2D simulations. The thin interface limit is supported by the chosen experimental case of the highly soluble mineral crystals, which can be considered as non-porous. This results in negligible solid-liquid thickness (at the mesoscopic scale) and it has been considered in this paper. The employed parameters are summarized in Table 2. For easier implementation of energy equations at microscale, all length dimensions were normalized by the solution radius and energy terms normalized by A (see Appendix B).

Summary of Modelling Assumptions
In the analytical solution and the PF model, following simplifying assumptions are made: 1.
The diffusion coefficients of aqueous species in solids and in solution are constants, respectively; 3.
The diffusion of all aqueous species is expressed in terms of a single ionic concentration;

4.
The solubility of NaCl in solution is independent of particle size.

Parameterization
The relationships between material properties (interface thickness l 0 and interface energy σ) and PF parameters (coefficient of PF gradient k and double-well potential ω) are discussed in Appendix A. The derivation of the interface mobility L, under the thininterface thickness condition, is shown in Appendix A. The curvature of the free energy density function A can be determined from the Gibbs free energy, i.e., ∆G (Appendix A).

Parameter Normalization
Normalization of the model parameters is one of the important steps of data preprocessing. A series of input values are normalized to the range [0, 1] according to Appendix B, in order to let models converge effectively.

Finite Element Implementation
Numerical implementation of the PF model is carried out by using finite element method in the framework of multiphysics object-oriented simulation (MOOSE) environment [28]. Transient solver with preconditioned Newton's method was used. In this case, the full and accurate Jacobian was calculated. The backward Euler algorithm was employed. Adaptive time stepping was used to improve computational efficiency. The time step would grow or shrink according to the number of iterations taken and needed to obtain a converged solution in the last converged step. The maximum number of nonlinear iterations per time step was also set to provide optimal solution efficiency. For 2D simulation, the triangular element type was chosen to mesh the geometry.
In addition, an adaptive mesh refinement (AMR) was used [97]. Based on the error estimated from the FEM results, the global and local mesh errors were calculated, and then the mesh size was automatically adjusted to the changing morphology of the grain boundaries at each time step. This is very effective for the numerical solution of partial differential equations in regions of arbitrary shape. In order to ensure numerical stability and simulation accuracy, at least 5 nodes on the diffusion interface were used to describe the boundary morphology, while coarser grids were used for solutions and solids that were far from the boundary (Figure 3). This does not only provide an accurate representation of the boundary evolution, but also improved the computational efficiency. The relative and absolute error tolerance was set to 10 −8 .

Central Processing Unit (CPU) Computation
The analytical experiments were performed using an Intel Core i7-6500U central processing unit (CPU) 2.5 GHz with 8 GB of RAM and MATLAB R2017a (64-bit). For the PF computation, a parallel computing was achieved by using Open Multi-Processing (OpenMP 42.0.51) on a High Performance Computer.

PF Validation against Analytical Solution for a Dissolution of Planar Mineral
Under diffusion-controlled dissolution conditions, if the initial thickness of the diffuse interface (l 0 ) is known, the interface mobility L can be determined using Equation (A10). Due to the lack of relevant experimental data upon the values of l 0 , a parametric study of L was carried out. Under the assumption of a thin interface condition [48], the value of l 0 should be taken much smaller than the minimum radial dimensions of the initial solid phase; however, from a computational point of view, it is expected that the thickness of the interface has to be as large as possible in order to keep the interface from being overly densely meshed, which increases the computational effort. Therefore, three cases of initial interface width, i.e., 1 × 10 −8 (PFM1); 1 × 10 −5 (PFM2) and 1 × 10 −4 (PFM3), were tested, corresponding to 0.0003%, 0.33% and 3.33% of the initial length of the solid phase. In addition, there must be at least 5 to 10 grid points in the interface area to ensure the stability of the numerical calculation and the reliability of the results [28]. Three cases of V 0 (1 × 10 −6 (PFM4); 1 × 10 −8 (PFM5) and 1 × 10 −10 (PFM6)) were tested. Figure 4 shows a comparison between the analytical (diffusion limited) model and six cases of PF models, where the reaction rates are slower than in case of diffusion control. As l 0 decreases, L increases (as they are inversely related by Equation (A10)), causing a faster dissolution reaction, till reaching a limit defined by a diffusion control. At 8 h, PFM1 dissolves at a thickness 1.15 times greater than that of PFM3. The result of PFM2 is in good agreement with that of the analytical model. It can also be seen that the slope of the dissolution curve becomes progressively smaller with dissolution time due to the diffusion-controlled dissolution, i.e., the overall rate of dissolution slows down as it is governed by the diffusion flux and thus dependent on the concentration gradient that reduces with saturation of the solution. However, the (slow) reaction-controlled mechanism (see curves PFM4, PFM5 and PFM6) approaches linearity as the V 0 decreases, and the lower the V 0 , the slower the dissolution speed. The above results show that the PF model developed is capable of describing both the diffusion-controlled and the reaction-controlled dissolution. Under the thin interface limit condition, the agreement between the analytical (diffusion-controlled) results and the PF model for the diffusion-controlled dissolution is in satisfactory agreement with the converged solution (PFM2). A slight overestimation by the PFM1 model may be attributed to the used approximation (Equation (A10)) to numerically approach the diffusion-limited case. In this approximation, V 0 is approximated as D L /l 0 (where D L is the diffusion coefficient of the solute in the solution). The slight difference could also be the result from small incompatibility issues between the employed thermodynamic parameters, namely the used NaCl interfacial energy in PF model, and the NaCl solubility constant (p). Overall, we argue that the obtained agreement is overwhelming considering that no fitting calibration of the parameters has been performed. Figure 5 shows the spatial distribution of solute normalised concentration over time. As solid dissolution starts, the concentration of solution c L at the diffusive solid-liquid interface rapidly reaches saturation concentration (i.e., c Le equilibrium state), while the solute normalised concentration in the solid phase keeps constant at 1 (i.e., at initial concentration). The solutes form a diffusion layer at the thin solid-liquid interface and continue to enter (diffuse) into the bulk solution. This results in a gradual decrease in the width of the solid phase. The concentration of solutes in the solution is gradually increasing, and the increase of the concentration near the solid-liquid interface is particularly significant due to diffusion limited transport through the solution. In turn the concentration of solute smoothly decreases to zero from the interface zone to the right end of the solution. This confirms also that the dissolution is carried out in a system that is regarded as semi-infinite space, which is only controlled by a diffusion mechanism without any significant effects of the imposed boundary conditions (which deviate from the idealized semi-infinitive case). Figure 6 shows the movement of the interface over time. As the dissolution proceeds, the interface gradually moves toward the solid phase. It can be seen from the reduced width of the solid phase that the speed of dissolution starts faster and then slows down, again due to the limited diffusion process of the solute through the exposure solution.

The Effect of Mineral Shape: Dissolution Simulation by Two-Dimensional (2D) PF Model
Irregular particles can be simplified by spherical shapes provided certain conditions are met. Numerical simulations are then performed using the spherical symmetry and coordinate system, reducing the model to only one space dimension (1D model). However, some studies have shown, both theoretically and experimentally, that the particle shape and surface roughness may affected the dissolution rates [98][99][100][101][102]. Therefore, before adopting the spherical simplified 1D model for NaCl crystals in this study, the effect of particle shape on the dissolution rate was analysed. In this way one can determine whether the simplification of the spherical shape for NaCl particle having some circularity factors (0.71 ± 0.06) as in used literature data [95] is reasonable. Therefore, Figure 7 shows the 2D dissolution simulation results over time for spherical, elliptical and irregular shapes but with the same area. Qualitatively, the sharp edges of irregular shapes gradually disappear during the initial stages of dissolution, and the curvature decreases until they are completely rounded. The ratio of the major axis to the minor axis of the ellipse gradually decreases and develops towards the circular direction. The radius of the circle decreases gradually and the curvature remains constant. It can be seen that the dissolution under all shapes follows the process of spheroidization. The edges of the particle become smoother during dissolution. As the dissolution reaction proceeds, the morphological differences between the particles with different shapes become smaller. The reason for this trend is that among closed geometries of equal area, circles have the smallest circumference. This means that the total interfacial free energy of the solid-liquid is minimal. Irregularly shaped solids have a high solid-liquid total interface free energy due to their uneven boundaries. A high interfacial energy means a high total free energy of the system. The system always tends to decrease the total free energy, this being an important factor in determining the mineral shape during dissolution. The area of the solid-liquid interface tends to decrease, which causes the flange at the solid-liquid interface to disappear and eventually to become round.
Since the solid-liquid interface is described in the PF model as a diffuse interface, with a certain width, it is difficult to describe the interface length in terms of the absolute perimeter of the solid phase. However, in order to characterize the change in the surface morphologies of crystals for different shapes with time, the contour length of φ = 0.5 is taken to approximate the solid-phase perimeter. Figure 8 shows that at the initial condition, the irregular mineral grain has the maximum interfacial length. As the dissolution reaction proceeds, the interfacial length gradually decreases and the curves of the ellipse and circle almost coincide. The normalized phase ratio φ * is defined as φ * = (φ t − φ min )/(φ max − φ min ), being φ t the integration of the phase at the time t, over the domain, while φ min and φ max are representing the min and max integration of the phase, respectively (for c as well). The profile of the normalized phase ratio physically represents the projected area of particles. It can be seen from Figure 9 that the projected area of the three shapes of particles decreases with time, while the normalized concentration ratio (c * ) keeps constant, which proves the conservation of mass for solute transport. The slope of φ * becomes progressively smaller. This is because as the solid phase dissolves, the interfacial area decreases. The contact area between the solute source and the diffusion-solution zone is getting smaller. This results in a decreasing solute flux to the solid surface, which leads to a progressively slower dissolution rate. It can also be seen from this figure that the φ * profile of the circles and ellipses basically overlap. The dissolution rate of irregular shapes before 200 s is slightly faster than that of circles and ellipses, while after 200 s, the three curves overlap and dissolve completely at the same time. This is because the perimeter of the irregular shape is much larger than that of the circles and ellipses, which exhibits a faster rate at the beginning of dissolution. As the dissolved shape tends to be round with the lowest interfacial energy, the circumference between the three shapes becomes similar ( Figure 8). Therefore, in the later stages of dissolution, the dissolution curves of the three shapes are coincident.    Figure 11. Phase profiles of single NaCl spherical particle along radial direction with time. Figure 12 shows a comparison of the analytical solution, the PF model results and experimental results regarding the dissolution rate of individual NaCl crystals. The undissolved area is calculated from the residual solid phase length in the 1D simulation. The analytical solution is slightly below the lower boundary of the experimental values. In order to ensure that the dissolution reaction rate is completely controlled by the diffusion, the length (volume) of the exposure solution must be large enough so that semi-infinitive conditions are met, corresponding to the analytical solution. In that case no increase in concentration of solute should occur at the system (solution) boundary point. The length of solution 0.3 times (PFM_D1), 0.4 times (PFM_D2) and the original length (PFM_D3, 4.92 × 10 −3 m) were tested using the PF method. From the results it can be seen that the dissolution rate slows down and the solute concentration increases significantly at the solution boundary point as the solution phase length becomes shorter ( Figure 13). The change in solute at the boundary point, calculated using the solution lengths in the experiment (PFM_D3), is almost zero. It was thus verified that the diffusive dissolution of individual NaCl crystals can be simulated well using our implementation of the PF method.  The experimental result in Figure 12 shows lower initial rates, due to the reaction controlled mechanism. Initially, the diffusion flux is very high, as the solute concentration of the initial exposure solution is zero (initial condition). Therefore the overall reaction rate is limited by the reaction rate which is lower than the initial high diffusion flux. Such reaction rates in the PF model can be considered in the interface mobility (L) that should vary with time and is fundamentally a function of the Gibbs free energy of the chemical reaction or the solute concentration. In future, it should be attempted to physically represent the interface mobility kinetics (L) as a function of the solute-under saturation, as commonly used in reaction rates expressions. Such a more fundamental approach is still missing in PF, due to the complexity in its mathematical derivation. Here, L is adjusted in a simplified way as a smaller value (4.11 × 10 −6 ) within 100 s and a larger value (4.11 × 10 −3 ) after 100 s whose result is represented by PFM_R1. As can be seen from the comparison, the dissolution rate is relatively flat at the early stage of dissolution when process controlled by the reaction. The curve of PFM_R1 is higher than that of PFM_D1, D2 and D3. However, after 100 s, PFM_R1 almost overlap with the other three due to the conversion of the dissolution rate control into the diffusion control mechanism.

Conclusions
Based on the results of this study, the following conclusions can be summarized: • by comparing with the results of the analytical method, it is verified that the PF model can accurately handle the dynamic evolution of the general diffusion-controlled phase transformation process; • using NaCl as an example, the PF model can successfully simulate the mesoscopic evolution of inorganic non-metallic materials caused by diffusion-controlled dissolution. Using the derived interfacial mobility, the PF numerical simulation results show accurate and consistent agreement with the analytical method results, as well as with the experimental ones derived with video-microscopy images analyses. It is worth mentioning that all the input parameters of the PF model have real physical meaning and are based on the experiments data; • an observed discrepancy was related to the dissolution mechanism, which was found to be initially limited by the reaction rate, being slower than the diffusion flux due to the rapid change of solute concentration. This change in dissolution mechanism was successfully captured by adjusting the PF interface mobility (L). • the dissolution characteristics of NaCl particles with different circularity factors were analysed by the 2D PF model. The simplification of spherical shape for NaCl particles was verified to hold.
In future studies, the reaction control and diffusion control mechanisms will be combined with the second law of thermodynamics and non-equilibrium thermodynamics with respect to L, so that L can be represented as the function of solute concentration or the Gibbs free energy of the reaction. In addition, numerical simulations need to be implemented at higher space dimensions to allow the introduction of complex microstructures in mineral particles, such as pores, grain structure and surface roughness, so that their impact on dissolution kinetics can be assessed. The above results confirm that the dissolution kinetics of mineral particles can be successfully simulated using the employed PF model and the irregular morphological evolution can be effectively simulated in 2D. It is worth noting that the surface morphology of irregular particles has a strong influence on the dissolution kinetics. The complex evolution of particle morphology in physicochemical processes can be accurately evaluated only in a full 3D system. Furthermore, a dynamic (apparent) diffusion coefficient should be explicitly taken into account, e.g., as a function of concentration, which is of critical importance in analyzing the diffusion-controlled dissolution through porous materials involving additional chemical interactions. However, highly soluble (NaCl) crystals, as investigated in this work, can be considered as non-porous which leads to assume that their diffusion coefficient can be kept constant. It is worth mentioning that the proposed PF approach can be also extended for simulating the opposite processes of those presented in this paper, namely the precipitation of NaCl. This is part of forthcoming research and will be undertaken by enriching the current free energy density function of the liquid phase through adding an extra precipitation term, i.e., ∆ r f . This should be considered for further development to develop a comprehensive mineral dissolution and precipitation modelling tool.

Appendix A.2. The Interface Mobility L
The interface velocity V is typically expressed as product of a factor involving the thermodynamic driving force for dissolution and a kinetic factor involving the interface mobility [103], i.e., where the kinetic factor V 0 corresponds to the limiting velocity under infinite driving force (forward reaction rate), ∆G is the Gibbs free energy difference, between the free energy of solid and liquid which is responsible for the interface displacement. The driving force of dissolution ∆F is given by ∆F = f L c e L − f S c e S − c e L − c e S f S c S (c S ), (−∆G = ∆F). Thus, Equation (A5) is obtained by Taylor series expansion and approximated as: Using the derivation in the KKS model, Equation (23) can be expressed in terms of ∆F as follows: Under 1D instantaneous steady state, φ is derived for the position as follows: Under the thin interface limit condition, Equation (A8) combined with the equilibrium phase expression (dx/dφ 0 = − √ k/ √ 2ω[1/φ 0 (1 − φ 0 )]) and Equation (A1) modifies into: It should be noted that this equation only holds if the diffusion potential of solute in the interface region is constant. In addition, if the dissolution process is controlled by diffusion, V 0 is usually approximated as D L /l 0 ; where D L is the diffusion coefficient of the solute in solution [31]. In general, the determination of V 0 is very difficult [104]. In summary, L can be derived from Equation (A9): The value of A can be derived from ∆G. Following [61], A makes the free energy of the solid phase is equal to ∆G when c is 0.4 (see Figure A1).

Appendix B. PF Parameters Normalization
In order to ensure the convergence of the model on the mesoscale and improve the calculation efficiency, all variables are normalized. All molar concentrations are normalized by the solid molar concentration c S , and all length scales are normalized by the length of the domain (l). According to the total energy function in Equation (9), three equations are proposed as follows: The normalized variables are presented as follows: where, p is the solubility of sodium chloride at 20 • C and 1 atm, which needs to be converted to units of mol/m 3 in the calculation.