Deciphering pore-level precipitation mechanisms

Mineral precipitation and dissolution in aqueous solutions has a significant effect on solute transport and structural properties of porous media. The understanding of the involved physical mechanisms, which cover a large range of spatial and temporal scales, plays a key role in several geochemical and industrial processes. Here, by coupling pore scale reactive transport simulations with classical nucleation theory, we demonstrate how the interplay between homogeneous and heterogeneous precipitation kinetics along with the non-linear dependence on solute concentration affects the evolution of the system. Such phenomena are usually neglected in pure macroscopic modelling. Comprehensive parametric analysis and comparison with laboratory experiments confirm that incorporation of detailed microscale physical processes in the models is compulsory. This sheds light on the inherent coupling mechanisms and bridges the gap between atomistic processes and macroscopic observations.

Interaction between fluids and solids, including the precipitation and the dissolution of minerals from and into aqueous solutions in porous media, has a complex feedback to the solute transport in the aqueous phase itself and vice versa 1,2 . There is clear evidence that the macroscopic properties of the medium depend on the physical and chemical transformations that occur at the micropore scale in a strongly non-linear way. For example the formation of nanometer to micrometer-sized precipitate layers around preexisting mineral grains, dramatically reduces the permeability of the medium, thus affecting the macroscopic flow through it, although the total porosity does not change significantly. Flow alteration has in turn an effect on mineral precipitation-dissolution processes. The result is a fully coupled hydro-geochemical reactive transport process. The involved physical mechanisms play a key role in geochemical processes and industrial applications that span from geothermal energy and oil industry to pharmaceutical products, catalysts and long-term nuclear waste containment [3][4][5][6] . The prediction of the evolution of such systems is very challenging and can be assisted by conducting dedicated laboratory experiments along with numerical reactive transport modelling [7][8][9][10] .
Macroscopic modelling and simulations use a simplistic description of the dynamic processes that actually take place at the microscale. The true evolution of the reactive surface area as well as the alteration of pore size and topology is usually unknown. The change in permeability and diffusivity of a porous medium, due to reactions, is correlated to a global change of the bulk porosity, using empirical laws with adjustable parameters 11 . When a relevant experiment is modelled using a pure macroscopic reactive transport code, the major macroscopic experimental observations can only be reproduced after fitting the adjustable parameters using the same experimental results that are intended to be modelled 9 .
In the case of macroscopic solvers the porosity changes are simply modelled by applying mass balance equations between the inlet and outlet flow rates of the domain of interest. This averaged porosity is then used in the parametric Kozeny-Carman equation to derive an overall permeability and in the parametric Archie's law to obtain the effective diffusivity. It is therefore not surprising, that smooth simplified relations frequently fail to predict accurately the temporal evolution of the system. It is the microscopic pore level physics per se that ultimately controls the macroscopic processes. Macroscopic models cannot circumvent the lack of microscopic physical description without empirical tuning to the experimental data. Conversely, power laws and correlations, of the porosity variation effect on the permeability and diffusivity, can be derived from microscopic simulations 12,13 and further up-scaled for use in macroscopic algorithms.

Reactive transport laboratory experiment, macroscopic modelling and characterization
In ref. 9 we have carried out a reference experiment to elucidate competitive dissolution-precipitation reactions in a granular porous medium. This kind of system is relevant to baryte scale formation in industrial applications at large scale geochemical systems 14 (hydrocarbon exploration and extraction, geothermal heat extraction, uranium mining, and technologically enhanced naturally occurring radioactive materials-TENORM waste). Further, an exceptionally complete set of thermodynamic and kinetic data is available. The experimental setup, consisted of a strip of reactive celestine (SrSO 4 ) with bimodal grain size distribution (<63 µm and 125-400 µm) embedded in a comparatively inert matrix (quartz sand). A barium chloride solution (0.3 M BaCl 2 ) was injected into the flow cell, leading to fast dissolution and partial replacement of Celestine (zero ionic strength solubility product constant = − . ) 15 thus modifying the hydraulic and structure properties of the medium (porosity, permeability, connectivity). A subsequent post mortem analysis of the reactor using optical and electron microscopy as well as synchrotron-based micro-X Ray Fluorescence (µ-XRF) and micro-X Ray Diffraction (µ-XRD) measurements identified two distinct baryte precipitates: i) nano-crystalline baryte filling the pore space and ii) thin baryte overgrowths coating large and medium-sized celestine crystals 16 . The distinction between the two phases was based on the analysis of Debye rings from the μ-XRD patterns after integrating over microscopic sample volumes. The formation of complete Debye rings from tiny volumes of a few μm 3 indicates that these regions contain a multitude of randomly oriented nanometer-size baryte particles. It should be noted that the absence of Ba-Sr solid solutions is confirmed by the same XRD analysis. This finding is in agreement with thermodynamic calculations and experimental observations 16 predicting stability of nearly pure baryte and celestine phases 17 at very low Sr/Ba ratios as used in this experiment. The microscopic characterization shows that two different precipitation mechanisms were operating during the experiments: homogeneous nucleation (HON -nanoparticle aggregate) and heterogeneous nucleation (HETepitaxially grown rim) both highly dependent on solution metastability (threshold supersaturation). In Fig. 1(a), a Back-Scattered Electron (BSE) microscopy image with the major characteristics of the reacted zone is shown. Large celestine crystals are covered by epitaxially grown baryte layers (HET). Baryte nano-crystalline phase (HON) is also distinguishable and black color corresponds to the remaining open pore-space after passivation of the zone.
The modelling of the experiment using a macroscopic approach could reproduce the evolution of the system only after fitting the kinetic and transport parameters in order to match the experimental results 9 . This level of modelling is not able to explain the presence and formation of two distinct baryte phases nor can provide detailed information of the system evolution, and therefore has limiting predictive capability. This analysis demonstrates the importance of pore-scale understanding of new phase precipitation in porous media. It is clear that in order to improve the predictive capability of the geochemical computational algorithms, such details and effects that span across length and temporal scales have to be incorporated in the modelling.

Cross scale modelling concept
The basis of our modelling is the simultaneous description of nucleation kinetics, a process that occurs at the atomic scale and of mass transport through realistic porous media. Pore geometry and topology dynamically change due to dissolution and precipitation requiring robust pore-level flow models and solvers. The nucleation mechanisms which lead to the precipitation of ionic solids and particularly of baryte are discussed in detail in ref. 18 and can be explained based on the Classical Nucleation Theory (CNT). We apply here this concept, in spite of current discussion on the general applicability of CNT [19][20][21] . As pointed out by ref. 22 , CNT cannot, for example, be applied to carbonate mineral precipitation due to the formation of amorphous precursor phases in this system. CNT is however still a valid concept to explain precipitation of sulfate minerals and particularly baryte, which unlike carbonate minerals does not form allotropic precursors.
In Fig. 1, our cross-scale modelling concept is illustrated. It links the nucleation kinetics and subsequent growth of nuclei that occurs at the atomic level, Fig. 1(b), with the solute concentration subjected to advection and diffusion through the porous structure, Fig. 1(c). In Fig. 1(b) a schematic plot of the total Gibbs free energy (ΔG) for HON and HET is shown as a function of the radius of spherical nuclei. As soon as a critical nucleus is formed and the energetic barrier is exceeded (at the maximum of the ΔG curve) further growth of the nucleus is more favorable than its disaggregation and precipitation starts. For HON nucleation to occur a higher activation barrier compared to HET has to be exceeded thus explaining why homogeneous nucleation requires higher levels of supersaturation to occur. When critical size nuclei are formed precipitation proceeds as growing spheres whose surface area and growth rate is traced throughout the simulation as depicted in the inset of Fig. 1(b). This process is controlled by the larger scale flow and solute features and vice versa. In Fig. 1(c) a generic pore-level flow simulation is depicted where colors represent the velocity magnitude (top) and signify for example the major preferential flow paths which change as the system evolves (the white streamlines in the right-bottom Fig. 1(c)). A model prediction for high saturation index is also visualized in Fig. 1(d) where celestine crystals are depicted in blue and total baryte precipitates are in yellow.
Mineral precipitation via HON from supersaturated solutions is inhibited in porous media compared to precipitation in free solutions 22,23 . This can be explained by a confinement effect of the solution, which makes, according to CNT, formation of critical size BaSO 4 nuclei more unlikely, the smaller the pore volume 22 . In larger pores the probability of forming supercritical nuclei is higher, and therefore induction times are much shorter. Because of variable pore sizes and non-uniform flow rate distributions, the local concentrations vary in the micrometer range. A Supersaturation-Nucleation-Time (SNT) diagram for homogeneous (HON) and heterogeneous (HET) nucleation optimized for the modelled reactive experiment is plotted in Fig. 2(a) to support the interpretation of the results of Fig. 2(b). In our modelling we dynamically assign, at every point of the domain, the value of its corresponding pore-size. The dependence of induction time on saturation index (SI), for pore sizes of 1 µm and 10 µm is also shown. SI is defined as the decimal logarithm of the ratio of ionic activity product in solution, IAP  For the description of simultaneous transport, dissolution and precipitation processes along with the resulting changes of the microscopic pore geometry we applied the lattice Boltzmann method (LB) 12,[24][25][26][27] . The LB method is a computational technique based on a special discretization of the kinetic Boltzmann equation and provides the necessary fluid dynamic framework (See Methods section). The combination of pore-level transport modelling with CNT can explain the evolution of the considered system and allows predictions under different conditions. For the present problem the discretization grid (lattice) of the porous domain is of the order of 1-10 μm. Sub-micrometer modelling as shown schematically in Fig. 1(b,c), allows to relate the solute transport and concentration with the parameters governing the HON and HET precipitation kinetics, such as the size of critical nuclei, the reactive surface area, the induction time, the local pore-volume and the growth rate. These quantities are locally defined and set the pace for baryte precipitation thus determining the evolution of the system.

Results and Discussion
The computational domain is setup using the BSE-SEM image of Fig. 1 as a template for the large celestine grains and filling the rest of the space stochastically with small celestine grains, until the initial porosity of the experiment was reached (ε = 0.33) as depicted in the inset of Fig. 2 The second case mimics the actual conditions of the experiment where 0.3 M BaCl 2 solution was constantly supplied via advection 9 . Simulations results are depicted on Fig. 2(b). For SI = 4.08 and SI = 3.96 baryte is precipitating mostly via the HON mechanism and the majority of the dissolved sulfate ions precipitate during the production of the nanocrystalline phase. All small celestine crystals are consumed. Further reduction of the SI increases the induction time of HON and thus favors the epitaxial growth of baryte rims (HET). The system then evolves at a slower pace since there is much less surface available compared to the HON case. In the plot of SI = 3.8 HON precipitation is more pronounced at the regions where larger pore-spaces exist (>50 μm). For SI = 3.6 nucleation proceeds strictly via the HET mechanism for the given reaction time.
The evolution of the system after 140 hours of reaction depends non-linearly on the SI and is depicted in Fig. 3(a). The black curves represent the four aforementioned cases. The plotted curves represent the percentage of conversion of small celestine crystals to baryte as a function of reaction time. The complexity of interconnected physical and chemical processes that take place at the pore scale as well as the effect of the degree of suspersaturation on the precipitation pathway is analyzed in Fig. 3(b), in terms of the relative amounts of HON and HET precipitation. We note that porosity changes are proportional to baryte precipitation. Macroscopic models that implement simple analytical relations would fail to predict the effect of these underlying processes.
Precipitation rates measured in the lab in batch experiments often greatly deviate from observations in the field. Reactive surface area, eventual pore topology and flow properties can only partially explain this deviation 10 , thus effects related to the size of the pores must be invoked. The importance and sensitivity of pore size dependent kinetics is illustrated with the dashed colored lines in Fig. 3(a). For the conditions of the experiment where SI = 3.96 we conduct an analysis to highlight the sensitivity of nucleation kinetics with relation to the pore sizes. We perform two additional calculations by keeping all parameters fixed and identical, except for the probability of forming critical nuclei with respect to the pore-size. For that, we increase the probability of forming the critical nuclei in all pores in such a way, that it corresponds to the probability associated to ten times larger (SI = 3.96(10x)-red curve) or ten times smaller pores (SI = 3.96(/10)-green curve). For example, in a 10 μm pore the induction time for the onset of precipitation is tuned as if the respective pore were 100 μm in the first case (10x), and 1μm in the second case (/10). This alteration of the induction times results in a completely different global reaction rate and distinct evolution paths. It demonstrates that kinetic parameters derived from batch experiments cannot always be transferred directly to reactive transport models of porous systems, especially when the average pore size is at the micrometer level. Thermodynamic data acquisition from laboratory experiments must be carefully correlated to the respective pore size distribution.
Our modelling allows monitoring also the evolution of the baryte rim thickness and its growth rate which is also strongly dependent on the SI as depicted in Fig. 4. At high SI, the HON mechanism is producing nanocrystalline precipitates at a much faster rate compared to the time needed for sulfate ions to diffuse and reach the growing rims, resulting in thinner baryte layers. Colored lines represent the state of the evolved system after the respective reaction time in hours. In order to compare our predictions with the actual experimental results we compare the epitaxially grown baryte rim thickness after the same reaction time (140 hours). The sample of the BSE-SEM image of Fig. 1(a) was analyzed digitally and the average rim thickness was calculated. The average rim thickness for SI = 3.96 is predicted from our model to be L sim = 4.5 μm and is in very good agreement with the microscopic post-mortem analysis of the experiment which measured L exp = 4.2 µm.

Summary and Outlook
Dissolution and precipitation processes and mass transport phenomena in porous media are coupled at pore level in a complex non-linear way. Macroscopic modelling approaches do not capture microscopic details of this process and therefore face a great difficulty in predicting evolution of such systems 9 . We introduced a cross-scale modelling approach that couples nucleation kinetics based on classical nucleation theory and pore level mass transport on the lattice Boltzmann framework. This coupled description of transport, dissolution and precipitation allows the simulation of complex systems, with a large amount of degrees of freedom at an unprecedented physical detail, provided the necessary parameters are known as in the case of baryte. Quantities such as the size of critical nuclei, the reactive surface area, the induction time, the local pore-volume and the local growth rate have been considered and a good agreement with experimental observations demonstrate the potential of this approach. Moreover, this improves the predictions of the evolution of systems where experiments are too difficult or even impossible, due to limitations in space and time. At the same time it can support and enhance the interpretation of experimental results.
The proposed method has a broad area of industrial and natural sciences applications that range from radioactive scale formation during exploitation of underground energy resources to pollutant dispersion and purification of water 28 . A feature challenge would be to apply this methodology in more complex geochemical systems including solid solutions. Further development would be needed to describe systems where precipitation proceeds along non-classical nucleation pathways.

Methods
where k is the Boltzmann constant, T is the absolute temperature and Γ a pre-exponential factor. The nucleation energy barrier is highly dependent upon the interfacial tension σ and saturation index SI. This barrier is higher for HON compared to HET as shown schematically in the lower part of Fig. 1 of the main text. This explains why HON appears only at relative high values of SI. The exact relation reads: where β is a geometry factor, set to 16.8 for a sphere, σ is set to 0.134 J m −2 for HON according to 29 , ν is the volume of a single spherical BaSO 4 monomer, and is set to 8.60 × 10 −29 m 3 22 . For HET, the interfacial tension σ HET , was calculated as a function of the contact angle θ between the substrate and the growing mineral, according to 18 : HET H ON 2 1 3 The pre-exponential factor in Eq. 1 represents the rate of attachment of monomers controlled by diffusion and is given as Γ = 2πZDN 1 N 0dc , where D is the diffusion coefficient of monomers set to 9.3 × 10 −9 m 2 s −1 and d c = 4σv/ kTln(SI). Parameters N 0 and N 1 are concentrations that represent the numbers of monomers per unit volume of fluid and the number of nucleation sites respectively. After 22 these parameters were set to N 0(HON) = 3.33 · 10 28 m −3 and N 0(HET) = 2.50 · 10 13 m −3 . Parameter N 1 depends on the supersaturation which has been evaluated using the Gibbs Energy Minimization Software -GEMS 30 . Z is the Zeldovich factor given as where V p is the volume of the local pore under consideration. Note that the implemented algorithm measures and updates the size of the respective pores at every time step, which changes due to dissolution (size increase) or precipitation (size decrease). This is necessary for the correct calculation of the necessary the local induction time needed to pass before the onset of precipitation. Baryte precipitation rates are calculated according to the empirical kinetic law proposed in ref. 32 . This law was derived from systematic precipitation experiments both in the presence and absence of precipitation inhibitors: = . × − − − , and A [m 2 m −3 ] is the relative surface area per unit volume. For HON, a sublattice model with respect to the transport solver is implemented as discussed in the main text. The sum of molecules that form the first nucleation clusters are converted to equivalent spheres. The resulting reactive surface area is calculated based on the evolving area of these virtual spherical objects. More than one critical nuclei are allowed to be simultaneously formed, and their growth is traced throughout the simulation (See Fig. 1(b) of main text).

Lattice Boltzmann modelling details (LB).
Pore-level processes can be simulated with a variety of numerical methods [33][34][35] . For the modelling of the advection-diffusion and precipitation processes a multi-component LB model is used. The model is composed of a basis fluid medium that recovers the Navier-Stokes equations at the macroscopic limit, plus several passive scalar coupled population sets that simulate the diffusion of ions. The isothermal guided equilibrium nine-velocity model (D2Q9 lattice) of ref. 36  The following population-moments correspond to the density of the solution ρ and the momentum j a in the direction a = x, y: The guided equilibrium populations f i eq are given in a closed form, where T 0 = 1/3: i eq a x y ia c ia ia a a , The Boltzmann BGK equation is solved: , where τ is the relaxation parameter and T 0 µ τρ = is the resulting macroscopic dynamic viscosity. BGK stands for the Bhatnagar-Gross-Krook collision model as depicted in the right hand side of the aforementioned equation and describes the relaxation of populations f i to their equilibrium state f i eq with relaxation time τ. For the advection-diffusion reaction equations of the reactive species a D2Q9 model is also implemented. Although the option of a D2Q5 model would be computationally more efficient, the D2Q9 lattice is preferred since the extra diagonal velocities offer a more physical description at the nodes where dissolution and precipitation take place.
The equilibrium populations ξ i eq for the reactive species i ξ are given: Experimental setup. The dissolution-precipitation experiments and post-mortem analysis and interpretation are described in detail in refs 9,16 , we therefore give only a summary of the work with details relevant for understanding the modeling. An acrylic flow cell of internal dimensions 10 cm × 10 cm × 1 cm was filled with a reactive porous medium sandwiched between two layers of an inert granular porous media as depicted in Fig. 5. The inert granular medium was a commercial acid washed and calcined quartz sand (SiO 2 ) of 99.1% purity. The reactive porous medium was created of naturally occurring celestine from Madagascar which comes in form of polished stones. The natural celestine was analyzed by X-ray diffraction which revealed a celestine content of 99.7% with 0.3% of anhydrous calcium sulfate. The stones were crushed and sieved to give batches of different grain size. The reactive medium was a mixture of grains: 30 wt. % with a size of less than 63 µm and 70 wt. % with a size of 125-400 µm. After compaction of the mixture a porosity of 0.33 ± 0.01 was measured. A stationary flow field was established by injection of a solution saturated with strontium sulphate to prevent strong dissolution of Celestine. This initial flow field, as well as flow field changes during the subsequent reaction experiments were visualized by injection of dye tracers.
For the reactive experiment a solution of 0.3 M barium chloride was injected at a constant flow rate of 20 µL min −1 , The initial strontium sulfate saturated solution has a lower density as the injected barium chloride solution; therefore at the start of the experiment the density differences enforced a transient flow and transport regime while the injected solution moved along the bottom of the tank. As the barium chloride solution reached the celestine region, the dissolution of celestine and the precipitation of baryte occurred, inducing changes in the pore space topology. We observed a strong decrease in permeability of the reactive layer. These changes were reflected in a rising level of barium chloride on the injection side of the tank which causes the ingress of barium chloride into increasingly higher parts of the reactive layer. The chemical composition of the effluents was measured (Cl −1 , SO 4 2− , Ba 2+ and Sr 2+ ) using ion chromatography. At the end of each experiment, but before dismantling the flow cell, isopropanol was injected in order to stop further chemical reactions. The reacted celestine was found to be solidified into a rectangular block. The reacted medium was dried and impregnated under vacuum with araldite XW936/XW397. An extensive post-mortem analysis of the reacted celestine layer was performed with different techniques including scanning electron microscopy, and synchrotron based X-ray micro-diffraction/micro fluorescence performed at the XAS beamline (Swiss light source).
The governing reaction describes the dissolution of solid celestine crystals and the precipitation of solid baryte. Baryte precipitates via two distinct mechanisms producing nano-crystalline (HON) and epitaxially grown rims (HET):