The dynamic nature of crystal growth in pores

The kinetics of crystal growth in porous media controls a variety of natural processes such as ore genesis and crystallization induced fracturing that can trigger earthquakes and weathering, as well as, sequestration of CO2 and toxic metals into geological formations. Progress on understanding those processes has been limited by experimental difficulties of dynamically studying the reactive surface area and permeability during pore occlusion. Here, we show that these variables cause a time-dependency of barite growth rates in microporous silica. The rate is approximately constant and similar to that observed on free surfaces if fast flow velocities predominate and if the time-dependent reactive surface area is accounted for. As the narrower flow paths clog, local flow velocities decrease, which causes the progressive slowing of growth rates. We conclude that mineral growth in a microporous media can be estimated based on free surface studies when a) the growth rate is normalized to the time-dependent surface area of the growing crystals, and b) the local flow velocities are above the limit at which growth is transport-limited. Accounting for the dynamic relation between microstructure, flow velocity and growth rate is shown to be crucial towards understanding and predicting precipitation in porous rocks.


Segmentation
Segmentation of the initial SiO2 porous structure from a 16-bit stack of images corresponding to the dry column before the experiment: 1) "Anisotropic diffusion 3D filter" (Threshold = 3000, Iterations = 10); 2) Threshold the fraction corresponding to the SiO2; 3) Mask the space outside the edge of the column. 4) "Fill holes 3D" function to remove air bubbles inside the beads.
Segmentation of the total volume of barite for each scan from 16-bit stacks of images: 1) Align the stack with the initial data set (fine registration function); 2) Mask the data with the initial SiO2 structure; 3) Threshold the brightest fraction that corresponds to barite. Particles of only one pixel were considered noise, thus not accounted as barite. This segmentation method can cause an overestimation of the volume if pixels that are accounted as barite are only partially filled with crystals, or an underestimation of the volume if the intensity of pixels that contain a small amount of crystals are below the applied threshold.
Calculation of Apore(t) for each scan: 1) five pixel dilation of the solid phase composed by silica and barite, 2) five pixel erosion of the result from 1, 3) overlap the result from 2 with a thresholded image of just barite crystals in each scan. The number of pixels of the overlapped image is multiplied by the surface area of the face of a voxel (i.e. a square), 1.54 m 2 . This sequence of operations identifies the pore-throats that are filled with crystals (blue in Supplementary Figure S1). The resulting area corresponds statistically to the surface area of pores that are occluded, thus no longer contribute to the growth rate. Finally, this surface area is subtracted from the initial Apore. Note that due to the lack of resolution it is not possible to account for the decrease of the reactive surface area caused by the lateral overlapping of two crystals. It is important to highlight that this is a statistical method aimed to minimize the error associated with overgrowth along restricted directions when applying equation (3) (see Supplementary Figure S1).  Figure 2: A distance map from the surface of the pore structure was generated using a 3-4-5 chamfer transform in Avizo 8, using the initial pore structure.

Graphic in
A threshold was applied to the color map in order to select the layer closest to the surface (LCS) of the beads that is 12.4 m thick. The resulting binary stack was multiplied by the stack of the flow field at t = 0 hours. The results consist of the flow velocities in the LCS, further called dataset A. The stack's histogram was divided in 100 groups by color intensity, each corresponding to an interval of fluid velocities (vx).
The total number of pixels in each group corresponds to ( ) in equation (S1).
The color scale was normalized using the maximum flow velocity corresponding to the original flow field (Table 1) if the overall Rpore is higher than Rfree, at 13.5 hours the volume precipitated is significantly lower than the amount expected on a free surface.
Supplementary Figure S5: Variation of the total volume of barite as a function of time.
Dashed line (X) corresponds to values measured in this experiment directly from sCT.
Continuous line correspond to the volume expected on free surfaces for the same solution composition using a nucleation density N=0.29 crystals/m 2 , the initial pore surface area S8 Apore=71.3 mm 2 , and an empirical equation to calculate the volume of single crystals V(t). 33 Dotted lines (o) correspond to volumes back calculated from Rpore and using Apore. Stage 2 when Rfree and Rcrystal are similar is zoomed in as an inset. Data corresponds to the total pore volume within a field of view, 3 mm 3 .

Methods
The mixing chamber was connected to the column by 15 cm of tygon tube 3 mm diameter. This is the minimum length technically possible that minimize precipitation before the solution enters the column, whilst allowing 360 degrees rotation of the sample.
BaCl2 and Na2SO4 solutions were prepared by dilution of 0.01 M stock solutions. Stock solutions were prepared from p.a. anhydrous Na2SO4 (Acros) and laboratory grade anhydrous BaCl2 (Fischer Scientific).
= 4 (S2) Supplementary Figure S6: Volume corresponding to each pore diameter at the beginning of the experiment measured by the method of inscribed spheres (function "thickness" in boneJ) 45 .

S9
Numerical simulation of flow velocities The Stokes problem was iteratively solved using the finite-difference scheme with pressure discretized at the center of the voxels and velocities at the boundaries of the voxels. This is known to produce accurate results while having low demand for computational resources. 46 An in-house C++ parallelized Stokes solver based on artificial compressibility method was used for all computations. 47 The code was tested on sphere packings and produce an accuracy comparable with D3Q19 TRT lattice-Botzmann simulations. 48 Model solutions only deviated significantly from analytical solutions for low resolved sphere packings where flow was controlled by connections between voxel edges, 49 which was not the case for all stages of our experiment (see pore-throat size distributions in Figure S5). All simulations were iterated until criteria based on imbalance for both continuity 43 (<10 -6 ) and motion equations have reached a plateau (less than 0.001% change in last 20 iterations). This required from 13 x 10 3 to 12 x 10 4 iterations depending on the subvolume, and ensured that the computerspecific limit, beyond which the difference between two consecutive iterations is small and random, was reached. After convergence of the velocity and pressure fields, the permeability was calculated using the Darcy's Law. 38 To study the correlation between fluid flow velocity and barite precipitation, permeability values are insufficient and the full velocity field is required. To calculate velocity fields representative of those in the sample at different times, the pressure boundary conditions were calculated using the permeability and initial flow rate to solve Darcy's law. First, equation (5) was solved for all five subvolumes to obtain permeability and velocity fields. 50,51 The permeability is a property of the pore geometry and the formulation of the Stokes problem depends only on the voxel resolution. The latter means that correct permeability values are obtained for all subvolumes with these pressure boundary conditions, which are different from experimental ones. Second, using the measured flow rate of 60 ml/h the unknown pressure drop conditions in the initial subvolume (t = 0 hours) were calculated. The continuity condition was applied, thus it was assumed that the flow rate is constant throughout the column. Finally, the velocity fields in the five subvolumes were rescaled to correct pressure drop in the initial subvolume.

Calculation of the Reynolds and Peclet numbers
To characterize transport conditions throughout the porous network, the Peclet number (Pe) and Reynolds number (Re) where calculated using equations (S3) and (S4). Pe is the ratio between diffusive and advective time scales, where D is the molecular diffusion coefficient of barium in solution, Uv is the Darcian velocity,  the total porosity measured experimentally by sCT and L is the characteristic length. For packings of mono and polydisperse spheres (our sample is close to such packings) it is common to use so-called Sauter diameter, 52 calculated from 6 ∑ / ∑ , where Vp is the volume of each particle and Ap is the particle surface area. In our calculation we did not divide the solid structure into separate particles (e.g. using a watershed algorithm), but rather used the total volume and total surface area of the solids measured from sCT. Calculated this way, the characteristic length is similar to that calculated in previous work for non-granular solid materials. 44