A semi-Lagrangian method for the direct numerical simulation of crystallization and precipitation at the pore scale

This article introduces a new efficient particle method for the numerical simulation of crystallization and precipitation at the pore scale of real rock geometries extracted by X-Ray tomography. It is based on the coupling between superficial velocity models of porous media, Lagrangian description of chemistry using Transition-State-Theory, involving underlying grids. Its ability to successfully compute dissolution process has been established in the past, and is presently generalized to precipitation and crystallization by means of adsorption modeling. Numerical simulations of mineral ${\text{CO}}_2$ trapping are provided, showing evidence of clogging/non-clogging regimes, and one of the main results is the introduction of a new non-dimensional number needed for this characterization.

1 Introduction

Studying reactive flows in porous media is essential to manage the geochemical effects arising from ${\text{CO}}_2$ capture and storage in natural underground reservoirs. While long-term predictions are commonly modeled at the field scale (Class et al., 2009), pore-scale approaches meanwhile provide insights into local geochemical interactions between the injected ${\text{CO}}_2$ and the aquifer structure (Payton et al., 2022). Through mathematical homogenization of the sub-micrometer porous medium and appropriate modeling, one can simulate the reactive processes that occur at the pore scale and predict their impact on the macro-scale properties, as developed in Allaire and Hutridurga (2012); Allaire et al. (2010). Geochemical processes are critical components for understanding the mineral trapping mechanisms and local evolving interfaces within the porous environment, and investigating the impact of such reactive processes provides insights into reservoir safety submitted to chemical interactions that may compromise the aquifer structure. Pore-scale modeling of reactive flow hence appears as a complementary mean to field scale studies wherein homogenization theory bridges the gap between these scales.

In this context, several geochemical mechanisms play a critical role in the ${\text{CO}}_2$ sequestration process and mainly involve precipitation, crystallization, and dissolution phenomena. On one side, carbonate precipitation and crystallization ensure efficient capture of the injected ${\text{CO}}_2$ in the form of minerals such as calcite, aragonite, or dolomite: this is referred to as mineral trapping, which informs about the storage capacities of the reservoir. These processes significantly impact the flow within the porous media at the pore scale, leading to restructuring of the flow path and morphological changes that alter, among other, the pore size distribution and the roughness of the interface due to partial or complete clogging of pore throats. Such alterations at the micro-scale subsequently alter the estimation of the macro-scale properties, namely, the porosity and permeability, and thereby require investigations to ensure wise management of the underground reservoir structures. On the other side, the reverse chemical process can also occur, resulting in carbonate mineral dissolution due to an acidification of the aqueous solution. This may compromise not only the efficiency of the trapping mechanisms, leading to an increase of both porosity and permeability, but also the integrity of the reservoir cap rock, and is, therefore, of great interest to prevent acute leakage issues. Consequently, one needs reliable estimations of the macro-properties changes due to these overall geochemical processes at the pore scale, to manage their impact on the reservoir scale modeling of ${\text{CO}}_2$ storage. This can be achieved through, first, efficient Direct Numerical Simulation (DNS) of reactive flows at the pore-scale and, subsequently, by embedding uncertainty concerns on the quantification of the petrophysical properties. In the present article, we address the first point with a focus on precipitation and crystallization modeling for ${\text{CO}}_2$ mineral storage into carbonate porous media.

Pore-scale investigation of reactive geochemical systems has garnered interest over the past decades based on imaging processes and laboratory experiments (Poonoosamy et al., 2023; Menke et al., 2015; Noiriel and Renard, 2022; Siena et al., 2021)), numerical simulations (Varzina et al., 2020; Patel et al., 2021; Payton et al., 2022; Soulaine et al., 2018), and a combination thereof (Molins et al., 2021; Noiriel and Soulaine, 2021). From this perspective, image-based DNS using microCT of a Representative Elementary Volume (REV) of a porous sample with efficient scientific computing and numerical method appears as a promising tool to query the impact of reactive processes on real rock geometries.

The present article focuses on the modeling aspects of ${\text{CO}}_2$ mineral trapping under the form of calcite crystal aggregates at the pore scale. Precipitation kinetics of calcite have been historically studied since the 1970s from the experimental and theoretical sides in Chou et al. (1989); Lasaga (1981); Plummer et al. (1978), and this has commonly established Transition State Theory (TST) as an efficient and straightforward way of predicting mineral reaction rate. Indeed, the deterministic TST is currently one of the most widely used models in reactive transport codes and DNS, detailed in Molins et al. (2012); Noiriel et al. (2021); Noiriel et al. (2016); Steefel et al. (2015). However, several doubts have risen in the research community about using such a deterministic approach for predicting crystal growth rates. In particular, comparison with experimentally determined growth rates has highlighted a wide range of discrepancies, querying the reliability of the TST model for overall crystallization processes, introduced in Hellevang et al. (2013); Pham et al. (2011). Meanwhile, probabilistic approaches, which find their origins in classical nucleation theory and the probabilistic nature of the precipitation and crystal growth mechanisms, have been developed in Masoudi et al. (2021); Nooraiepour et al. (2021a); Wolthers et al. (2012). These models make it possible to incorporate the effects of induction time characterizing the onset of the nucleation, the ionic affinities of the growing sites, and attachment frequencies of the ionic species involved in the reaction. Such attachment frequencies are, especially, significant for modeling surface adsorption and crystal aggregation that largely hinges on the surrounding porous structure in the sense that observable kinks and corners (above the voxel scale), for instance, are experimentally identified as preferential growing sites. However, such a geometrical dependency of the crystal aggregation is commonly neglected in most models, which makes it difficult to predict the spatial distribution of the new crystals.

Therefore, the main contribution of the present article is to build a novel numerical method dedicated to the precipitation and aggregation into crystal at the pore scale of microCT based geometry. Therefore, we formalize a two-step crystallization process wherein nuclei generation relies on a deterministic TST model before considering the probabilistic mineral aggregation—crystal growth—into the pore interface. The latter accounts for adsorption frequencies of the precipitate to the growth sites, which is weighted by a non-uniform probability of attachment rate depending on local mineral volume fraction. The current numerical method relies on a semi-Lagrangian approach, which handles a Lagrangian description of the chemistry with underlying grid methods for the hydrodynamic, based on the superficial velocity formalism introduced in the 1980s in Quintard and Whitaker (1988). The latter makes it possible to account for the involvement of the porous matrix in the overall flow process through a micro-continuum description of the medium. In this sense, one considers an intermediate state between the full resolution of each individual solid grain and the completely averaged continuum representation of the porous media at Darcy’s scale. This establishes two-scale models that are widely used in hydrodynamics pore-scale modeling and $\mu$CT image-based DNS, as shown in Molins et al. (2021); Panga et al. (2005); Soulaine et al. (2017). The present semi-Lagrangian formalism has been successfully employed in the context of carbonate dissolution at the pore scale (Etancelin et al., 2020) and extensively benchmarked against state-of-the-art numerical alternatives, as detailed in Molins et al. (2021). The present work originality is to blend together the micro-continuum hydrodynamics and the two-step process involving precipitation with TST and aggregation with attachment probability. The result is an efficient numerical tool able to adress clogging studies in real geometries.

Moreover, this article investigates a case study involving a pure calcite sample in which alcalin solution saturated with calcium and dissolved ${\text{CO}}_2$ flow through, applying the numerical method developed here and allowing to compute whether this leads to material clogging or not. We show that the cohabitation of precipitation and aggregation define naturally two kinds of Damköhler numbers ${\mathrm{D}\mathrm{a}}_{\text{II}}$. This case study exhibits that as long as the precipitation Damköhler number is high enough, the crystallization Damköhler number drives whether there is clogging or not, for several Peclet numbers.

The present manuscript is structured as follows. The section 2 introduces reactive flow model at the pore-scale of rocks. First, the Darcy-Brinkman-Stokes model including to the Kozeny-Carman correlation terms in order to model highly heterogeneous medium at its pore scale, naturally meaningful for evolving fluid-solid interface defining the pores structure. Second, the transport-reaction-diffusion of chemical species. In Section 3, we detail our particle method for highly heterogeneous diffusion arising from Archie’s law in the two-scale description of the medium. The application of such models to the precipitation and the crystallization for the ${\text{CO}}_2$ mineral trapping is described in Section 4. It includes the implementation and its HPC features on hybrid architectures, i.e., coupling CPU and GPU devices, for the present application, detailed in Section 4.4. The related numerical results are discussed in Section 5 in terms of clogging or non-clogging regimes of crystallization.

2 Models in reactive microfluidics

The present section focuses on the modeling of reactive hydrodynamics in the context of ${\text{CO}}_2$ mineral storage and presents the mathematical model used to simulate reactive processes at the pore scale. We first introduce the so-called Darcy-Brinkman-Stokes formulation for microfluidic flows based on superficial velocity formalism. We subsequently incorporate transport-reaction-diffusion equations modeling the geochemical interactions between the different species involved. Finally, we present an alternative formulation in velocity-vorticity for the hydrodynamics equation, which ensures the fluid incompressibility condition.

2.1 Darcy-Brinkman-Stokes: a superficial velocity formalism at the pore scale

We introduce a spatial domain $\mathrm{\Omega }\subset {\mathbb{R}}^n$, $n=\mathrm{1,2,3}$ which corresponds to the porous medium described at its pore scale. This sample description involves a pure fluid region ${\mathrm{\Omega }}_F$, also called void-space and assumed to be a smooth connected open set, and a surrounding solid matrix ${\mathrm{\Omega }}_S$ itself considered as a porous region. This region is seen as complementing the full domain $\mathrm{\Omega }$, which in practice represents the computational box of the numerical simulations such that ${\mathrm{\Omega }}_F=\mathrm{\Omega }\backslash {\bar{\mathrm{\Omega }}}_S$, and the internal fluid/solid interface is denoted $\mathrm{\Sigma }$. We denote the computational domain boundary by $\partial \mathrm{\Omega }$ and use ${\mathrm{\Gamma }}_F=\partial \mathrm{\Omega }\cap {\mathrm{\Omega }}_F$ and ${\mathrm{\Gamma }}_S=\partial \mathrm{\Omega }\cap {\mathrm{\Omega }}_S$ to refer to the fluid and solid parts of the computational domain boundary, respectively, such that $\partial \mathrm{\Omega }={\mathrm{\Gamma }}_F\cup {\mathrm{\Gamma }}_S$ (see Figure 1 for instance).

Figure 1

Figure 1. Upscaling from the pore-scale. Schematic representation of a reservoir or material scale structure, on the left, with its inherent averaged macro-properties $\varphi$ and $K$ computed on a representative elementary volume (REV). Local micro-continuum description of the pore-scale heterogeneity in this REV, on the right, along with its intrinsic micro-scale properties. These properties are the local micro porosity field $\epsilon$ and the micro-scale permeability $K_{\epsilon }$, based on the Kozeny-Carman relationship developed in Equation 2.

The boundary conditions at the inlet and outlet faces, typically for a cubic computational domain $\mathrm{\Omega }$ =]$0,l$[³ but not exclusively, either impose a prescribed flow rate $\bar{u}$ on the velocity or satisfy periodic boundary conditions for a prescribed driving force $f$. The boundary conditions on the other lateral faces are systematically periodic since rock samples are commonly constrained in an impermeable solid casing when $\mu$CT experiments are conducted. This ensures a consistent numerical representation of the sample compared to the experiments. This also guarantees ${\mathcal{C}}^{\infty }$ regularity on the boundary even if the domain exhibit corners, since the problem can be formalized by considering the equivalence relationship with the quotient space $\mathrm{\Omega }\equiv \mathcal{Q}/G$ where $\mathcal{Q}={\mathbb{R}}^2\times$]$0,l$[and $G=l{\mathbb{Z}}^2\times \{0\}$ (e.g., see Sanchez et al. (2019) for detailed configurations of acceptable domains).

From the $\mu$CT images, we can also characterize the static pore-space structure, corresponding to the sample’s initial state before any geochemical interactions. we denote by $\epsilon ={\epsilon }_f=1-{\epsilon }_s$ the micro-porosity field defined on $\mathrm{\Omega }$, given ${\epsilon }_f$ and ${\epsilon }_s$ respectively the volume fractions of void and solid according to usual notations from Soulaine et al. (2017). This defines a micro continuum description of the porous medium such that $\epsilon =1$ in the pure fluid region ${\mathrm{\Omega }}_F$ and takes a small value in the surrounding matrix ${\mathrm{\Omega }}_S$. In fact, the local micro-porosity $\epsilon$ is assumed to have a strictly positive lower bound $\epsilon (x,t)⩾{\epsilon }_0>0$ for all $(x,t)$ in the spatiotemporal domain $\mathrm{\Omega }\times [0,T_f]$ for a final real-time $T_f>0$ in the reactive process. This lower bound ${\epsilon }_0$ characterizes the residual porosity of the porous matrix, potentially unresolved due to X-ray $\mu$CT imaging limitations (as discussed in Perez et al. (2022), see also Figure 1). In practice, we assume throughout this work ${\epsilon }_0=5\%$, but this value can range from 1% to 10%.

Such a two-scale description of the local heterogeneities in the carbonate rocks is appropriate to simulate the pore-scale physics and establish the governing flow and transport equations in each distinct region. Indeed, although the hydrodynamic of a viscous flow in a pure fluid region is commonly quantified through the Navier-Stokes equation, we can formulate the problem on the whole domain $\mathrm{\Omega }$ based on the two-scale micro continuum description of the medium. We, therefore, consider the model on the superficial velocity $u$ introduced and derived rigorously by Quintard and Whitaker in the late 80 s (Quintard and Whitaker, 1988) and commonly used until nowadays (Molins et al., 2021; Soulaine et al., 2017; Wood et al., 2007):

$\begin{array}{ll}\hfill & {\epsilon }^{-1}\frac{\partial \rho u}{\partial t}+{\epsilon }^{-1}\mathrm{d}\mathrm{i}\mathrm{v}\left({\epsilon }^{-1}\rho u\otimes u\right)-{\epsilon }^{-1}\mathrm{d}\mathrm{i}\mathrm{v}\left(2\mu D\left(u\right)\right)+{\mu }^{*}{K}_{\epsilon }^{-1}u=f-\nabla p\hfill \end{array}(1)$

along with the divergence-free condition $\mathrm{d}\mathrm{i}\mathrm{v}u≔\nabla \cdot u=0$. It is noticeable that this incompressibility condition should be changed when considering evolving porous structures to account for density variations, especially in the context of fast dissolution or nucleation (Soulaine et al., 2018). Indeed, this only depicts that crystal nucleation within a liquid volume, for instance, drastically increases the density and induces divergence effects in its neighborhood. Nevertheless, the divergence-free condition can be assumed to remain if the characteristic time of fluid/solid interface changes is way larger than the hydrodynamics time scale (Soulaine et al., 2017), which is the case for our study. In Equation 1, the notation $D(u)$ refers the shear-rate tensor $D(u)=(\nabla u+\nabla u^T)/2$, $\mu$ is the dynamic viscosity, $p$ is the volumic pressure, $f$ the volumic driving force and $\rho$ the fluid density. The related viscosity $\mu *$ usually coincides with the fluid viscosity $\mu$ but may be different in order to account for viscous deviations.

The quantities $\rho$, $\mu$, $\mu *$ and $f$ are assumed to be constant. The quantities $\epsilon$, $\rho$, $\mu$ and $p$ are scalar fields; $u$, its rotational $\omega =\nabla \times u$, $f$ and $\nabla p$ are vector fields, while $D(u)$ and $K_{\epsilon }$ are matrices.

The permeability $K_{\epsilon }$ refers to the micro-scale permeability and depends on the local micro-porosity field $\epsilon$. This permeability of the micro-porous domain can be modeled by the empirical Kozeny-Carman relationship, historically introduced in Kozeny (1927) and Carman (1937), generalized in Quintard and Whitaker (1993) and now used meaningfully at the micrometer subscale (see Soulaine et al. (2017); Molins and Knabner (2019) and their next publications). This relation is given by

$K_{\epsilon }^{-1}={\kappa }_b^{-1}\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^3}(2)$

where ${\kappa }_b$ is the bulk permeability, which can be taken as a coarse estimation of the reference macro-scale permeability ${\kappa }_0$. For instance, Soulaine et al. (2017) estimated that four orders of magnitude below the permeability are required for ${\kappa }_b$ to ensure adherent boundary conditions at the pore interface. In this article, we consider both $K_{\epsilon }$, ${\kappa }_b$ and ${\kappa }_0$ as scalars, meaning we restrict ourselves to the isotropic case although this formalism can be extended to anisotropic porous media. Despite the validation of this correlation law for porous media at the sub-millimeter scale and its applicability to dissolution (Molins et al., 2021), we assume that it is also valid for crystal growth at pore walls as a reverse mechanism of dissolution. The dissolution-crystallization (whether it is one or two steps process) is not reversible since the resulting crystal may not reach the same tortuosity than the initial material (hence a different permability), but the similarity in the material is sufficiently strong to work with this assumption. The superficial velocity formulation Equation 1 defines a two-scale model that can be solved on the overall domain $\mathrm{\Omega }$ — using, for instance, penalization principles—and retrieves the usual Navier-Stokes equation in the pure fluid region ${\mathrm{\Omega }}_F$ (since $K_{\epsilon }^{-1}=0$ for $\epsilon =1$). Depending on the flow regime hypothesis, one can also encounter simplified versions of Equation 1 wherein some terms can be neglected. In the context of pore-scale simulations, in particular, the inertial effects become negligible compared to viscous forces due to low Reynolds number, denoted $\mathrm{R}\mathrm{e}$. The latter is a characteristic dimensionless number defined as:

$\mathrm{R}\mathrm{e}=\rho \bar{u}L/\mu ,$

where $\bar{u}$ and $L$ are respectively the characteristic velocity and length of the sample. The characteristic length $L$ can be related to average pore throat diameters and, therefore, we typically fall within the assumption $\mathrm{R}\mathrm{e}\ll 1$ throughout this manuscript.

At low Reynolds numbers and for highly viscous Darcian flows, Equation 1 hence reduces to the following Darcy-Brinkman-Stokes (DBS) model:

$-\mathrm{d}\mathrm{i}\mathrm{v}\left(2\mu D\left(u\right)\right)+\mu {\kappa }_b^{-1}\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^2}u=\epsilon \left(f-\nabla p\right),\text{in}\mathrm{\Omega }(3)$

where $\mu *=\mu$ for sake of readability. In the present work, we consider this DBS Equation 3, which is adequate in the flow regime hypothesis of low Reynolds number representative in pore-scale modeling. The DBS equation based on the superficial velocity is an efficient formalism to model the hydrodynamic in multi-scale porous media, and account for heterogeneous porosity levels.

2.2 Reactive flow model: general formulation

The DBS flow model Equation 3 needs to be complemented by transport-reaction-diffusion equations of the different species involved in the geochemical processes. These equations are derived from the mass balance of the chemical species (Soulaine et al., 2017), and can be written under the form:

$\frac{\partial \epsilon {\tilde{C}}_k}{\partial t}+\mathrm{d}\mathrm{i}\mathrm{v}\left(u{\tilde{C}}_k\right)-\mathrm{d}\mathrm{i}\mathrm{v}\left({\alpha }_k\left(\epsilon \right)\epsilon \nabla {\tilde{C}}_k\right)={\dot{m}}_k/M_k,(4)$

where ${\tilde{C}}_k={\rho }_f{\bar{\omega }}_{f,k}/M_k$ is a concentration per unit of fluid with $M_k$ the molar mass and ${\dot{m}}_k$ the rate of mass transfer for the $k^{\text{th}}$ species. We follow here the notations introduced by Quintard and Whitaker in Quintard and Whitaker (1988), and afterward commonly used, where ${\rho }_f$ is the fluid density and ${\bar{\omega }}_{f,k}$ is the mass fraction of the $k^{\text{th}}$ component averaged on the fluid phase. The term ${\alpha }_k(\epsilon )$ is a space-variable effective diffusion coefficient and accounts for a reduced diffusion in the surrounding porous matrix due to the tortuosity effect, which is usually quantified using Archie’s law, historically introduced in Archie (1942):

${\alpha }_k\left(\epsilon \right)=D_{m,k}{\epsilon }^{\eta }.$

In this empirical relationship, $\eta$ refers to the tortuosity index and $D_{m,k}$ to the molecular diffusion of the $k^{th}$ considered species, quantified at milimeter scale in Wakao and Smith (1962), generalized and downscaled in Coindreau and Vignoles (2005) and Glover (2016). We finally introduce $C_k=\epsilon {\tilde{C}}_k$, so that the Equation 4 is written:

$\frac{\partial C_k}{\partial t}+\mathrm{d}\mathrm{i}\mathrm{v}\left({\epsilon }^{-1}u{C}_k\right)-\mathrm{d}\mathrm{i}\mathrm{v}\left(D_{m,k}{\epsilon }^{1+\eta }\nabla \left({\epsilon }^{-1}{C}_k\right)\right)=R_k\left(\mathbf{C}\right),(5)$

which is no more than a superficial modeling of the chemistry, that is to say $C_k$ is the amount of moles per unit of volume while ${\tilde{C}}_k$ is the amount of moles per unit of fluid volume. The notation $R_k(\mathbf{C})$ refers to a function (without differential operators) that models the rate contribution of the chemical reactions for the $k^{\text{th}}$ component, where we denote by $\mathbf{C}\in {\mathbb{R}}^{N_s}$ the vector of the concentrations $C_k$ of all the $N_s$ chemical species. We distinguish $N_m$ mobile species, and $N_i$ immobile species such that $N_s=N_m+N_i$. The $k^{\text{th}}$ rate contribution $R_k(\mathbf{C})$ is, practically, the balance of kinetics of all reactions involving the $k^{\text{th}}$ species. The sign of individual reaction rates lies in the nature of the species $k$ considered, either positive for a chemical product or negative for a reactant.

The model Equation 5 is the formalism that we retain for the aqueous species in the liquid phase. In particular, this model highlights a superficial gradient operator denoted ${\nabla }^{\epsilon }≔\epsilon \nabla {\epsilon }^{-1}$ involved in the heterogeneous diffusion arising from the Archie’s law. One should notice that the superficial gradient can become highly sensitive at the mineral boundary, mainly due to jumps in the porosity levels on either side of the interface, and thus will require special considerations to adequately manage evolving medium under reactive processes.

Concerning the solid phase of concentration $C_{(s)}$ (e.g., the $k^{\text{th}}$ component in vector $\mathbf{C}$), which we assume contains only one chemical species of molar volume $\upsilon$ $(N_i=1)$, it is not subject to transport or diffusion, so that

$\frac{\partial C_{\left(s\right)}}{\partial t}=R_k\left(\mathbf{C}\right).(6)$

This solid concentration is subsequently linked to the micro-porosity $\epsilon$ by the relation $C_{(s)}=(1-\epsilon )/\upsilon$, so one gets

$\frac{\partial \epsilon }{\partial t}=-\upsilon R_k\left(\mathbf{C}\right).$

In the case of a typical reaction involving a unique solid $\text{X}\left(\text{s}\right)$ of molar volume $\upsilon$, and two aqueous species $\text{Y}$ and $\text{Z}$ in the liquid phase, with their respective positive stoechiometric coefficients ${\chi }_i$ and following, for instance, the general chemical reaction:

${\chi }_1{\text{X}}_{\left(s\right)}+{\chi }_2\text{Y}\rightleftharpoons {\chi }_3\text{Z},(7)$

we define the vector of concentrations $\mathbf{C}≔{(C_1,C_2,C_3)}^T={\left([{\text{X}}_{\left(\text{s}\right)}],[\text{Y}],[\text{Z}]\right)}^T\in {\mathbb{R}}^3$. Since there is only one reaction, one gets a unique kinetic balance written as $R_i(\mathbf{C})=\pm {\chi }_{i}R(\mathbf{C})$, with $R(\mathbf{C})$ the kinetic rate. By default, we assume a positive sign for the solid species, so that one follows the convention $R(\mathbf{C})<0$ for the forward reaction corresponding to the solid ${\text{X}}_{\left(\text{s}\right)}$ dissolution, while $R(\mathbf{C})>0$ for the reverse reaction, e.g., precipitation and crystallization processes. The sign for the aqueous species subsequently depends on its interaction with the solid $\text{X}\left(\text{s}\right)$: we get a positive sign for species $\text{Y}$, which is either consumed or produced in the same way as the solid, and a negative sign for species $\text{Z}$, which behaves oppositely. The reaction rate $R(\mathbf{C})$ can involve many concentrations, specific areas, chemical activities, equilibrium constants, etc. (see Section 4 thereafter for practical examples and further details).

Along with its boundary and initial conditions, the model for reaction (Equation 7) defines a set of partial differential equations modeling reactive flows at the pore scale:

$\left\{\begin{array}{c}-\mathrm{d}\mathrm{i}\mathrm{v}\left(2\mu D\left(u\right)\right)+\mu {\kappa }_b^{-1}\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^2}u=\epsilon \left(f-\nabla p\right),\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ \frac{\partial C_1}{\partial t}={\chi }_{1}R\left(\mathbf{C}\right),\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ \frac{\partial C_2}{\partial t}+\mathrm{d}\mathrm{i}\mathrm{v}\left({\epsilon }^{-1}u{C}_2\right)-\mathrm{d}\mathrm{i}\mathrm{v}\left(D_{m,2}{\epsilon }^{1+\eta }\nabla \left({\epsilon }^{-1}{C}_2\right)\right)={\chi }_{2}R\left(\mathbf{C}\right),\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ \frac{\partial C_3}{\partial t}+\mathrm{d}\mathrm{i}\mathrm{v}\left({\epsilon }^{-1}u{C}_3\right)-\mathrm{d}\mathrm{i}\mathrm{v}\left(D_{m,3}{\epsilon }^{1+\eta }\nabla \left({\epsilon }^{-1}{C}_3\right)\right)=-{\chi }_{3}R\left(\mathbf{C}\right),\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ \epsilon =1-\upsilon C_1,\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ +\text{adequate boundary and initial conditions, along with}\mathrm{d}\mathrm{i}\mathrm{v}u=0\hfill \end{array}\right.(8)$

which is strongly coupled, since $u$ and $\mathbf{C}$ depend on each other by means of the micro-porosity changes $\epsilon$ during the chemical process. It is also possible to straightforwardly substitute $C_1$ with $\epsilon$ in this system (Equation 8). Finally, one can notice that the reactive system is valid on the whole domain $\mathrm{\Omega }$, whether the local state is fluid or not. In the pure fluid region, this system indeed converges toward a Stokes hydrodynamic model coupled with a standard transport-diffusion equation. Mathematical modeling of reactive hydrodynamics at the pore-scale can be expressed under the general form of the PDE system (Equation 8) coupling DBS with transport-diffusion-reaction equations. It is noticeable that the micro-porosity $\epsilon$ remains in the range $[{\epsilon }_0,1]$ which provides a well-posed Darcy-Brinkman-Stokes equation for the flow due to the expressions of Kozeny-Carman term and reaction formula.

This can be extended naturally to systems of reactions involving as many aqueous species in the liquid phase as needed, and involving potentially several solids: in this case the porosity is a linear combination of solid species. Most of the configurations studied in this article involve solid calcite—or calcium carbonate—whose concentration is denoted $C_{{\text{CaCO}}_{3\left(\text{s}\right)}}$ or $[{\text{CaCO}}_3]$, and whose molar volume is given by $\upsilon =36.93\times 10^{-3}\mathrm{L}.\mathrm{m}\mathrm{o}{\mathrm{l}}^{-\mathrm{1}}$.

2.3 A velocity-vorticity formulation

Two distinct approaches are successfully used in the literature to solve numerically the DBS Equation 3, namely, the velocity-pressure or velocity-vorticity formulations (Hume and Poncet, 2021; Lamichhane, 2013; Molins et al., 2021; Angot, 2018). The latter, inherited from vortex methods (Chatelain et al., 2008; Cottet et al., 2000; El Ossmani and Poncet, 2010; Hejlesen et al., 2015) introduces the vorticity field $\omega$ which is intrinsically related to the fluid velocity $u$ through the relation:

$\omega =\nabla \times u.$

Several advantages arise when considering the velocity-vorticity formulation that regards the PDE unknowns $(u,\omega )$ and can be interpreted as describing the local spinning motions generated by the flow constraints. First of all, one can benefit from the velocity projection on divergence-free fields, and thereby analytically ensures the incompressibility condition $\mathrm{d}\mathrm{i}\mathrm{v}u=0$. Secondly, this formalism can be effectively coupled with splitting strategies that sequentially separate the resolution of distinct physical phenomena, such as convection and diffusion. Finally, this also makes it possible to eliminate the pressure unknown from the momentum equation by applying the curl operator on the DBS equation, which reads as follows:

$-\mu \mathrm{\Delta }\omega +\mu {\kappa }_b^{-1}\nabla \times \left(\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^2}u\right)=\nabla \epsilon \times \left(f-\nabla p\right)(9)$

given the assumption $\nabla \times f=0$. Developing the curl of the Kozeny-Carman term, one gets the following expression:

$\nabla \times \left(\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^2}u\right)=\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^2}\omega +2\left(\epsilon -1\right){\epsilon }^{-3}\nabla \epsilon \times u$

which, in practice, exhibits terms that become dominant compared to $\nabla \epsilon \times (f-\nabla p)$. Consequently, the right-hand side in the vorticity formulation of the DBS Equation 9 is usually neglected (Etancelin et al., 2020; Molins et al., 2021).

Equation 9 is then supplemented with an equation that retrieves the velocity field from the related vorticity, and results in the relation:

$-\mathrm{\Delta }u=\nabla \times \omega (10)$

Based on the incompressibility condition. In practice, the previous Poisson Equation 10 is not straightforwardly considered, and one relies on an alternative using a stream function $\psi :\mathrm{\Omega }\subset {\mathbb{R}}^3\to {\mathbb{R}}^3$ (a vector potential) solution of:

$\left\{\begin{array}{c}-\mathrm{\Delta }\psi =\omega ,\text{in}\mathrm{\Omega }\hfill \\ +\text{boundary}\text{conditions}\text{such}\text{that}\mathrm{d}\mathrm{i}\mathrm{v}\psi =0\text{\hspace{0.17em}on\hspace{0.17em}}\partial \mathrm{\Omega }.\hfill \end{array}\right.$

The condition $\mathrm{d}\mathrm{i}\mathrm{v}\psi =0$ on $\partial \mathrm{\Omega }$ is essential to ensure the overall incompressibility condition of the stream function on $\mathrm{\Omega }$ and thereby identify $u=\nabla \times \psi$. This requires satisfying appropriate boundary conditions, namely, the following combination of homogeneous Dirichlet/Neumann conditions for a computational cubic domain $\mathrm{\Omega }$ = ]$x_{\text{min}},x_{\text{max}}$[$\times$]$y_{\text{min}},y_{\text{max}}$[$\times$]$z_{\text{min}},z_{\text{max}}$[:

$\left\{\begin{array}{cc}x=x_{\text{min}}\text{\hspace{0.17em}or\hspace{0.17em}}x=x_{\text{max}}\hfill & :{\psi }_y={\psi }_z=\frac{\partial {\psi }_x}{\partial n}=0,\hfill \\ y=y_{\text{min}}\text{\hspace{0.17em}or\hspace{0.17em}}y=y_{\text{max}}\hfill & :{\psi }_x={\psi }_z=\frac{\partial {\psi }_y}{\partial n}=0,\hfill \\ z=z_{\text{min}}\text{\hspace{0.17em}or\hspace{0.17em}}z=z_{\text{max}}\hfill & :{\psi }_x={\psi }_y=\frac{\partial {\psi }_z}{\partial n}=0.\hfill \end{array}\right.$

Such boundary conditions ensure $(\nabla \times \psi )\cdot n=0$, $\mathrm{d}\mathrm{i}\mathrm{v}\psi =0$ and $\psi \times n=0$ at the same time, where $n$ is the normal field at the interface, and consequently lead to a zero average velocity field lifted by a prescribed flow rate $\bar{u}$ oriented in the flow direction. Indeed, given $u=\bar{u}+\nabla \times \psi$, one gets:

$<{u>}_{\mathrm{\Omega }}=\frac{1}{|\mathrm{\Omega }|}{\int }_{\mathrm{\Omega }}u\mathrm{d}\upsilon =\bar{u}+\frac{1}{|\mathrm{\Omega }|}{\int }_{\partial \mathrm{\Omega }}\psi \times n\mathrm{d}s=\bar{u},$

which replaces the setting of the driving force $f$ by a prescribed flow rate, managed through the lifted vector $\bar{u}$.

Finally, using the stream function $\psi$ analytically ensures the divergence-free condition on the velocity as $\mathrm{d}\mathrm{i}\mathrm{v}u=\mathrm{d}\mathrm{i}\mathrm{v}(\nabla \times \psi )=0$. Overall, the velocity-vorticity formulation (13) of the DBS equation is subsequently coupled with the transport-reaction-diffusion PDE system developed in Section 2.2 to model reactive hydrodynamics at the pore-scale.

3 Method: hybrid grid-particle scheme (semi-Lagrangian)

The present work relies on a semi-Lagrangian numerical method, mixing Eulerian and Lagrangian formalism, to tackle dynamically evolving porous media due to reactive micrometric processes. Such a semi-Lagrangian approach has been successfully used for simulations of calcite dissolution at the pore scale in Etancelin et al. (2020) and extensively benchmarked against state-of-the-art numerical alternatives in Molins et al. (2021).

A Lagrangian formalism consists of describing the flow motion through the observation along time of a large number of fluid particles, with their attached intrinsic properties and spatially varying positions (Cottet and Koumoutsakos, 2000). Each particle is then tracked throughout the evolving mechanism (transport, diffusion, … ) to measure variations in the main properties (velocity, concentration, … ). On the contrary, from the Eulerian point of view, the previous property changes are characterized on a predetermined spatial grid along the dynamical process. This section is dedicated to presenting this hybrid formalism, which is subsequently improved to account for the heterogeneous diffusion of the chemical reactants through the porous matrix.

The particle formulation is especially well-suited for transport-dominant phenomena as it avoids the explicit discretization of convective terms and alleviates the consideration of their related stability constraints—namely, the CFL conditions which constrains the time step for a given spatial discretization. The lack of regularity in the particle distributions throughout the dynamic process is, however, a recurring problem of Lagrangian methods. Indeed, as the particle positions move according to the flow field gradients, accumulation or scarcity issues in the particle distribution commonly occur. This, thereby, requires periodic remeshing steps to avoid this problem and not to lose information: namely, one proceeds two successive interpolations from the disorganized particle structure to a regular grid and subsequently from the grid to the new particle distribution (Cottet and Poncet, 2004; Magni and Cottet, 2012).

This is particularly suitable for hybrid approaches, wherein dedicated solvers can be straightforwardly implemented in the Eulerian context before performing the remeshing step. This also allows a representation of the quantities of interest on the grid, which can be coupled with domain decomposition or mesh adaptation methods. Hybrid grid-particle formalism has, thereby, garnered considerable interest in addressing multiple complex phenomena in Computational fluid Dynamics (CFD) and geosciences (Beaugendre et al., 2012; Chatelain et al., 2008; Chatelin and Poncet, 2013; El Ossmani and Poncet, 2010; Hume and Poncet, 2021). Besides, incorporating high-order and compact support interpolation kernels makes it possible to reduce the overall computational complexity of the remeshing steps while keeping accurate predictions of the numerical solution. The choice of the interpolation kernels is, however, important to ensure a robust numerical method and guarantee properties such as mass conservation and sign-preservation of the interpolated quantities (Magni and Cottet, 2012). Improvements of the interpolation kernels, especially for applications to dissolution processes at the pore scale, have been investigated by Etancelin et al. (2020). Such improvements focused on sign preservation and accurate high-order interpolation through a correction step of the potential over-diffusive effects resulting from the remeshing step. This provides a well-established hybrid grid-particle framework that can robustly address pore-scale reactive flows.

In the present work, we aim to benefit from the main advantages of both approaches to model reactive hydrodynamics at the pore scale. We will, thereby, use a Lagrangian description for the chemical equations—including the heterogeneous diffusion operator and reactant transport—with an underlying regular grid for solving the DBS equation in its velocity-vorticity formulation. To do so, we describe in this sections this original numerical scheme, by detailing its components and how they are linked together:

•Section 3.1 describes the Lagrangian formalism of reaction-diffusion-transport equations and its resulting dynamical system (the particle formalism): it requires a velocity field to transport particles and the computation of a heterogeneous diffusion,

•Section 3.2 shows how this velocity, solution to the Darcy-Brinkman-Stokes, can be computed by an operator-splitting strategy using the algorithm 1/Brinkman penalization 2/Diffusion using improved PSE 3/Projection on divergence-free fields,

•Section 3.3 details the computational method for the required heterogeneous diffusion (that is to say space and time variable diffusion, including the space mapping based on the porosity field), and shows how its original formalism from Degond and Mas-Gallic (1989) is improved to become intrinsically second order.

3.1 Reactive dynamical system using Lagrangian formulation

In this section, we present the Lagrangian formulation dedicated to the resolution of the reactive dynamical system (Equation 8). detailed in Section 2.2, and more specifically to the transport-reaction-diffusion Equation 5.

We define a set of $N_p$ fluid particles covering the computational domain $\mathrm{\Omega }$ and characterize as triplets ${(C_i,x_i,{\upsilon }_i)}_{i=1..N_p}$ of species concentrations $C_{i,k}$ ($k=1..N_m$ indexing the $N_m$ chemical species), positions $x_i\in \mathrm{\Omega }$ and volumes ${\upsilon }_i$, where $i$ refers to the particle index. This mathematically introduces the particle description, denoted $C^h$, of the concentration fields as follows:

$C^h\left(t\right)={\displaystyle \sum _{i=1}^{N_p}}{C}_i\left(t\right){\upsilon }_i\left(t\right){\delta }_{x_i}\left(t\right)(11)$

which only depends on time $t$, and where $\delta$ refers to the Dirac function. The Lagrangian formulation of Equation 5 can then be written using the particle description Equation 11:

$\left\{\begin{array}{c}\frac{\mathrm{d}{C}_{i,k}}{\mathrm{d}t}=R_k\left({\mathbf{C}}_i\left(t\right)\right)+{\left[\mathrm{d}\mathrm{i}\mathrm{v}\left({\alpha }_k\left(\epsilon \right){\nabla }^{\epsilon }{C}_k\right)\right]}_{x_i\left(t\right)}\forall k=1..N_m\hfill \\ \frac{\mathrm{d}{x}_i}{\mathrm{d}t}={\left[{\epsilon }^{-1}u\right]}_{x_i\left(t\right)}\hfill \\ \frac{\mathrm{d}{\upsilon }_i}{\mathrm{d}t}=0\hfill \end{array}\right.$

given the incompressibility condition $\mathrm{d}\mathrm{i}\mathrm{v}u=0$ and the notations introduced in Section 2.2. This results in a dynamical system over the particles whose positions are controlled by the field ${\epsilon }^{-1}u$, and volumes remain constant under divergence-free conditions. The main advantage of such a Lagrangian formulation (20) is that the transport term $\mathrm{d}\mathrm{i}\mathrm{v}({\epsilon }^{-1}u{C}_k)$ vanishes along with its stability condition and, thereby, the method presents the ability to use arbitrary large time steps. This is, especially, of great interest when the Courant–Friedrichs–Lewy (CFL) condition on the transport term induces a stronger constraint on the time step compared to the diffusion stability condition.

The velocity field $u$ in (20) arises from the solution of the DBS equation which is solved on an underlying Cartesian grid and coupled with the Lagrangian formulation of the chemical PDE system. Regarding such a strong coupling between these equations, one needs to interpolate on the grid the particle description of the solid chemical species—namely, $C_{{\text{CaCO}}_{3\left(\text{s}\right)}}^h$ — which is related to the micro-porosity field $\epsilon$ and consequently involved in the DBS model. Similarly, the velocity field subsequently needs to be interpolated on the particles to solve the Lagrangian set of chemical equations. This requires a convolution with high-order remeshing kernels with compact supports (Etancelin et al., 2020; Magni and Cottet, 2012). The dynamical system (20) is finally integrated using standard numerical methods for Ordinary Differential Equations (ODE), such as explicit Runge-Kutta, while the diffusion term $\mathrm{d}\mathrm{i}\mathrm{v}({\alpha }_k(\epsilon ){\nabla }^{\epsilon }{C}_k)$ is approximated through the Particle-Strength- Exchange (PSE) method, described in Section 3.3 and detailed in Supplementary Appendix.

In the present work, we incorporate in the semi-Lagrangian workflow the consideration of robustly estimating Archie’s law term through such a particle-based PSE method. Thanks to this method, we make it possible to fully address the superficial gradient ${\nabla }^{\epsilon }$ approximation with heterogeneous diffusion throughout the porous matrix.

3.2 Splitting operator strategy

The semi-Lagrangian formalism introduced in the previous Section 3.1 intrinsically relies on splitting strategies. Time-splitting methods, also known as fractional time-step algorithms, arise in many fields of computational science related to physics-based models and have been developed by Chorin (1973) in the context of vortex methods for the Navier-Stokes equation. Such methods have been widely investigated afterward to separate the resolution of distinct physical phenomena and render more efficient algorithms (Beale and Majda, 1981; Hume and Poncet, 2021; Faragó, 2008). Indeed, one of the main advantages of splitting strategies is one can use different approaches for the distinct parts of the overall model, namely, either a Lagrangian or Eulerian formulation. This also straightforwardly extends to the choice of the numerical solver available for each component, allowing the use of the most efficient, accurate, and robust schemes independently.

The first natural splitting arising in workflow, thereby, lies in the semi-Lagrangian formalism itself, wherein we do not consider solving the overall Partial Differential Equation (PDE) system at once. We rather separate the transport-diffusion-reaction dynamics in its Lagrangian formulation from the pore-scale hydrodynamic resolved on the underlying Cartesian grid. The hydrodynamic part, composed of the DBS equation in the velocity-vorticity formulation, is also solved through a time-splitting method. Indeed, we approximate the solution of Equation 9 by the limit in time of the evolution equation

$\frac{\partial \omega }{\partial t}-\mu \mathrm{\Delta }\omega +\mu {\kappa }_b^{-1}\nabla \times \left(\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^2}u\right)=0,$

together with $\omega =\nabla \times u$, using a three-step operator splitting strategy coupled with a fixed-point algorithm. This means that considering a sequence $(u_m,{\omega }_m)$ of velocity-vorticity we aim to successively perform, over a time interval $[t_m,t_{m+1}]$ with $t_m=m\delta t$, Brinkman penalization, diffusion, and projection on divergence-free fields. The latter is achieved through the reconstruction of the velocity field $u$ based on the stream function $\psi$ (see Section 2.3). In practice, these three steps are specifically defined as:

•The Brinkman iteration given by the ordinary differential equation $\frac{\partial u}{\partial t}+\mu \lambda (\epsilon )(u+\bar{u})=0$ with prescribed flow rate $\bar{u}$ and $\lambda (\epsilon )≔{\kappa }_b^{-1}{(1-\epsilon )}^2{\epsilon }^{-2}$, whose exact solution after a $\delta t$ is generated by

$\mathrm{\Lambda }\left(u\right)≔e^{-\mu \lambda \left(\epsilon \right)\delta t}\left(u+\bar{u}\right)-\bar{u}$

•The vorticity diffusion iteration, $\frac{\partial \omega }{\partial t}-\mu \mathrm{\Delta }\omega =0$, solved using an implicit Euler scheme given by the operator

${\mathcal{D}}_{\omega }\left(u\right)≔{\left[I-\mu \delta t\mathrm{\Delta }\right]}^{-1}\left(\nabla \times u\right)$

•The projection step $\mathrm{\Pi }(\zeta )=\nabla \times \left((-{\mathrm{\Delta }}^{-1})\zeta \right)$ which takes as $\zeta$ the right-hand side of the Poisson equation $-\mathrm{\Delta }\psi =\zeta$ satisfying the boundary conditions (17), and followed by $u=\nabla \times \psi$.

Overall, this leads to the full iteration of the Brinkman-Diffusion-Projection splitting $\mathrm{\Pi }◦{\mathcal{D}}_{\omega }◦\mathrm{\Lambda }$ over a time step $[t_m,t_{m+1}]$, which reads as follows:

$u_{m+1}=\mathrm{\Pi }◦{\mathcal{D}}_{\omega }◦\mathrm{\Lambda }\left(u_m\right)$

and whose consistency has been theoretically discussed in Hume and Poncet (2021). One should notice that this projection step is not a projection by pressure gradient correction as in Chatelin and Poncet (2013), but an operator that takes the vorticity field $\omega$ and retrieves a divergence-free velocity field whose mean velocity is zero. The final velocity field is subsequently given by

$u=\bar{u}+\underset{m\to \infty }{\lim }{u}_m$

whose average is $\bar{u}$ and satisfies the Kozeny-Carman relation inside the solid region.

From a numerical perspective, we consider the exact treatment of the Brinkman term, a fourth-order finite difference scheme for the discrete curl operator, and FFT solvers for the vorticity diffusion and stream function recovery. Using FFT avoids matrix assembly procedure and, therefore, consists of efficient solvers in terms of computational time and memory storage requirements. Besides, the complexity of such algorithm scales as $\mathcal{O}(N_p\log (N_p))$, where $N_p$ is recalled to refer to the number of particles. Finally, a stopping criterion on this fixed point algorithm is also defined based on the relative residual norm on the velocity, which manages the convergence of the pore-scale hydrodynamics toward a stationary state.

The updated velocity $u_{m+1}$ is subsequently interpolated from the grid to the particles and used for solving the Lagrangian reactive system, which is split into convective/remeshing and diffusive/reactive steps. Regarding the convection, the particle trajectories are pushed to the next step through a directional advection, given the field ${\epsilon }^{-1}u$ according to the Lagrangian formulation (20), and are then remeshed to avoid stagnation issues. The purpose of such a directional splitting is to reduce the dimensionality of the overall advection problem by considering several one-dimensional equations, and is developed in Cottet et al. (2014); Magni and Cottet (2012). The particle remeshing is also addressed using directional treatment of the multidimensional convolution. This means that within a time step $[t_m,t_{m+1}]$, the joint step of advection/remeshing of the particles is successively performed by alternating the spatial directions (Etancelin, 2014; Keck, 2019). This presents the advantage of significantly improving the computational efficiency of the method and allows the use of high-order remeshing kernels with large stencils while maintaining a minimal cost compared to multidimensional cases. This is also well-suited to parallel implementation on GPU architecture. The dimensional splitting is addressed, in practice, by a second-order Strang formula (Strang, 1968) and coupled with a second-order Runge-Kutta for time integration. The diffusion/reaction step is finally solved by means of a second-order explicit Runge-Kutta scheme along with PSE approximation of the heterogeneous diffusion operator. Once the Lagrangian formulation of the chemistry has been fully updated and remeshed on the underlying grid, one starts pushing the DBS hydrodynamics to the next sequential step of these temporal iterations.

Such an operator splitting strategy, in the context of a semi-Lagrangian approach, has been applied to the modeling of dissolution processes on a 2D synthetic calcite core and validated against state-of-the-art alternatives and experiments in Molins et al. (2021). This has also been used in Etancelin et al. (2020) on real porous media structures at the pore scale to investigate the dissolution of 3D carbonate rocks arising from $\mu$CT scans. Nonetheless, these previous works assumed that the superficial gradient ${\nabla }^{\epsilon }$ involved in the heterogeneous diffusion could be approximated by the gradient operator, and subsequently addressed this Archie’s law term with standard finite differences schemes. In Section 3.3, we intend to alleviate this assumption and, therefore, improve the semi-Lagrangian method by incorporating in the workflow a PSE approximation of the heterogeneous diffusion.

3.3 Particle-strength-exchange method and Archie’s law approximation

The PSE method consists in the approximation of a diffusion operator $\mathrm{d}\mathrm{i}\mathrm{v}(\mathbf{L}\nabla f)(x)$ with $x\in \mathrm{\Omega }\subset {\mathbb{R}}^d$ and $\mathbf{L}$ a positive symmetric matrix, accounting for heterogeneous diffusion for instance. The main idea is then to approximate the diffusion by an integral operator, more suitable for particle methods:

$Q^{\xi }\cdot f\left(x\right)={\int }_{\mathrm{\Omega }}{\sigma }^{\xi }\left(x,y\right)\left(f\left(y\right)-f\left(x\right)\right)dy$

where the kernel ${\sigma }^{\xi }$ is supposed to be symmetric and satisfies some moment conditions, described thereafter and detailed in Supplementary Appendix. In the Lagrangian formulation, a particle approximation of the function $f$, denoted $f^h$, is also introduced based on the particle triplet $(f_i,x_i,{\upsilon }_i)$, such that:

$f^h={\displaystyle \sum _{i=1}^{N_p}}{f}_i{\upsilon }_i{\delta }_{x_i}\text{where}{f}_i=f\left(x_i\right)$

where $x_i$ and ${\upsilon }_i$ are respectively the particle positions and volumes, while ${\delta }_{x_i}$ refers to the Dirac measure at position $x_i$. With such a $N_p$-particle representation of the function $f$, a discrete version of the operator $Q^{\xi }$ is obtained by using the particles as quadrature points where $h$ refers to the characteristic distance between particles. This results in the following quadrature expression, called Particle-Strength-Exchanges:

$Q^{\xi }\cdot f^h\left(x_k\right)={\displaystyle \sum _{x_l\in \mathcal{S}\left(x_k\right)}}{\sigma }^{\xi }\left(x_k,x_l\right)\left(f_l-f_k\right){\upsilon }_l.$

where $\mathcal{S}(x_k)≔\text{Supp}\left({\sigma }^{\xi }(x_k,\cdot )\right)$ refers to the set of points in the support of the kernel function ${\sigma }^{\xi }$.

In Supplementary Appendix, we detail that we can explicitly compute two constants ${\gamma }_1$ and ${\gamma }_2$ based on a given scale $\xi$ and a given function $\mathrm{\Theta }$ (named stencil generator) such as

${\sigma }^{\xi }\left(x,y\right)=\frac{1}{{\xi }^{d+4}}\mathrm{\Theta }\left(\frac{y-x}{\xi }\right)\frac{\mathbf{m}\left(x\right)+\mathbf{m}\left(y\right)}{2}:{\left(x-y\right)}^{\otimes 2}$

with

$\mathbf{m}=c_0\mathbf{L}-c_1\mathrm{T}\mathrm{r}\left(\mathbf{L}\right){\mathbf{I}\mathbf{d}}_{\mathbf{3}}+\mathbf{H}$

$\begin{array}{ll}\hfill c_0& =\frac{2\left({\gamma }_1+2{\gamma }_2\right)}{{\gamma }_1^2+{\gamma }_1{\gamma }_2-2{\gamma }_2^2},c_1=\frac{2{\gamma }_2}{{\gamma }_1^2+{\gamma }_1{\gamma }_2-2{\gamma }_2^2},\hfill \\ \hfill H_{ij}& =\left(\frac{{\gamma }_1^2-{\gamma }_1{\gamma }_2-6{\gamma }_2^2}{{\gamma }_2\left({\gamma }_1^2+{\gamma }_1{\gamma }_2+2{\gamma }_2^2\right)}\right)\left(1-{\delta }_{ij}\right)L_{ij}\hfill \end{array}$

In Supplementary Appendix, we detail that we can explicitly compute two constants ${\gamma }_1$ and ${\gamma }_2$ based on a given scale $\xi$ and a given function $\mathrm{\Theta }$ (named stencil generator) such as

${\sigma }^{\xi }\left(x,y\right)=\frac{1}{{\xi }^{d+4}}\mathrm{\Theta }\left(\frac{y-x}{\xi }\right)\frac{\mathbf{m}\left(x\right)+\mathbf{m}\left(y\right)}{2}:{\left(x-y\right)}^{\otimes 2}$

with

$\mathbf{m}=c_0\mathbf{L}-c_1\mathrm{T}\mathrm{r}\left(\mathbf{L}\right){\mathbf{I}\mathbf{d}}_{\mathbf{3}}+\mathbf{H}$

and

$\begin{array}{ll}\hfill c_0& =\frac{2\left({\gamma }_1+2{\gamma }_2\right)}{{\gamma }_1^2+{\gamma }_1{\gamma }_2-2{\gamma }_2^2},c_1=\frac{2{\gamma }_2}{{\gamma }_1^2+{\gamma }_1{\gamma }_2-2{\gamma }_2^2},\text{and}\hfill \\ \hfill H_{ij}& =\left(\frac{{\gamma }_1^2-{\gamma }_1{\gamma }_2-6{\gamma }_2^2}{{\gamma }_2\left({\gamma }_1^2+{\gamma }_1{\gamma }_2+2{\gamma }_2^2\right)}\right)\left(1-{\delta }_{ij}\right)L_{ij}\hfill \end{array}$

where ${\delta }_{ij}$ is the Kronecker symbol, ${\gamma }_1$, and ${\gamma }_2$ are the coefficients based on the fourth moments of $\mathrm{\Theta }$ defined in 3D by

$\begin{array}{ll}\hfill & {\gamma }_1={\displaystyle \sum _{x\in \mathbb{J}}}{x}_k^4\mathrm{\Theta }\left(x\right)h^3,k\in \left[\mathrm{1,3}\right]\hfill \\ \hfill & {\gamma }_2={\displaystyle \sum _{x\in \mathbb{J}}}{x}_k^2{x}_l^2\mathrm{\Theta }\left(x\right)h^3,k\ne l\in \left[\mathrm{1,3}\right]\hfill \end{array}$

for $\mathbb{J}\subset h{\mathbb{Z}}^3$ a three-dimensional lattice, including at least one neighborhood of the current mesh point. The convergence of the original PSE method introduced in Degond and Mas-Gallic (1989) depends on the scale $\xi$ as $\mathcal{O}[{(h/\xi )}^2]$, which becomes $\mathcal{O}(1)$ when $\xi$ is adapted linearly to $h$. The present discrete corrected version developed in Poncet (2006) and Schrader et al. (2010) is intrinsically second order convergence $\mathcal{O}(h^2)$ and is detailed in Supplementary Appendix.

This computational method of diffusion is used to compute the Archie’s law, involving a tortuosity index $\eta$ in the heterogeneous diffusion operator

$\mathcal{D}\left(\epsilon ,C\right)≔\mathrm{d}\mathrm{i}\mathrm{v}\left({\epsilon }^{1+\eta }\nabla \left({\epsilon }^{-1}C\right)\right)$

Involved in the PDE system (11). In this case, one gets a particular expression for the diffusion matrix $\mathbf{L}={\epsilon }^{1+\eta }{\mathbf{I}}_3$ with ${\mathbf{I}}_3$ the identity matrix in ${\mathbb{R}}^3$. This finally leads to the following discretized Archie’s law diffusion operator (30) using the renormalized PSE scheme:

$\begin{array}{ll}\hfill Q^{\xi }\cdot f^h\left(x_k\right)& =\frac{1}{{\xi }^7}{\displaystyle \sum _{l\sim k}}\left(f_l-f_k\right)\mathrm{\Theta }\left(\frac{x_l-x_k}{\xi }\right)\frac{\left(c_0-3c_1\right)}{2}\hfill \\ \hfill & \times \left({\epsilon }^{1+\eta }\left(x_k\right)+{\epsilon }^{1+\eta }\left(x_l\right)\right)|x_l-x_k{|}^2{\upsilon }_l\hfill \end{array}(12)$

where $|.|$ is the Euclidean norm in ${\mathbb{R}}^3$ and $f_{•}≔{\epsilon }^{-1}(x_{•})C(x_{•})$. The overall Formula 12 accounts for the heterogeneous diffusion and ensures the accurate management of the chemical reactant penetration, given by its concentration field $C$, within the porous matrix. Indeed, one of the main advantages of the PSE scheme is that it includes all the lattice neighborhoods in the computation of the heterogenous diffusion, unlike crossed finite differences scheme. Finally, this guarantees the strict conservation of the reactant exchanges between the fluid portion and the porous matrix, and is numerically validated at the second order in Section 4 of Supplementary Appendix.

4 Application to precipitation and crystallization

${\text{CO}}_2$ mineral storage in natural underground reservoirs, such as saline aquifers, involves competing geochemical phenomena occurring at a large variety of scales. Among them, mineral dissolution and precipitation play crucial roles. On one side, studying the dissolution of native carbonate species, already present in the aquifers, provides insight into potential leakage issues and queries the reservoir safety. On the other side, ${\text{CO}}_2$ trapping under the form of carbonate precipitates and crystals informs on the storage capacities of the reservoir. These geochemical processes also induce changes in the macro-scale properties of the subsurface material, including permeability and porosity evolutions, that need to be investigated to ensure sustainable management of the reservoir structures.

In this section, we develop mathematical models for calcite precipitation and crystallization at the pore scale, with special considerations on the reaction rate expressions arising in the PDE system (11). Numerical simulations are performed within the HySoP platform (Etancelin et al., 2022) along with PSE treatment of the heterogeneous diffusion on accelerated GPU devices and address porous sample arising from X-ray $\mu$CT observations. This enables the investigations of macro-scale property changes along the ${\text{CO}}_2$ mineral trapping on real 3D rock geometries, which is an important component in the overall study of ${\text{CO}}_2$ storage.

4.1 The transition state theory: from dissolution to precipitation modeling

Dissolution of the injected ${\text{CO}}_2$ in the aqueous phase of deep underground reservoirs will affect the pH of the formation water through the following series of chemical reactions:

${\mathrm{C}\mathrm{O}}_{2(\mathrm{g})}\rightleftharpoons {\mathrm{C}\mathrm{O}}_{2(\mathrm{a}\mathrm{q})}(\mathrm{13a})$

${\mathrm{H}}_2\mathrm{O}+{\mathrm{C}\mathrm{O}}_2(\mathrm{a}\mathrm{q})\rightleftharpoons {\mathrm{H}}_2{\mathrm{C}\mathrm{O}}_3(\mathrm{13b})$

${\mathrm{H}}_2{\mathrm{C}\mathrm{O}}_3\rightleftharpoons {\mathrm{H}}^{+}+{\text{HCO}}_3^{-}(\mathrm{13c})$

${\mathrm{H}\mathrm{C}\mathrm{O}}_3^{-}-\rightleftharpoons {\mathrm{H}}^{+}+{\mathrm{C}\mathrm{O}}_3^{\mathrm{2-}}.(\mathrm{13d})$

Indeed, once the ${\text{CO}}_2$ has dissolved into water and established a first equilibrium under the form of the weak acid ${\text{H}}_2{\text{CO}}_3$, the ${\text{H}}_2{\text{CO}}_3$ species dissociates successively to bicarbonate ${\text{HCO}}_3^{-}$ and carbonate ${\mathrm{C}\mathrm{O}}_3^{\mathrm{2-}}$ ions as the pH increases. These chemical reactions are pH-dependent, and the distribution evolutions of all these carbonate species are displayed in Figure 2 against the pH of the solution. In alkaline media, the chemical reactions (Equation 13c) and (Equation 13d) can, therefore, join together to read as follows:

${\text{H}}_2{\text{CO}}_3+2{\text{OH}}^{-}\rightleftharpoons {\text{CO}}_3^{2-}+2{\text{H}}_2\text{O}$

such that carbonate ions are the main carbonate species present in the solution. Such transformations in the ionic species composition of the aquifer water will considerably impact the original mineral structure through chemical rock-water interactions such as carbonate dissolution and precipitation.

Figure 2

Figure 2. Distribution diagram of aqueous carbonate species against pH solution from Bohn et al. (1980) and Turk et al. (2015). Species distributions are represented as a fraction of total dissolved carbonate. The grey dotted lines highlight the transition pH of the chemical equilibria from (Equation 13).

Historically, the dissolution and precipitation kinetics of calcite in the context of ${\text{CO}}_2$ injection have been studied since the 1970s, both from the experimental and theoretical sides. Plummer et al. (1978) investigated the influence of several parameters on the forward reaction rates of calcite dissolution under far-from-equilibrium conditions. Among these parameters, one finds the partial pressure of ${\text{CO}}_2$ denoted $P_{{\text{CO}}_2}$, the hydrogen ions activities denoted $a_{{\text{H}}^{+}}$ — directly related to the $\text{pH\hspace{0.17em}}$ — and the temperature. Their experimental work was subsequently extended by Chou et al. (1989) using a fluidized bed reactor to compare the dissolution kinetics mechanisms between various carbonate minerals—involving inter alia calcite, aragonite, and dolomite—at $25^{◦}\mathrm{C}$. One should notice that these experiments were conducted under laboratory conditions in terms of pressures and temperatures, in opposition to the current abilities of in situ experiments to manage realistic reservoir conditions (Andrew et al., 2013; Wigand et al., 2008). Nonetheless, these experimental studies have highlighted three kinetic mechanisms occurring simultaneously in the process of calcite dissolution due to ${\text{CO}}_2$ injection. Such mechanisms are given by the following chemical reactions:

${\text{CaCO}}_{3\left(\text{s}\right)}+{\text{H}}^{+}\stackrel{{\text{K}}_1}{\rightleftharpoons }{\text{Ca}}^{2+}+{\text{HCO}}_3^{-}(\mathrm{14a})$

${\text{CaCO}}_{3\left(\text{s}\right)}+{\text{H}}_2{\text{CO}}_3\stackrel{{\text{K}}_2}{\rightleftharpoons }{\text{Ca}}^{2+}+2{\text{HCO}}_3^{-}(\mathrm{14b})$

${\text{CaCO}}_{3\left(\text{s}\right)}\underset{{\text{K}}_{\mathrm{-3}}}{\overset{{\text{K}}_3}{\rightleftharpoons }}{\text{Ca}}^{2+}+{\text{CO}}_3^{2-}(\mathrm{14c})$

where the notations $K_i,i=1\dots 3$ refer to forward reaction rate constants, depending on the temperature (Busenberg and Plummer, 1986; Plummer et al., 1978; Plummer and Busenberg, 1982), and $K_{-3}$ is the backward reaction rate corresponding to the reverse calcite precipitation process in Equation 14c. They experimentally identified both the forward and backward reaction rates and established the validity of kinetic models for carbonate dissolution and precipitation in comparison to thermodynamics theoretical considerations. In the present study, we assume the production of pure calcite without any of its polymorphs like aragonite or vaterite.

Meanwhile, mineral reaction rates were, indeed, theoretically investigated by Lasaga (1981) using the Transition State Theory (TST), originally formulated by Eyring (1935). Since then, this formalism has successfully been extended (Aagaard and Helgeson, 1982; Lasaga, 1984; Steefel and Lasaga, 1994) and widely accepted in current kinetic geochemical models (Etancelin et al., 2020; Molins et al., 2012; Steefel et al., 2015). In this context, the reaction rates are commonly expressed as the product of far-from-equilibrium terms, involving the activities of the chemical species in solution, with an affinity term written as a function of the Gibbs free energy change $\mathrm{\Delta }G$ for close to the equilibrium conditions. Considering the chemical model of calcite dissolution (34) suggested by Plummer et al. (1978) and Chou et al. (1989) and the vector of concentrations $\mathbf{C}$, the reaction rate arising from TST writes:

$R\left(\mathbf{C}\right)=A_s\left(K_1{a}_{{\text{H}}^{+}}+K_2{a}_{{\text{H}}_2{\text{CO}}_3}+K_3\right)\left(\frac{a_{{\text{Ca}}^{2+}}{a}_{{\text{CO}}_3^{2-}}}{K_{\text{eq}}}-1\right)(15)$

where $K_{\text{eq}}$ is the equilibrium constant of the reaction, also called the solubility product, $A_s$ is the reactive surface area of the mineral—in ${\mathrm{m}}^{-\mathrm{1}}$. The notations $a_{•}={\gamma }_{•}{C}_{•}$ refer to the dimensionless species activities with ${\gamma }_{•}$ and $C_{•}$, respectively, their activity coefficients and molar concentrations—whose unit is $\mathrm{m}\mathrm{o}\mathrm{l}.{\mathrm{m}}^{-\mathrm{3}}$. It follows that the micro-porosity changes reads as:

$\frac{\partial \epsilon }{\partial t}=-\upsilon R\left(\mathbf{C}\right),$

given Equation 6 from Section 2.2 and the relation $C_{{\text{CaCO}}_{3\left(\text{s}\right)}}=(1-\epsilon )/\upsilon$. Denoting by $Q=a_{{\text{Ca}}^{2+}}{a}_{{{\text{CO}}_3}^{2-}}$ the ion activity product, one obtains the following relation between $Q$ and the Gibbs energy change (Aagaard and Helgeson, 1982; Lasaga, 1984):

$\mathrm{\Delta }G=RT\ln \left(\frac{Q}{K_{\text{eq}}}\right)$

with $T$ the temperature in Kelvin $\mathrm{K}$, and $R$ the universal gas constant in $\mathrm{J}.\mathrm{m}\mathrm{o}{\mathrm{l}}^{-\mathrm{1}}.{\mathrm{K}}^{-\mathrm{1}}$. The sign of the reaction term $R(\mathbf{C})$ in (Equation 15) is driven by the sign of $\ln \left(Q/K_{\text{eq}}\right)$ that is negative for dissolution and positive for precipitation, which is consistent with the convention from Section 2.2.

From now on, we focus on the concern of calcite precipitation and crystallization resulting from ${\text{CO}}_2$ injection based on the series of homogeneous reactions (Equation 13) along with the mineral-solute interaction given by Equation 14c. In practice, we enforce a pH greater than 10.33 such that the carbonate ions ${\text{CO}}_3^{2-}$ are the predominant species (see Figure 2). This enables the restriction of the overall set of chemical reactions (Equation 13) to merely consider the Equation 14c in the sense that we assume the intermediate reactions as instantaneous and conservative—without loss of quantity of matter. Such an assumption is acceptable, in practice, since fluid-mineral reaction rates are usually slower than intra-aqueous reaction rates. Therefore, the initial concentration of carbonate ions, denoted $C_{{{\text{CO}}_3}^{2-}}(x,t=0)$ for $(x,t)\in \mathrm{\Omega }\times [0,T_f]$ following the notations introduced in Section 2, is directly related to the partial pressure of injected ${\text{CO}}_2$ by means of the Henry law. The latter states, at a constant temperature, the relation between the amount of dissolved gas in a solute and its partial pressure based on Henry’s law constant denoted $K_H$, which depends on the gas and temperature, such that for the ${\text{CO}}_2$ at $25^{◦}\mathrm{C}$ one gets:

$C_{\text{C}\text{O}2\left(\text{a}\text{q}\right)}=\frac{P_{{\text{CO}}_2}}{K_H}\simeq C_{{\text{CO}}_3^{\mathrm{2-}}}\left(x,t=0\right)(16)$

where $K_H=29.41\mathrm{L}.\mathrm{a}\mathrm{t}\mathrm{m}.\mathrm{m}\mathrm{o}{\mathrm{l}}^{-\mathrm{1}}$. Considering such an alkaline medium—with $\text{pH\hspace{0.17em}}>10.33$ — also results in the treatment of the chemical reaction (Equation 14c) as completely irreversible which corresponds to far-from-equilibrium conditions modelling the calcite precipitation chemical reaction:

$\mathrm{C}{\mathrm{a}}^{\mathrm{2+}}+\mathrm{C}{\mathrm{O}}_{\mathrm{3}}^{\mathrm{2-}}\stackrel{{\mathrm{K}}_{-\mathrm{3}}}{\to }\mathrm{C}\mathrm{a}\mathrm{C}{\mathrm{O}}_{\mathrm{3}\left(\mathrm{p}\right)}(17)$

where the subscript $p$ here refers to the precipitate form of the calcite product. In this case, the rate constants $K_1=K_2=0$ in (35) and the affinity term dependent on the Gibbs energy satisfies the condition $\ln \left(Q/K_{\text{eq}}\right)\gg 0$ — corresponding to a supersaturated solution—so that we obtain an overall reaction rate for the calcite precipitation which reads as:

$R_{\text{prec}}\left(\mathbf{C}\right)=K_{-3}{A}_s{a}_{{\text{Ca}}^{2+}}{a}_{{\text{CO}}_3^{2-}}(18)$

where $K_{-3}=K_3/K_{\text{eq}}$, which theoretically results from the TST law in Equation 15 and has been experimentally validated, inter alia, by Chou et al. (1989).

Therefore, in the following, we rely on the kinetic formulation of the mineral precipitation given by Equation 18, considering the rate laws determined by laboratory experiments (Chou et al., 1989) and normalized by the reactive surface area of the mineral $A_s$ (Plummer et al., 1978). As the geometry evolves, the micro-porosity $\epsilon$ and the reactive specific area $A_s$, associated with the porous structure, also change. These evolutions are taken into account in the reaction rates management and the hydrodynamic modeling of the reactive process (see the overall PDE system (Equation 8) in Section 2.2). The calcite precipitation reaction is subsequently supplemented with a crystallization model which is elaborated in the next Section 4.2.

4.2 Crystal growth modeling: a two-step process

Crystal growth kinetics involves complex mechanisms occurring simultaneously and depending, inter alia, on the concentration of the constituent ions in the solute, but also on attachment frequencies of the ions or precipitates to lattice growth sites (Nielsen and Toft, 1984; Wolthers et al., 2012). Indeed, the growth rate is first controlled by advection and diffusion of the ${\text{Ca}}^{2+}$ and ${\text{CO}}_3^{2-}$ ions to the crystal surface coupled with a surface adsorption process that largely hinges on the crystal lattice shape. For instance, the growth of crystal aggregates is more likely to occur near kinks or corners (Nielsen and Toft, 1984; Yoreo and Vekilov, 2003). Mineral heterogeneity of the pore interface is also an important factor that influences the crystal growth location and morphology, providing preferential sites (Lioliou et al., 2007). This first process is commonly called primary heterogeneous nucleation, for which the crystallization reaction is catalyzed by the solid surface of the porous medium. In the absence of solid interface, crystal clusters can also form spontaneously in the solute, which is known as primary homogeneous nucleation and is closely related to the supersaturation state of the solution in order to initiate the nucleation—namely, satisfying the condition $\ln \left(Q/K_{\text{eq}}\right)\gg 0$. Finally, secondary nucleation occurs in the presence of existing crystals and is more likely to generate large crystal aggregates at the mineral surface. Overall, calcite crystallization results from a combination of all these previous phenomena.

In this section, we consider a two-step crystallization process wherein calcite precipitates, also referred to as nuclei and denoted ${\text{CaCO}}_{3\left(\text{p}\right)}$, are first generated within the solute during the so-called nucleation stage according to the chemical Equation 17. These precipitates are subsequently aggregated at the mineral surface through adsorption phenomena during the crystal growth step. This sequential crystallization process is described in Figure 3, where the notation ${\text{CaCO}}_{3\left(\text{s}\right)}$ stands for the calcite crystal. In the applications developed in Section 5, we consider that the solid matrix of the 3D porous sample has a similar carbonate nature to the calcite crystal generated, though rock mineral heterogeneities can be integrated into the numerical framework as prospects. From now on, we refer to precipitation as the primary homogeneous nucleation and we investigate the surface attachment of the calcite precipitate, ${\text{CaCO}}_{3\left(\text{p}\right)}$, based on an autocatalytic process to model calcite crystal growth, referred to as the secondary surface nucleation in Figure 3. In the reaction scheme from Figure 3, we account for the precipitate diffusion and advection until the solid boundary where it leads to crystal growth through adsorption phenomena. Besides, we neglect the direct crystallization process induced by the solute diffusion to the solid matrix and the so-called primary heterogeneous nucleation.

Figure 3

Figure 3. Overall reactive diagram of the reversible chemical Equation 14c and two-step calcite crystallization process. This diagram represents the chemical interactions between the different ionic species (${\text{Ca}}^{2+}$ and ${\text{CO}}_3^{2-}$), the calcium carbonate precipitate ${\text{CaCO}}_{3\left(\text{p}\right)}$ and mineral crystal ${\text{CaCO}}_{3\left(\text{s}\right)}$ and accounts for both dissolution, precipitation crystal growth processes. In the numerical applications (see Section 5), we mainly address the green part of this diagram that highlights the two-step modeling of the calcite crystallization process.

In the literature, two distinct approaches are mainly developed when considering precipitation and crystal growth modeling altogether. The former can be regarded as “deterministic models” relying on the TST modeling developed in Section 4.1. Noiriel et al., for instance, investigated the effects of pore-scale precipitation on permeability through a combination of X-ray $\mu$CT experiments and “deterministic” modeling (Noiriel et al., 2021). They also derived crystal growth rates directly from the $\mu$CT through an imaging comparison between the beginning and end of the precipitation experiment. In this case, only two $\mu$CT scans were performed, and therefore, the process should be understood as distinct from 4D $\mu$CT experiments incorporating time dynamics. Such experimental identification of crystal growth rates can, however, be prone to intrinsic imaging limitations (Perez et al., 2022) and result in wide discrepancies in the reaction rate estimations. Nonetheless, their results showed satisfactory agreement between the experiments and numerical experiments for precipitation processes into fractures (Noiriel et al., 2021). Alternative modeling approaches lie in the probabilistic nature of nucleation or crystal growth and are referred to as “probabilistic models” (Fazeli et al., 2020; Masoudi et al., 2021; Nooraiepour et al., 2021a). Wolthers et al., for instance, developed a probabilistic approach for calcite crystal growth based on the nature of the kink sites depending on their ionic affinities and attachment frequencies of the constituent ions (Wolthers et al., 2012). Estimations of such adsorption frequency ranges can also be found in the literature (Christoffersen and Christoffersen, 1990; Nielsen, 1984). Finally, while it is commonly established experimentally that crystal growth occurs preferentially at kinks and corners of the surface lattice, few models incorporate the geometrical dependency of the crystal aggregation in the reaction rates (von Wolff et al., 2021).

In the present article, we propose a new approach coupling a deterministic model for the precipitation, which directly depends on the supersaturation ratio following the TST formalism, and a probabilistic formulation of the crystal growth process. The latter accounts for the adsorption frequencies of the precipitate to the growth sites with a coefficient, quantifying the physical probability of attachment rate, denoted $P_{\text{ad}}$, which relies on a locally averaged mineral volume fraction. In this formalism, one obtains the crystal growth reaction rate, which is expressed as:

$R_{\text{crys}}\left(\mathbf{C}\right)=K_c{P}_{\text{ad}}{C}_{{\text{CaCO}}_{3\left(\text{p}\right)}}(19)$

with $K_c$ the adsorption frequency (expressed in ${\mathrm{s}}^{-\mathrm{1}}$). It possible then to include the geometrical dependency in the crystal growth reaction rate by means of involving the neighborhood of the solid using a convolution of the porosity (the size of the neighborhood is provided by the support of the convolution kernel). Practically, this leads to the relation

$P_{\text{ad}}=C\left(1-\epsilon *W_h\right){\epsilon }^m,(20)$

where $W_h(x)=h^{-d}{W}^{\otimes d}(x/h)$ is the rescaled kernel based on a local averaging kernel $W$, so that $(1-\epsilon *W_h)\epsilon$ is enhanced in the layer close to the fluid/solid interface and depending on the solid proportion in the neighborhood (see Figures 4, 5). The calibration of $C$ and kernel support size $h$ can be established by comparison with experiments (Poonoosamy et al., 2023) or numerical simulations in standardized geometry (Varzina et al., 2020; Patel et al., 2021), or estimated in the same spirit as the crystallites used in the LBM schemes in Masoudi et al. (2024).

Figure 4

Figure 4. Probability of attachment rate $P_{\text{ad}}=C(1-\epsilon *W_h){\epsilon }^m$ for crystal growth modeling in 1D different contexts. Synthetic continuous and discrete (at the voxel size) representation of a 1D porous medium, with a residual micro-porosity in the porous matrix ${\epsilon }_0=5\%$ for $x>0$, for sharp and smooth interface. The probability values are computed based on Equation 20. (A) Sharp interface, continuous. (B) Sharp interface, discrete. (C) Smooth interface, continuous. (D) Smooth interface, discrete.

Figure 5

Figure 5. Impact of the probability of attachment rate $P_{\text{ad}}$ for crystal growth modeling. Synthetic representation of a discrete (at the voxel size) porous medium with smooth interface (t be compared with Figure 4), with a rock matrix micro-porosity ${\epsilon }_0=5\%$. The probability values are computed based on Equation 20 and include the one-neighborhood (e.g., the red square).

In the case of $K_c$ being a production of crystal volume per unit of volume and time, then the normalization coefficient $C$ satisfies $1/C={⟨(1-\epsilon *W_h){\epsilon }^m⟩}_{\mathrm{\Omega }}$ so that $P_{ad}$ is a probability distribution, and the index $m$ describes the ability to decrease strongly the reaction rate inside the solid in the spirit of the Archie’s law (in practice $m=2$ is suitable). In the case of $K_c$ being the local production of crystal at the fluid/solid interface, which is the case considered in the present study, then $C$ satisfies

$1/C=\left(1-{\epsilon }_0\right)/2$

so that $P_{ad}$ is a point-wise probability of capture by the crystal. This is based on the jump between ${\epsilon }_0$ and 1 through a sharp interface, with a factor $1/2$ (half the integral of the kernel $W$, which can be slightly adjusted for concave interface as shown in Poncet (2007)). This convolution-based formulation is appropriate for a crystallization process and is inspired by the gradient-based technique from Luo et al. (2012) and Soulaine et al. (2017) that locates the first layer on the solid side and is suitable for dissolution processes. In practice, the local averaging kernel $W$ can be as simple as $W={\mathbb{1}}_{[-\mathrm{1,1}]}/2$ or Gaussian in a continuous description, or approximated by its discrete formulation $W=({\delta }_{-1}+{\delta }_0+{\delta }_1)/3$, or even with weighting in order to modulate the length of capture with respect to the grid size.

The Figure 4 shows different probability distributions in 1D with a residual micro-porosity of ${\epsilon }_0=5\%$ (in the area $x>0$), for sharp and smooth fluid/solid interface, and for the continuous and discrete formulations described above. The resulting probability values in 2D are represented in Figure 5 for a synthetic example with smooth interface and discrete formulation (to be compared to the 1D version of Figure 4D), where the residual micro-porosity in ${\mathrm{\Omega }}_S$ is estimated to ${\epsilon }_0=5\%$ as well.

4.3 Final system of PDEs

Finally, we define the vector of concentrations

$\mathbf{C}=\left(C_{{\text{CaCO}}_{3\left(\text{s}\right)}},C_{{\text{CaCO}}_{3\left(\text{p}\right)}},C_{{{\text{CO}}_3}^{2-}},C_{{\text{Ca}}^{2+}}\right)$

and consider the reactions rates $R_{\text{prec}}(\mathbf{C})$ and $R_{\text{crys}}(\mathbf{C})$ respectively given by Formula 40, 41. Overall, the calcite crystallization modeled as a two-step process of precipitation and crystal growth, according to the reaction scheme from Figure 3, writes:

$\left\{\begin{array}{c}-\mathrm{d}\mathrm{i}\mathrm{v}\left(2\mu D\left(u\right)\right)+\mu {\kappa }_b^{-1}\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^2}u=\epsilon \left(f-\nabla p\right),\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ \frac{\partial C_{{{\text{CO}}_3}^{2-}}}{\partial t}+\mathrm{d}\mathrm{i}\mathrm{v}\left(F\left(C_{{{\text{CO}}_3}^{2-}}\right)\right)=-R_{\text{prec}}\left(\mathbf{C}\right),\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ \frac{\partial C_{{\text{Ca}}^{2+}}}{\partial t}+\mathrm{d}\mathrm{i}\mathrm{v}\left(F\left(C_{{\text{Ca}}^{2+}}\right)\right)=-R_{\text{prec}}\left(\mathbf{C}\right),\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ \frac{\partial C_{{\text{CaCO}}_{3\left(\text{p}\right)}}}{\partial t}+\mathrm{d}\mathrm{i}\mathrm{v}\left(F\left(C_{{\text{CaCO}}_{3\left(\text{p}\right)}}\right)\right)=R_{\text{prec}}\left(\mathbf{C}\right)-R_{\text{crys}}\left(\mathbf{C}\right),\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ \frac{\partial C_{{\text{CaCO}}_{3\left(\text{s}\right)}}}{\partial t}=R_{\text{crys}}\left(\mathbf{C}\right),\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ \epsilon =1-v{C}_{{\text{CaCO}}_{3\left(\text{s}\right)}},\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ F\left(C_{•}\right)={\epsilon }^{-1}u{C}_{•}-{\alpha }_{•}\left(\epsilon \right){\nabla }^{\epsilon }{C}_{•},{\alpha }_{•}\left(\epsilon \right)={\epsilon }^{\eta }{D}_{m,•}\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\times ]0,T_f[\hfill \\ +\text{adequate boundary and initial conditions, along with}\mathrm{d}\mathrm{i}\mathrm{v}u=0\hfill \end{array}\right.(21)$

where the advective and diffusive flux $F(C_{•})$ is taken from expressions in Section 2.2, using the lifted gradient ${\nabla }^{\epsilon }(\cdot )=\epsilon \nabla ({\epsilon }^{-1}\cdot )$. The reactive hydrodynamic model (44) ensures that part of the precipitate, generated in the solute through homogeneous nucleation, is transferred to the mineral surface by adsorption. One should notice that in this model the precipitate $C_{{\text{CaCO}}_{3\left(\text{p}\right)}}$ is both advected and diffused. Such diffusion enables to account for the interaction of the precipitates with potential unresolved roughness or features in the porous matrix ${\mathrm{\Omega }}_S$.

In the next Section 5, we apply this two-step crystallization model to DNS of ${\text{CO}}_2$ mineral trapping into a real rock geometry at the pore-scale. We investigate the morphological changes in the porous matrix structure, the clogging of pore throats, and the evolution of the macro-scale properties, namely, the porosity and permeability.

4.4 High-performance-computing aspects

One of the major constraints when dealing with a semi-Lagrangian formulation lies in the ability of the computational framework to handle an overall hybrid approach in terms of grid-particle formalism, numerical methods, multi-grid resolutions, and hardware devices. Indeed, Cottet et al. (2009) suggested a semi-Lagrangian method coupled with hybrid grid resolutions to address a multi-scale transport problem of a passive scalar. The scalar is discretized on a sub-grid compared to the velocity and vorticity fields and enables the accurate prediction of the small-scale effects. Finally, considering hybrid computing methodologies makes it possible to distribute the distinct parts of an overall problem to different types of hardware architectures. This formalism, therefore, exploits the advantages of each method individually according to the characteristics of the problem.

Nonetheless, implementing such a hybrid approach requires a highly flexible computational framework that gathers a wide range of numerical methods and, therefore, benefits from their intrinsic strengths. One also needs libraries incorporating effective parallel computing tools and able to address, inter alia, hybrid CPU-GPU programming. In this section, we present the HPC framework considered to address this DNS of pore-scale reactive flows for ${\text{CO}}_2$ mineral trapping into carbonate rocks.

The library HySoP is a high-performance computing platform (Etancelin et al., 2022), jointly developed at LMAP (Laboratoire de Mathématiques et de leurs Applications, UMR 5142 CNRS, UPPA), LJK (Laboratoire Jean Kuntzmann, Alpes-Grenoble University, UMR 5224 CNRS), and M2N (Laboratoire Modélisation mathématique et numérique, Conservatoire National des Arts et Métiers–CNAM, Paris, EA 7340 CNRS). The library, originally developed to address flow simulations based on remeshed particle methods on hybrid muti-CPU and multi-GPU architectures, was initiated by the work of Etancelin (2014) and has been successfully extended to a larger scope of HPC applications including dissolution at pore-scale (Etancelin et al., 2020).

The entire code is structured around the operator splitting strategy that defines the different operators involved in a problem (at a high level of abstraction), and afterward, enables the discretization of these operators, which are solved using the most appropriate numerical method (at a lower level of abstraction). The overall problem is formally described through an acyclic graph that expresses the operators interactions in the splitting formulation through data dependencies (Etancelin, 2014; Etancelin et al., 2020; Keck, 2019). Even if the code is mainly written in Python, the numerical methods are either implemented using compiled language that enables threads parallelism (using OpenMP on CPU and OpenCL on accelerators) or taken from external libraries. User interface enables building together both numerical methods provided with HySoP on cartesian grid and user defined Python functions. In the latter, performances are obtained provided the usage of interfaces to compiled language (i.e., among others f2py, swig, cython, numba, …). A complementary distributed parallelism is naturally available using domain decomposition implemented with a Message Passing Interface library (MPI).

One of the core features of HySoP is to target hybrid computing using both CPU and GPU. The latter was made possible by the emergence, in the 2000s, of the so-called GPGPU concept, which integrates the GPU as a CPU co-processing partner targeting accelerated performances. The choice of OpenCL was made to cope with portability of the performances as a generic multicore artichecture programming standard (Stone et al., 2010). The initial GPU computing feature was latter enhanced with code generation from templates or symbolic mathematical expressions together with auto-tuning tanks to micro-benchmarks at run time (Keck, 2019). Just-in-time compiling is extensively used to achieve performance portability and lazy specification of kernel parameters which known as challenging problem widely dependent on the hardware architecture (Dolbeau et al., 2013).

The previous numerical method have been implemented on a hybrid computing strategy. Darcy-Brinkmann-Stokes equation is solved on CPU architecture using full MPI parallelization technique while all the reactive and transport part is computed using OpenCL on GPUs.

Overall profiling is represented on Figure 6 for two cases of the next section. Results are showing that recent hardware is capable of better performances without any changes in the code. We benefit here from OpenCL code generation and micro-benchmarks on accelerator and usual compiler optimizations for CPU part. Significant variation appears in computational time associated to flow steps because larger ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}$ implies more intense geometrical changes and then reorganization of main flowpath, as described later in Section 5. One can notice the quite large proportion of computational time spend in data management specially for reactive transport part. It corresponds to ghosts layers that handle boundary conditions and whole data transposition in order to get the first varying index identical to the current direction of directional splitted advection and remeshing operator. Ghosts layer widths are quite large due to large CFL numbers handled here. Data management on GPU has been identified as a known bottleneck of the code and would be investigated in future high performance computing engineering efforts.

Figure 6

Figure 6. Performances in various hardware configuration and discretization: The considered hardware is NVIDIA P100 and H100 GPUs and 16-cores Intel Xeon 6130 cpus and 24-cores AMD EPYC 9254. The computational time is splitted according to the algorithm steps, including operator splitting and data management. Both cases where run with $Pe=4.47$ (see Section 5 for details).

5 Results and discussion on clogging/non-clogging regimes

5.1 Pore-scale and reactive setup

The present section focuses on the effects of calcite crystallization on changes in the pore geometry, macro-properties, and flow at the pore scale in the context of ${\text{CO}}_2$ mineral trapping. We consider a pore-scale geometry obtained by microtomography measurements from Sheppard and Prodanovic (2015) and freely available on the Digital Rocks data portal, which includes $\mu$CT datasets of limestone, glass bead pack, and Castlegate sandstone. The numerical simulations are performed on the Castlegate geometry at a resolution of $128^3$ with a voxel size of $5.6\mu$m, which represents a numerical sub-sample of about $L=0.7168$mm. We assume, as previously introduced in Section 4.2, that the porous matrix is of identical mineral nature as the crystal generated along the reactive process and, thereby, consider that the sub-sample is composed of calcite. While including mineral heterogeneities as perspectives, we here hypothesize the homogeneity of the mineral structure within the sample. Finally, the specific area is numerically estimated for this sample to get, at the initial state, $A_s=8300{\mathrm{m}}^{-\mathrm{1}}$, and is afterward updated along the reactive process.

Numerical simulations are performed under atmospheric conditions in terms of pressure and temperature and rely on the experimental identification of the reaction rate constants, arising from the literature (Chou et al., 1989). We consider isothermal conditions with a temperature of $T=25^{◦}\mathrm{C}$ and an injection of ${\text{CO}}_2$ with a partial pressure of $P_{{\text{CO}}_2}=3.15\times 10^{-2}\mathrm{b}\mathrm{a}\mathrm{r}=2.96\times 10^{-2}\mathrm{a}\mathrm{t}\mathrm{m}$ — which is about 100 times greater than the partial pressure of ${\text{CO}}_2$ in the atmosphere. Given Henry’s law constant for the ${\text{CO}}_2$ at $25^{◦}\mathrm{C}$, and under the assumption of a highly alkaline medium—with $\text{pH\hspace{0.17em}}>10.33$ — we estimate from Equation (16) that the initial concentration of carbonate ions is given by $C_{{{\text{CO}}_3}^{2-}}(x,t=0)=10^{-3}\mathrm{m}\mathrm{o}\mathrm{l}.{\mathrm{L}}^{-\mathrm{1}}$.

The calcium initial concentration is subsequently determined based on the equilibrium constant $K_{\text{eq}}=10^{-8.48}$ (Chou et al., 1989; Plummer and Busenberg, 1982) to ensure a far-from-equilibrium precipitation regime given by the supersaturation condition $Q\gg K_{\text{eq}}$. We assume that the medium pore space is initially filled with a saturated solution wherein the initial concentration of calcium ions is $C_{{\text{Ca}}^{2+}}(x,t=0)=0.1\mathrm{m}\mathrm{o}\mathrm{l}.{\mathrm{L}}^{-\mathrm{1}}$. Therefore, in our case, the saturation in calcium ions ${\text{Ca}}^{2+}$ initially present in the domain is not a limiting factor of the precipitation reaction, and we consider a continuous calcium injection that maintains the supersaturation constraint. Actually, in order to maintain this supersaturation, we will assume that the concentration in ${\text{Ca}}^{2+}$ remains constant at its initial value. Comparable initial conditions and assumptions have been employed in investigating dissolution experiments by Maes et al. (2022), wherein the sample core was initially flooded with a brine solution that had previously reached equilibrium with supercritical ${\text{CO}}_2$.

Regarding the reaction rate constant for the precipitation, we rely on the experimental identification from Chou et al. (1989) such that:

$K_{-3}=\frac{K_3}{K_{\text{eq}}}=\frac{\mathrm{6.6·}{10}^{-7}}{10^{-8.48}}=199\mathrm{m}\mathrm{o}\mathrm{l}.{\mathrm{m}}^{-\mathrm{2}}.{\mathrm{s}}^{-\mathrm{1}}$

while adsorption frequencies $K_c$ commonly encountered in the literature range from $10^3$ to $10^8{\mathrm{s}}^{-\mathrm{1}}$ (Christoffersen and Christoffersen, 1990; Nielsen, 1984; von Wolff et al., 2021; Wolthers et al., 2012). In practice, we set for the numerical simulations $K_c=10^3{\mathrm{s}}^{-\mathrm{1}}$, the molecular diffusion $D_m=10^{-9}{\mathrm{m}}^{\mathrm{2}}.{\mathrm{s}}^{-\mathrm{1}}$ for all the species and the prescribed flow rate $\bar{u}=1.{10}^{-3}\mathrm{m}.{\mathrm{s}}^{-\mathrm{1}}$.

5.2 Pertinent non-dimensional numbers

Investigating the different precipitation patterns and regimes relies on the definition of well-established characteristic dimensionless numbers, namely, the Peclet and (second or catalytic) Damköhler numbers respectively denoted $\mathrm{P}\mathrm{e}$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}$ (Noiriel et al., 2021; Soulaine et al., 2017; Steefel and Lasaga, 1990). However, in the context of the joint precipitation and crystal growth modeling, we define two distinct Damköhler numbers characterizing each process and respectively denoted ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}$. These dimensionless numbers are subsequently defined as:

$\mathrm{P}\mathrm{e}=\frac{\bar{u}L}{D_m},{\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}=\frac{K_{-3}{\gamma }_{{\text{CO}}_3^{\mathrm{2-}}}{A}_s{L}^2}{D_m}\text{\hspace{0.17em}and\hspace{0.17em}}{\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=\frac{K_c{L}^2}{D_m}$

where $\bar{u}$ and $L$ are respectively the characteristic velocity and length of the sample, and the activity coefficient of the carbonate ions is $\gamma {\text{CO}}_3^{2-}=10^{-3}{\mathrm{m}}^{\mathrm{3}}.\mathrm{m}\mathrm{o}{\mathrm{l}}^{-\mathrm{1}}$. The characteristic length $L$ can be related to pore size (Steefel and Lasaga, 1990), though it is commonly set as $L=\sqrt{{\kappa }_0}$ provided an experimental or numerical estimation of ${\kappa }_0$ (Soulaine et al., 2017). The latter alternative is applied here, with an estimation of ${\kappa }_0=$ 2$\times$10⁻¹¹ ${\mathrm{m}}^{\mathrm{2}}$ for the porous sample considered.

5.3 Evidence of two clogging regimes at same Da_II^prec and two different Da_II^crys

The first crystallization regime we investigate is characterized by the following dimensionless numbers $\mathrm{P}\mathrm{e}=4.47$, ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}=33.034$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=20$. Precipitation and crystallization of calcite lead to a significant decrease in the macro-scale permeability and porosity, resulting from flow path disruptions at the micro-scale through partial or complete clogging of pore throats. This can also affect the roughness of the mineral interface and the pore-size distribution of the sample and, thereby, contribute to influencing the sample hydrodynamics properties. In particular, we observe these effects at the pore scale in Figure 7 along the reactive process and for several physical times $t$ going from $2\mathrm{h}45$ to $T_f=6\mathrm{h}56$. On the right side of Figure 7, we depict partial views of the porous sample’s morphology, illustrating the changes in pore structure over the reaction time, along with the micro-porosity field $\epsilon$ within the porous matrix ${\mathrm{\Omega }}_S$. On the left side, we represent along a slice in the main flow direction (taken at $z=-0.0168$mm), the local variations on the micro-porosity with respect to the initial state—before the reaction process—given by $\epsilon (t)-\epsilon (0)$. Initially, we notice that higher micro-porosity variations are more likely localized at the mineral interfaces but also near thin pore throats.

Figure 7

Figure 7. Time evolution of the sample geometry and micro-porosity at the pore scale, illustrating pore-clogging effects. Slice at $z=-16.8\mu$m of the porosity variations $\epsilon (t)-\epsilon (0)$ in the fluid region of the pore space for various times $t$, on the left. Partial view of the pore space structure as an isosurface of $\epsilon (t)$ for several times $t$, on the right. The clogging occurs in the red areas, mainly upstream (flow direction is given by the arrow on picture (A) Time $t$ = 2h45. (B) Time $t$ = 5h30. (C) Time $t$ = $T_f$ = 6h56.

These variations subsequently lead to pore-clogging and reorganization of the main flow pathway (see Figure 7C). In Figure 8B, we investigate the effects of such micro-scale changes on the evolution of the petrophysical properties at the macro-scale, namely, the porosity $\varphi$ and permeability ${\kappa }_0$ (upscaled quantities from Figure 1). The results are consistent with the expected decrease along the ${\text{CO}}_2$ mineral trapping process but also highlight sharp permeability drops, which characterize the pore-clogging phenomena. Finally, in order to identify more clearly the crystallization pattern in this particular regime, we display in Figure 8A an isosurface of the micro-porosity variation between the final and initial times. This illustrates that micro-scale variations occur preferentially in a compact and non-uniform manner in the first inlet part of the domain while following the individual ramifications in the pore structure.

Figure 8

Figure 8. Evolution of the sample properties along the reactive process: (A) Micro-porosity evolution represented by an isosurface of the porosity variation $\epsilon (T_f)-\epsilon (0)$ at half of its maximum value. Results after nearly 7 h of precipitation and crystal growth, illustrating a non-uniform compact regime following the natural ramification of the sample. (B) Evolution of the macro-scale properties, porosity $\varphi$ and permeability ${\kappa }_0$, along the two-step crystallization process from Figure 3.

We subsequently investigate the impact of different dominant regimes on the overall two-step crystallization process. To do so, we consider both transport dominant cases with $\mathrm{P}\mathrm{e}=4.47>1$ and diffusion dominant cases with $\mathrm{P}\mathrm{e}=0.447<1$, coupled with two different crystal growth regimes characterized by ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=2$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=20$. One should notice that the effect of precipitation Damköhler ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}$ changes are not analyzed since this number, characterizing the first homogeneous nucleation regime, is a limiting factor of the crystallization process from Figure 3. Here, we assume in all the previous cases that ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}<{\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}$ which guarantees a supersaturation state suitable to the development of crystal aggregates on the mineral surface.

To the best of our knowledge, considering that the crystallization regime can be driven by three distinct dimensionless numbers is one of the contributions of the present article. Indeed, most of the regime diagrams presented in the literature mainly characterize precipitation patterns according to the $\mathrm{P}\mathrm{e}$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}$ dimensionless numbers, which implies neglecting the effects of nuclei adsorption at the mineral surface in the different regime configurations (Tartakovsky et al., 2007; Yang et al., 2021). However, our results indicate that both homogeneous calcite nucleation (the precipitation step in Figure 3) and growth stages are important in the development of precipitation and crystallization patterns. In particular, we establish that the crystal growth Damköhler number ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}$ has a non-negligible impact on precipitation pattern and porosity variations along the reactive process, regardless of the other dimensionless numbers $\mathrm{P}\mathrm{e}$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}$ commonly investigated.

Indeed, Figures 9B, 10B highlight two distinct crystallization regimes at similar $\mathrm{P}\mathrm{e}$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}$, but with a ten times smaller adsorption frequency $K_c$ in Figure 9B which impact the ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}$. Figure 9B shows that, for a small adsorption frequency $K_c$, the calcite precipitate is uniformly generated and advected along the main flow path direction (due to the transport dominant regime with $\mathrm{P}\mathrm{e}>1$) while the micro-porosity changes are minimal. This illustrates that the main flow path is maintained since few calcite nuclei aggregate to the mineral surface. By increasing the adsorption frequency in Figure 10B, the precipitate formation concentrates on the inlet of the domain and is less subject to advection, while the micro-porosity changes become significant. On the one hand, this highlights pore-clogging that prevents further transport of the calcite nuclei. On the other hand, we also notice a backward increase in calcite precipitates that accumulate behind the clogging after some time. The same analysis holds for Figures 9A, 10A, except we consider a diffusion dominant regime with $\mathrm{P}\mathrm{e}<1$, which means that the precipitate transport is reduced so that the nuclei formation and micro-porosity changes are even more constrained to the inlet boundary of the domain. The results presented in Figures 7 and 8 correspond to the crystallization regime identified by $\mathrm{P}\mathrm{e}>1$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=20$ in Figure 10B. This confirms that the permeability drops observed in Figure 8B are characteristics of pore-clogging effects.

Figure 9

Figure 9. Crystallization regimes at ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=2$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}=33$ for Peclet numbers below 1 (A) and above 1 (B): no clogging. Slice integrals of the precipitates and macro-porosity—computed over 2D YZ directional slices—plotted with respect to the main flow path direction coordinates x (in mm mm) and where each curve (using segmented colors) represents a distinct time in the reactive process. The porosity integral remains close to its initial value, which shows that pores remain open with the same flow, hence no clogging effect.

Figure 10

Figure 10. Crystallization regimes at ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=20$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}=33$ for Peclet numbers below 1 (A) and above 1 (B): pore-clogging. Slice integrals of the precipitates and macro-porosity—computed over 2D YZ directional slices—plotted with respect to the main flow path direction coordinates x (in mm mm) and where each curve represents a distinct time in the reactive process. The porosity integral goes from 30% to 10% (close to the cristal porosity at 5%) in the upstream part of the sample, which shows a pore-clogging effect.

To conclude, at ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=2$, ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}=33$ and for $\mathrm{P}\mathrm{e}=0.447<1$ as well as for $\mathrm{P}\mathrm{e}=4.47>1$, Figure 9 shows that there is no clogging measured. At ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=20$, ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}=33$ and for the two same Peclet numbers, Figure 10 displays evidence of pore-clogging. The clogging is consequently driven principally by the non-dimensional number ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}$ in this case study, the Peclet number affecting the depth of the clogging effect. This is valid as long as the ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}$ is sufficiently high in order to provide precipitate able to turn into crystal.

5.4 Discussion

This section aims at positioning the final set of Equation 21 with the existing literature on nucleation. These models involve the physics of fluid precipitation (homogeneous nucleation) and its subsequent crystallization (heterogeneous nucleation), while as shown on Figure 3, the direct crystallization from reactants is also another possible kind of nucleation. The later direct crystallization is not considered in the present study, but can be compared to our two-step modeling as long as ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}\gg {\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}$ since we don’t model the precipitate rheology.

Moreover, the existing literature considers either internal flows (porous medium as a pore collection) or external flows (porous medium as a grain collection, also referred to as pore-scale). Every study has its own configuration and Dahmköhler/Peclet setup that we cannot investigate one by one in this study.

More in detail concerning recent literature, confined internal flows have been studied experimentally in Poonoosamy et al. (2023); Noiriel et al. (2021) and numerically in Molins and Knabner (2019), but most geometries of porous media described at their pore-scale for systematic analysis are grain-shaped and involve the hydrodynamics of low Reynolds external flow. Such external flow analysis are numerical (Starchenko, 2022; Nooraiepour et al., 2021b; Yang et al., 2021; Varzina et al., 2020; Fazeli et al., 2020) or experimental (Nooraiepour et al., 2021a), some of them focusing on the wall as an initially flat fluid/solid interface Deng et al. (2022). Also, some configuration are in between confined and grain-shaped geometries (Masoudi et al., 2024), with some focus on cement material (Tong et al., 2024; Patel et al., 2021).

Among these studies, qualitative comparison with Yang et al. (2021), and quantitative comparison with Masoudi et al. (2024) are of direct interest for the present investigation.

First, in Yang et al. (2021), a diagram showing the expected nucleation regime with respect to ${\mathrm{D}\mathrm{a}}_{\text{II}}$ and $Pe$ is available on their Figure 13. They consider crystal growth around a grain, studying the hydrodynamic and reaction balance and the resulting grain shape. Their model rely on Kozeny-Carman correlation law and TST so only the numerical method significantly differs (they use LBM while we use semi-Lagrangian with PSE method). Their results are compatible to ours, as for $1⩽{\mathrm{D}\mathrm{a}}_{\text{II}}⩽100$ and $\mathrm{P}\mathrm{e}⩽10$ we don’t observe dendrite formation nor uniform distribution of crystallization.

Second, the article Masoudi et al. (2024) considers only heterogeneous nucleation with both spherical grains and pore structures at several crystallites levels and links porosity and permeability during crystallization process. They consider Lattice Boltzmann simulations of crystal growth from aqueous material with a rate $R_G=k_{G}S$ in mol/s which is no more than our relation (41). Despite the formulation of crystallites events fitting well to LBM expression of reaction formula but far from our deterministic PDE formulation, we are comparable to their heterogeneous geometry with high crystallite level. Our configuration at ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=20$ contains region of high crystallization rate and some not reached by reactants so the comparison is pertinent to the cases B and D of the Figure 3 from Masoudi et al. (2024). On Figure 11 we plot the similar petro-physical diagram showing the renormalized permeability $K/K_0$ with respect to renormalized global porosity $\varphi /{\varphi }_0$. In order to get a pertinent measurement of the crystallization process, we restrict the computation of $K$ and $\varphi$ to the upstream quarter of the computational domain, so that the porosity displayed on Figure 10 for $\mathrm{P}\mathrm{e}>1$ evolves from 30% to 13.9% in 3h40, that is to say $\varphi /{\varphi }_0=$46.5% of the initial porosity remains.

Figure 11

Figure 11. Renormalized $K-\varphi$ comparison with Masoudi et al. (2024) indicators for the regimes displayed on Figures 9, 10. Kozeny-Carman relationship developed in Equation 2. (A) K-$\varphi$ relation, long time scale. (B) K-$\varphi$ relation, zoom on small effects

With this setup, the renormalized $K-\varphi$ diagram is displayed on the Figure 11A, along with the power laws $K/K_0={(\varphi /{\varphi }_0)}^n$ for $n=\mathrm{2,3,8,64}$, the case $n=3$ being related to the Kozeny-Carman correlation far from $\varphi =1$. We exhibits on this diagram our 4 configurations (two different Dahmköhler and Peclet number) for the whole simulation and for the strongest crystallization events, corresponding to individual pore clogging time window and identified by a sharp decrease of the permeability (see Figure 8 for instance). At a global level we observe clogging events with exponents $n$ between 64 and 8 as expected [see Figure 3D in Masoudi et al. (2024)], but a lower exponent around 1.5 at long time scale [2–4 was expected from Figure 3B in Masoudi et al. (2024)]. Concerning the time evolution of our regime, the Figure 11B shows a zoom at high $K$ and $\varphi$ values, exhibiting short time scale (the maximum renormalized values 1 corresponding to the initial condition). We can then read that for short time scale at ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=20$, the Kozeny-Carman law is quite well followed but deviates strongly after the last individual pore clogging $(t=1h25)$. Moreover, the low Dahmköhler case ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{crys}}=2$ that exhibits no clogging (see Figure 9 for instance) has very limited effect on permeability and porosity (see Figure 11B), but shows nevertheless a small crystallization event at low Pecler number that is also coherent in term of exponent.

Moreover, to conclude with this comparison, such a $K-\varphi$ criterion can be applied only when the full sample is involved in the crystallization process, otherwise the porosity will reach a limit and the permeability will still decrease, which is not compatible with a power-law. Consequently, a special care has to be taken for the computational domain of this criterion when considering a preferential spot of cristallization (such as our case of upstream furring).

Finally, in Starchenko (2022), nucleations occur around a single grain but can be compared to our model (Equation 21) if reduced to 1D (our ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}$ values are similar, 34 and 33 respectively). Deng et al. (2022) shows growth on a wall, also comparable if reduced to 1D. These two growth rates at wall scale will be the topic of future investigation, but are currently too far, at a geometrical level, from our pore-scale configuration.

6 Concluding remarks

This article focused on developing an efficient DNS framework to address reactive flows at the pore scale in the context of ${\text{CO}}_2$ mineral storage. Indeed, the injected ${\text{CO}}_2$ will interact with the aquifer structure and eventually lead to mineral trapping in the form of calcite precipitates and crystals. These processes are interesting to study at the pore scale to ensure a comprehensive analysis of the local rock-fluid interactions and evolving pore structures. This can subsequently translate into meaningful estimations of the macro-scale properties changes and measure the impact of the geochemical processes on the natural underground reservoirs. In particular, precipitation and crystallization lead to a significant reduction in the macro-scale permeability and porosity, which result from partial or complete pore clogging and thus from a reorganization of the flow path at the micro-scale.

From a conceptual perspective, we developed a new crystallization model that efficiently combines a classical deterministic TST approach of the nucleation process with a probabilistic view of the crystal aggregation to the pore surface. This enables us to account for spatial and geometrical dependency in the crystal growth modeling through a probabilistic attachment rate depending on local mineral volume fraction. In this sense, we integrate the modeling of preferential growing sites that largely hinge on the surrounding pore arrangement. To the best of our knowledge, such considerations are here accounted for the first time to model crystallization processes in complex 3D geometries at the pore scale. Investigating probabilistic attachment rates based on the surrounding pore structure also ensures reliable prediction of pore-clogging at the pore scale. Finally, we demonstrate that the proper characterization of crystallization regimes both depends on the nucleation process and crystal aggregation. Indeed, we exhibit that the two commonly considered dimensionless numbers, $\mathrm{P}\mathrm{e}$ and ${\mathrm{D}\mathrm{a}}_{\text{II}}^{\text{prec}}$, are not sufficient to explain clogging effects and precipitation patterns. A novelty of the present manuscript is, therefore, that the crystallization regimes are characterized by three dimensionless numbers that include the effects of nuclei adsorption to the pore surface.

This reactive hydrodynamic model consistently couples a Lagrangian formulation for the reaction equations with a grid-based approach for the flow using the DBS equation with the superficial velocity formalism. This semi-Lagrangian method is addressed through a splitting operator strategy coupled with high-order remeshing steps for grid-particle interpolations. This original numerical method introduced to solve this coupled model has been efficiently incorporated into the hybrid numerical framework HySoP and results in a CPU-GPU implementation of the method. It includes a robust estimation of the heterogeneous diffusion operator arising from Archie’s law term in the reactive system.

At the same time, this article also demonstrated strong implications in the overall reactive system of several parameters that can be subject to a wide range of discrepancies. In particular, morphological features and kinetic parameters, such as the micro-porosity $\epsilon$, specific area $A_s$, rate constants $K_i$, and adsorption frequencies $K_c$, have a significant impact on the reaction rates and the dynamical patterns.

Experimental determination of these parameters can range over several orders of magnitude and result in highly different regimes that drastically affect the estimation of the macro-scale properties: the adsorption frequencies $K_c$ commonly found in the literature can range from $10^3$ to $10^8{\mathrm{s}}^{-\mathrm{1}}$ as shown in Christoffersen and Christoffersen (1990); Nielsen (1984); von Wolff et al. (2021); Wolthers et al. (2012). Our present method has been shown to compute accurately the phenomena of crystallization and precipitation in different regimes, and will be used intensively in future works for inverse problems in order to get a robust and accurate estimation of such adsorption frequencies from experiments imaging. New methods relying of AI and BPINNs for the estimation of reaction rates at the pore scale are already available for dissolution (Perez and Poncet, 2024) and could be adapted to nucleation process in the near future.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

SP: Conceptualization, Formal Analysis, Investigation, Methodology, Resources, Software, Validation, Writing–original draft, Writing–review and editing. J-ME: Investigation, Resources, Software, Validation, Visualization, Writing–original draft, Writing–review and editing. PP: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Visualization, Writing–original draft, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was partially supported by ANR Grant ANR-20-CE45-0022, Carnot Institute ISIFoR Grant P450902ISI and E2S-UPPA project MicroMineral.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2025.1493305/full#supplementary-material

References

Aagaard, P., and Helgeson, H. C. (1982). Thermodynamic and kinetic constraints on reaction rates among minerals and aqueous solutions; I, Theoretical considerations. Am. J. Sci. 282, 237–285. doi:10.2475/ajs.282.3.237

CrossRef Full Text | Google Scholar

Allaire, G., Brizzi, R., Mikelić, A., and Piatnitski, A. (2010). Two-scale expansion with drift approach to the Taylor dispersion for reactive transport through porous media. Chem. Eng. Sci. 65, 2292–2300. doi:10.1016/j.ces.2009.09.010

CrossRef Full Text | Google Scholar

Allaire, G., and Hutridurga, H. (2012). Homogenization of reactive flows in porous media and competition between bulk and surface diffusion. IMA J. Appl. Math. 77, 788–815. doi:10.1093/imamat/hxs049

CrossRef Full Text | Google Scholar

Andrew, M., Bijeljic, B., and Blunt, M. J. (2013). Pore-scale imaging of geological carbon dioxide storage under in situ conditions. Geophys. Res. Lett. 40, 3915–3918. doi:10.1002/grl.50771

CrossRef Full Text | Google Scholar

Angot, P. (2018). Well-posed Stokes/Brinkman and Stokes/Darcy coupling revisited with new jump interface conditions. ESAIM M2AN 52, 1875–1911. doi:10.1051/m2an/2017060

CrossRef Full Text | Google Scholar

Archie, G. E. (1942). The electrical resistivity log as an aid in determining some reservoir characteristics. Petroleum Trans. AIME 146, 54–62. doi:10.2118/942054-g

CrossRef Full Text | Google Scholar

Beale, J. T., and Majda, A. J. (1981). Rates of convergence for viscous splitting of the Navier-Stokes equations. Math. Comput. 37, 243–259. doi:10.1090/s0025-5718-1981-0628693-0

CrossRef Full Text | Google Scholar

Beaugendre, H., Huberson, S., and Mortazavi, I. (2012). Coupling particle sets of contours and streamline methods for solving convection problems. Appl. Math. Lett. 25, 11–19. doi:10.1016/j.aml.2011.06.031

CrossRef Full Text | Google Scholar

Bohn, H., McNeal, B., and O’Connor, G. (1980). Soil chemistry. Soil Chem. Soil Sci. 129, 389. doi:10.1097/00010694-198006000-00010

CrossRef Full Text | Google Scholar

Busenberg, E., and Plummer, L. N. (1986). A comparative study of the dissolution and crystal growth kinetics of calcite and aragonite. Stud. diagenesis 1578, 139–168. doi:10.3133/b1578

CrossRef Full Text | Google Scholar

Carman, P. C. (1937). Fluid flow through granular beds. Trans. Institution Chem. Eng. 15, 150–166.

Google Scholar

Chatelain, P., Curioni, A., Bergdorf, M., Rossinelli, D., Andreoni, W., and Koumoutsakos, P. (2008). Billion vortex particle direct numerical simulations of aircraft wakes. Comput. Methods Appl. Mech. Eng. 197, 1296–1304. doi:10.1016/j.cma.2007.11.016

CrossRef Full Text | Google Scholar

Chatelin, R., and Poncet, P. (2013). A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity. SIAM J. Sci. Comput. 35, B925–B949. doi:10.1137/120892921

CrossRef Full Text | Google Scholar

Chorin, A. J. (1973). Numerical study of slightly viscous flow. J. Fluid Mech. 57, 785–796. doi:10.1017/S0022112073002016

CrossRef Full Text | Google Scholar

Chou, L., Garrels, R. M., and Wollast, R. (1989). Comparative study of the kinetics and mechanisms of dissolution of carbonate minerals. Chem. Geol. 78, 269–282. doi:10.1016/0009-2541(89)90063-6

CrossRef Full Text | Google Scholar

Christoffersen, J., and Christoffersen, M. R. (1990). Kinetics of spiral growth of calcite crystals and determination of the absolute rate constant. J. Cryst. Growth 100, 203–211. doi:10.1016/0022-0248(90)90623-S

CrossRef Full Text | Google Scholar

Class, H., Ebigbo, A., Helmig, R., Dahle, H. K., Nordbotten, J. M., Celia, M. A., et al. (2009). A benchmark study on problems related to CO2 storage in geologic formations. Comput. Geosci. 13, 409–434. doi:10.1007/s10596-009-9146-x

CrossRef Full Text | Google Scholar

Coindreau, O., and Vignoles, G. L. (2005). Assessment of geometrical and transport properties of a fibrous c/c composite preform using x-ray computerized micro-tomography: Part i. image acquisition and geometrical properties. J. Mater. Res. 20, 2328–2339. doi:10.1557/jmr.2005.0311

CrossRef Full Text | Google Scholar

Cottet, G.-H., Balarac, G., and Coquerelle, M. (2009). Subgrid particle resolution for the turbulent transport of a passive scalar. Adv. Turbul. XII, 779–782. doi:10.1007/978-3-642-03085-7_188

CrossRef Full Text | Google Scholar

Cottet, G.-H., Etancelin, J.-M., Pérignon, F., and Picard, C. (2014). High order Semi-Lagrangian particle methods for transport equations: numerical analysis and implementation issues. ESAIM Math. Model. Numer. Analysis 48, 1029–1060. doi:10.1051/m2an/2014009

CrossRef Full Text | Google Scholar

Cottet, G.-H., and Koumoutsakos, P. (2000). Vortex methods: theory and practice. Cambridge University Press. doi:10.1017/CBO9780511526442

CrossRef Full Text | Google Scholar

Cottet, G.-H., Koumoutsakos, P., and Salihi, M. L. O. (2000). Vortex Methods with spatially varying cores. J. Comput. Phys. 162, 164–185. doi:10.1006/jcph.2000.6531

CrossRef Full Text | Google Scholar

Cottet, G.-H., and Poncet, P. (2004). Advances in direct numerical simulations of 3d wall-bounded flows by vortex-in-cell methods. J. Comput. Phys. 193, 136–158. doi:10.1016/j.jcp.2003.08.025

CrossRef Full Text | Google Scholar

Degond, P., and Mas-Gallic, S. (1989). The weighted particle method for convection-diffusion equations. Math. Comput. 53, 485–526. doi:10.2307/2008716

CrossRef Full Text | Google Scholar

Deng, H., Poonoosamy, J., and Molins, S. (2022). A reactive transport modeling perspective on the dynamics of interface-coupled dissolution-precipitation. Appl. Geochem. 137, 105207. doi:10.1016/j.apgeochem.2022.105207

CrossRef Full Text | Google Scholar

Dolbeau, R., Bodin, F., and de Verdière, G. C. (2013). “One OpenCL to rule them all?,” in 2013 IEEE 6th international workshop on multi-/many-core computing systems (MuCoCoS), 1–6. doi:10.1109/MuCoCoS.2013.6633603

CrossRef Full Text | Google Scholar

El Ossmani, M., and Poncet, P. (2010). Efficiency of multiscale hybrid grid-particle vortex methods. Multiscale Model. and Simul. 8, 1671–1690. doi:10.1137/090765006

CrossRef Full Text | Google Scholar

Etancelin, J.-M. (2014). Couplage de modèles, algorithmes multi-échelles et calcul hybride. Grenoble, France: Université de Grenoble. Ph.D. thesis.

Google Scholar

Etancelin, J.-M., Moonen, P., and Poncet, P. (2020). Improvement of remeshed Lagrangian methods for the simulation of dissolution processes at pore-scale. Adv. Water Resour. 146, 103780. doi:10.1016/j.advwatres.2020.103780

CrossRef Full Text | Google Scholar

[Dataset] Etancelin, J.-M., Mimeau, C., Keck, J.-B., Picard, C., Cottet, G.-H., Perignon, F., et al. (2022). HySoP: hybrid simulation with particles. Available at: https://hal.science/hal-04606382.

Google Scholar

Eyring, H. (1935). The activated complex in chemical reactions. J. Chem. Phys. 3, 107–115. doi:10.1063/1.1749604

CrossRef Full Text | Google Scholar

Faragó, I. (2008). A modified iterated operator splitting method. Appl. Math. Model. 32, 1542–1551. doi:10.1016/j.apm.2007.04.018

CrossRef Full Text | Google Scholar

Fazeli, H., Masoudi, M., Patel, R. A., Aagaard, P., and Hellevang, H. (2020). Pore-scale modeling of nucleation and growth in porous media. ACS Earth Space Chem. 4, 249–260. doi:10.1021/acsearthspacechem.9b00290

CrossRef Full Text | Google Scholar

Glover, P. W. J. (2016). Archie’s law – a reappraisal. Solid earth. 7, 1157–1169. doi:10.5194/se-7-1157-2016

CrossRef Full Text | Google Scholar

Hejlesen, M. M., Koumoutsakos, P., Leonard, A., and Walther, J. H. (2015). Iterative brinkman penalization for remeshed vortex methods. J. Comput. Phys. 280, 547–562. doi:10.1016/j.jcp.2014.09.029

CrossRef Full Text | Google Scholar

Hellevang, H., Pham, V. T., and Aagaard, P. (2013). Kinetic modelling of CO2–water–rock interactions. Int. J. Greenh. Gas Control 15, 3–15. doi:10.1016/j.ijggc.2013.01.027

CrossRef Full Text | Google Scholar

Hume, L., and Poncet, P. (2021). A velocity-vorticity method for highly viscous 3D flows with application to digital rock physics. J. Comput. Phys. 425, 109910. doi:10.1016/j.jcp.2020.109910

CrossRef Full Text | Google Scholar

Keck, J.-B. (2019). Numerical modelling and high performance computing for sediment flows. Theses.

Google Scholar

Kozeny, J. (1927). Umber kapillare Leitung des Wassers im Boden. Sitzungsber Akad. Wiss., Wien 136, 271–306.

Google Scholar

Lamichhane, B. P. (2013). A new Finite Element method for Darcy-Stokes-Brinkman equations. ISRN Comput. Math. 2013, 1–4. doi:10.1155/2013/798059

CrossRef Full Text | Google Scholar

Lasaga, A. C. (1981). Transition state theory. Rev. Mineral. 8.

Google Scholar

Lasaga, A. C. (1984). Chemical kinetics of water-rock interactions. J. Geophys. Res. Solid Earth 89, 4009–4025. doi:10.1029/JB089iB06p04009

CrossRef Full Text | Google Scholar

Lioliou, M. G., Paraskeva, C. A., Koutsoukos, P. G., and Payatakes, A. C. (2007). Heterogeneous nucleation and growth of calcium carbonate on calcite and quartz. J. Colloid Interface Sci. 308, 421–428. doi:10.1016/j.jcis.2006.12.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, H., Quintard, M., Debenest, G., and Laouafa, F. (2012). Properties of a diffuse interface model based on a porous medium theory for solid–liquid dissolution problems. Comput. Geosci. 16, 913–932. doi:10.1007/s10596-012-9295-1

CrossRef Full Text | Google Scholar

Maes, J., Soulaine, C., and Menke, H. P. (2022). Improved volume-of-solid formulations for micro-continuum simulation of mineral dissolution at the pore-scale. Front. Earth Sci. 10. doi:10.3389/feart.2022.917931

PubMed Abstract | CrossRef Full Text | Google Scholar

Magni, A., and Cottet, G.-H. (2012). Accurate, non-oscillatory, remeshing schemes for particle methods. J. Comput. Phys. 231, 152–172. doi:10.1016/j.jcp.2011.09.005

CrossRef Full Text | Google Scholar

Masoudi, M., Fazeli, H., Miri, R., and Hellevang, H. (2021). Pore scale modeling and evaluation of clogging behavior of salt crystal aggregates in CO2-rich phase during carbon storage. Int. J. Greenh. Gas Control 111, 103475. doi:10.1016/j.ijggc.2021.103475

CrossRef Full Text | Google Scholar

Masoudi, M., Nooraiepour, M., Deng, H., and Hellevang, H. (2024). Mineral precipitation and geometry alteration in porous structures: how to upscale variations in permeability–porosity relationship? Energy and Fuels 38, 9988–10001. doi:10.1021/acs.energyfuels.4c01432

CrossRef Full Text | Google Scholar

Menke, H. P., Bijeljic, B., Andrew, M. G., and Blunt, M. J. (2015). Dynamic three-dimensional pore-scale imaging of reaction in a carbonate at reservoir conditions. Environ. Sci. and Technol. 49, 4407–4414. doi:10.1021/es505789f

PubMed Abstract | CrossRef Full Text | Google Scholar

Molins, S., and Knabner, P. (2019). Multiscale approaches in reactive transport modeling. Rev. Mineralogy Geochem. 85, 27–48. doi:10.2138/rmg.2019.85.2

CrossRef Full Text | Google Scholar

Molins, S., Soulaine, C., Prasianakis, N. I., Abbasi, A., Poncet, P., Ladd, A. J. C., et al. (2021). Simulation of mineral dissolution at the pore scale with evolving fluid-solid interfaces: review of approaches and benchmark problem set. Comput. Geosci. 25, 1285–1318. doi:10.1007/s10596-019-09903-x

CrossRef Full Text | Google Scholar

Molins, S., Trebotich, D., Steefel, C. I., and Shen, C. (2012). An investigation of the effect of pore scale flow on average geochemical reaction rates using direct numerical simulation. Water Resour. Res. 48. doi:10.1029/2011WR011404

CrossRef Full Text | Google Scholar

Nielsen, A. E. (1984). Electrolyte crystal growth mechanisms. J. Cryst. Growth 67, 289–310. doi:10.1016/0022-0248(84)90189-1

CrossRef Full Text | Google Scholar

Nielsen, A. E., and Toft, J. M. (1984). Electrolyte crystal growth kinetics. J. Cryst. Growth 67, 278–288. doi:10.1016/0022-0248(84)90188-X

CrossRef Full Text | Google Scholar

Noiriel, C., and Renard, F. (2022). Four-dimensional X-ray micro-tomography imaging of dynamic processes in geosciences. Comptes Rendus. Géoscience 354, 255–280. doi:10.5802/crgeos.137

CrossRef Full Text | Google Scholar

Noiriel, C., Seigneur, N., Le Guern, P., and Lagneau, V. (2021). Geometry and mineral heterogeneity controls on precipitation in fractures: an X-ray micro-tomography and reactive transport modeling study. Adv. Water Resour. 152, 103916. doi:10.1016/j.advwatres.2021.103916

CrossRef Full Text | Google Scholar

Noiriel, C., and Soulaine, C. (2021). Pore-scale imaging and modelling of reactive flow in evolving porous media: tracking the dynamics of the fluid–rock interface. Transp. Porous Media 140, 181–213. doi:10.1007/s11242-021-01613-2

CrossRef Full Text | Google Scholar

Noiriel, C., Steefel, C. I., Yang, L., and Bernard, D. (2016). Effects of pore-scale precipitation on permeability and flow. Adv. Water Resour. 95, 125–137. doi:10.1016/j.advwatres.2015.11.013

CrossRef Full Text | Google Scholar

Nooraiepour, M., Masoudi, M., and Hellevang, H. (2021a). Probabilistic nucleation governs time, amount, and location of mineral precipitation and geometry evolution in the porous medium. Sci. Rep. 11, 16397. doi:10.1038/s41598-021-95237-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Nooraiepour, M., Masoudi, M., Shokri, N., and Hellevang, H. (2021b). Probabilistic nucleation and crystal growth in porous medium: new insights from calcium carbonate precipitation on primary and secondary substrates. ACS Omega 6, 28072–28083. doi:10.1021/acsomega.1c04147

PubMed Abstract | CrossRef Full Text | Google Scholar

Panga, M. K. R., Ziauddin, M., and Balakotaiah, V. (2005). Two-scale continuum model for simulation of wormholes in carbonate acidization. AIChE J. 51, 3231–3248. doi:10.1002/aic.10574

CrossRef Full Text | Google Scholar

Patel, R. A., Churakov, S. V., and Prasianakis, N. I. (2021). A multi-level pore scale reactive transport model for the investigation of combined leaching and carbonation of cement paste. Cem. Concr. Compos. 115, 103831. doi:10.1016/j.cemconcomp.2020.103831

CrossRef Full Text | Google Scholar

Payton, R. L., Sun, Y., Chiarella, D., and Kingdon, A. (2022). Pore scale numerical modelling of geological carbon storage through mineral trapping using true pore geometries. Transp. Porous Media 141, 667–693. doi:10.1007/s11242-021-01741-9

CrossRef Full Text | Google Scholar

Perez, S., Moonen, P., and Poncet, P. (2022). On the deviation of computed permeability induced by unresolved morphological features of the pore space. Transp. Porous Media 141, 151–184. doi:10.1007/s11242-021-01713-z

CrossRef Full Text | Google Scholar

Perez, S., and Poncet, P. (2024). Auto-weighted Bayesian physics-informed neural networks and robust estimations for multitask inverse problems in pore-scale imaging of dissolution. Comput. Geosci. 28, 1175–1215. doi:10.1007/s10596-024-10313-x

CrossRef Full Text | Google Scholar

Pham, V. T., Lu, P., Aagaard, P., Zhu, C., and Hellevang, H. (2011). On the potential of CO2–water–rock interactions for CO2 storage using a modified kinetic model. Int. J. Greenh. Gas Control 5, 1002–1015. doi:10.1016/j.ijggc.2010.12.002

CrossRef Full Text | Google Scholar

Plummer, L. N., and Busenberg, E. (1982). The solubilities of calcite, aragonite and vaterite in CO2-H2O solutions between 0 and 90°C, and an evaluation of the aqueous model for the system CaCO3-CO2-H2O. Geochimica Cosmochimica Acta 46, 1011–1040. doi:10.1016/0016-7037(82)90056-4

CrossRef Full Text | Google Scholar

Plummer, L. N., Wigley, T. M. L., and Parkhurst, D. L. (1978). The kinetics of calcite dissolution in CO2-water systems at 5 degrees to 60 degrees C and 0.0 to 1.0 atm CO2. Am. J. Sci. 278, 179–216. doi:10.2475/ajs.278.2.179

CrossRef Full Text | Google Scholar

Poncet, P. (2006). Finite difference stencils based on particle strength exchange schemes for improvement of vortex methods. J. Turbul. 7, N23. doi:10.1080/14685240600595586

CrossRef Full Text | Google Scholar

Poncet, P. (2007). Analysis of direct three-dimensional parabolic panel methods. SIAM J. Numer. Analysis 45, 2259–2297. doi:10.1137/050625849

CrossRef Full Text | Google Scholar

Poonoosamy, J., Obaied, A., Deissmann, G., Prasianakis, N. I., Kindelmann, M., Wollenhaupt, B., et al. (2023). Microfluidic investigation of pore-size dependency of barite nucleation. Commun. Chem. 6, 250. doi:10.1038/s42004-023-01049-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Quintard, M., and Whitaker, S. (1988). Two-phase flow in heterogeneous porous media: the method of large-scale averaging. Transp. Porous Media 3, 357–413. doi:10.1007/BF00233177

CrossRef Full Text | Google Scholar

Quintard, M., and Whitaker, S. (1993). Transport in ordered and disordered porous media: volume-averaged equations, closure problems, and comparison with experiment. Chem. Eng. Sci. 48, 2537–2564. doi:10.1016/0009-2509(93)80266-S

CrossRef Full Text | Google Scholar

Sanchez, D., Hume, L., Chatelin, R., and Poncet, P. (2019). Analysis of the 3D non-linear Stokes problem coupled to transport-diffusion for shear-thinning heterogeneous microscale flows, applications to digital rock physics and mucociliary clearance. ESAIM Math. Model. Numer. Analysis 53, 1083–1124. doi:10.1051/m2an/2019013

CrossRef Full Text | Google Scholar

Schrader, B., Reboux, S., and Sbalzarini, I. F. (2010). Discretization correction of general integral PSE Operators for particle methods. J. Comput. Phys. 229, 4159–4182. doi:10.1016/j.jcp.2010.02.004

CrossRef Full Text | Google Scholar

[Dataset] Sheppard, A., and Prodanovic, M. (2015). Network generation comparison forum. doi:10.17612/P7059V

CrossRef Full Text | Google Scholar

Siena, M., Bussetti, G., Recalcati, C., Riva, M., Duò, L., and Guadagnini, A. (2021). Statistical characterization of heterogeneous dissolution rates of calcite from in situ and real-time AFM imaging. Transp. Porous Media 140, 291–312. doi:10.1007/s11242-021-01624-z

CrossRef Full Text | Google Scholar

Soulaine, C., Roman, S., Kovscek, A., and Tchelepi, H. A. (2017). Mineral dissolution and wormholing from a pore-scale perspective. J. Fluid Mech. 827, 457–483. doi:10.1017/jfm.2017.499

CrossRef Full Text | Google Scholar

Soulaine, C., Roman, S., Kovscek, A., and Tchelepi, H. A. (2018). Pore-scale modelling of multiphase reactive flow: application to mineral dissolution with production of CO₂. J. Fluid Mech. 855, 616–645. doi:10.1017/jfm.2018.655

CrossRef Full Text | Google Scholar

Starchenko, V. (2022). Pore-scale modeling of mineral growth and nucleation in reactive flow. Front. Water 3. doi:10.3389/frwa.2021.800944

CrossRef Full Text | Google Scholar

Steefel, C. I., Appelo, C. A. J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., et al. (2015). Reactive transport codes for subsurface environmental simulation. Comput. Geosci. 19, 445–478. doi:10.1007/s10596-014-9443-x

CrossRef Full Text | Google Scholar

Steefel, C. I., and Lasaga, A. C. (1990). “Evolution of dissolution patterns,” in Chemical Modeling of aqueous systems II (American chemical society), vol. 416 of ACS symposium series, 212–225. doi:10.1021/bk-1990-0416.ch016

CrossRef Full Text | Google Scholar

Steefel, C. I., and Lasaga, A. C. (1994). A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems. Am. J. Sci. 294, 529–592. doi:10.2475/ajs.294.5.529

CrossRef Full Text | Google Scholar

Stone, J. E., Gohara, D., and Shi, G. (2010). OpenCL: a parallel programming standard for heterogeneous computing systems. Comput. Sci. and Eng. 12, 66–73. doi:10.1109/MCSE.2010.69

PubMed Abstract | CrossRef Full Text | Google Scholar

Strang, G. (1968). On the construction and comparison of difference schemes. SIAM J. Numer. Analysis 5, 506–517. doi:10.1137/0705041

CrossRef Full Text | Google Scholar

Tartakovsky, A. M., Meakin, P., Scheibe, T. D., and Wood, B. D. (2007). A smoothed particle hydrodynamics model for reactive transport and mineral precipitation in porous and fractured porous media. Water Resour. Res. 43. doi:10.1029/2005WR004770

PubMed Abstract | CrossRef Full Text | Google Scholar

Tong, L.-y., Cai, Y., and Liu, Q.-f. (2024). Carbonation modelling of hardened cementitious materials considering pore structure characteristics: a review. J. Build. Eng. 96, 110547. doi:10.1016/j.jobe.2024.110547

CrossRef Full Text | Google Scholar

Turk, M. C., Shi, X., Gonyer, D. A. J., and Roy, D. (2015). Chemical and mechanical aspects of a Co-Cu planarization scheme based on an alkaline slurry formulation. ECS J. Solid State Sci. Technol. 5, P88–P99. doi:10.1149/2.0271602jss

CrossRef Full Text | Google Scholar

Varzina, A., Cizer, O., Yu, L., Liu, S., Jacques, D., and Perko, J. (2020). A new concept for pore-scale precipitation-dissolution modelling in a lattice Boltzmann framework – application to portlandite carbonation. Appl. Geochem. 123, 104786. doi:10.1016/j.apgeochem.2020.104786

CrossRef Full Text | Google Scholar

von Wolff, L., Weinhardt, F., Class, H., Hommel, J., and Rohde, C. (2021). Investigation of crystal growth in enzymatically induced calcite precipitation by micro-fluidic experimental methods and comparison with mathematical modeling. Transp. Porous Media 137, 327–343. doi:10.1007/s11242-021-01560-y

CrossRef Full Text | Google Scholar

Wakao, N., and Smith, J. M. (1962). Diffusion in catalyst pellets. Chem. Eng. Sci. 17, 825–834. doi:10.1016/0009-2509(62)87015-8

CrossRef Full Text | Google Scholar

Wigand, M., Carey, J. W., Schütt, H., Spangenberg, E., and Erzinger, J. (2008). Geochemical effects of CO2 sequestration in sandstones under simulated in situ conditions of deep saline aquifers. Appl. Geochem. 23, 2735–2745. doi:10.1016/j.apgeochem.2008.06.006

CrossRef Full Text | Google Scholar

Wolthers, M., Nehrke, G., Gustafsson, J. P., and Van Cappellen, P. (2012). Calcite growth kinetics: modeling the effect of solution stoichiometry. Geochimica Cosmochimica Acta 77, 121–134. doi:10.1016/j.gca.2011.11.003

CrossRef Full Text | Google Scholar

Wood, B. D., Radakovich, K., and Golfier, F. (2007). Effective reaction at a fluid-solid interface: applications to biotransformation in porous media. Adv. water Resour. 30, 1630–1647. doi:10.1016/j.advwatres.2006.05.032

CrossRef Full Text | Google Scholar

Yang, F., Stack, A. G., and Starchenko, V. (2021). Micro-continuum approach for mineral precipitation. Sci. Rep. 11, 3495. doi:10.1038/s41598-021-82807-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoreo, J. J. D., and Vekilov, P. G. (2003). Principles of crystal nucleation and growth. Berlin: De Gruyter. doi:10.1515/9781501509346-008

CrossRef Full Text | Google Scholar