2. Mathematical formulation

We consider a two-dimensional porous medium with square geometry composed mainly of carbonate rocks, or $\text{CaCO}_{3}$ (figure 1). The porous domain has impermeable boundaries at the top and bottom, and we assume periodicity in the lateral direction. The height of the square domain is denoted $H$ . There are two participating solutions in the system; solution 1 consists of brine saturated with $\text{CO}_{2}$ , and it enters the domain through the top boundary; solution 2 is the ambient brine that initially fully saturates the porous medium. An important assumption we make in this formulation is that both solutions are initially at chemical equilibrium with the carbonate formation. This assumption is justified because (i) the pH for solution 1 is around 4, which indicates fast reaction with the carbonate; (ii) solution 2 is the ambient brine, which has already reached chemical equilibrium with the carbonate during its geologic residence time. This assumption distinguishes our formulation from the classical kinetic reaction behaviour. Under this set-up, neither solution 1 nor solution 2 alone will induce reactions; rather it is the mixing of the two solutions, which disturbs local chemical equilibrium, that triggers geochemical reactions and rock dissolution. We assume Boussinesq, incompressible Darcy flow through porous media (Nield & Bejan Reference Nield and Bejan2006):

(2.1) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad \boldsymbol{u}=-\frac{k({\it\phi})}{{\it\mu}}(\boldsymbol{{\rm\nabla}}P-{\it\rho}g\boldsymbol{z}),\end{eqnarray}$$

where $\boldsymbol{u}$ is the Darcy velocity, $k$ is permeability, ${\it\mu}$ is dynamic viscosity, $P$ is pressure, ${\it\rho}$ is fluid density, $g$ is gravitational acceleration and ${\it\phi}$ is porosity. We adopt a simple cubic law for the permeability–porosity relationship, $k=k_{0}({\it\phi}/{\it\phi}_{0})^{3}$ , where $k_{0}$ and ${\it\phi}_{0}$ are the initial permeability and porosity, respectively. This relationship captures the evolution in permeability induced by porosity increase observed in laboratory experiments (Carroll et al. Reference Carroll, Hao, Smith and Sholokhova2013).

Figure 1. The problem is set up in a two-dimensional square homogeneous porous medium composed mainly of carbonate rocks. The domain is impermeable at the top and bottom. Solution 1, the equilibrated $\text{CO}_{2}$ -rich brine, enters via diffusion at the top of the formation, which is initially saturated with solution 2, the equilibrated aquifer brine.

We denote the density of solution 1 as ${\it\rho}_{1}$ and that of solution 2 as ${\it\rho}_{2}$ , where ${\it\rho}_{1}>{\it\rho}_{2}$ , and the density difference as ${\rm\Delta}{\it\rho}={\it\rho}_{1}-{\it\rho}_{2}>0$ . We define the mixing ratio ${\it\alpha}$ as the volumetric ratio of solution 1 in the mixture:

(2.2) $$\begin{eqnarray}{\it\alpha}(x,z,t)=\frac{V_{solution\,1}}{V_{total}}.\end{eqnarray}$$

The density of the fluid mixture, ${\it\rho}$ , increases linearly with the mixing ratio:

(2.3) $$\begin{eqnarray}{\it\rho}={\it\rho}_{2}+{\rm\Delta}{\it\rho}\,{\it\alpha}.\end{eqnarray}$$

When the two solutions mix at ${\it\alpha}$ , chemical equilibrium is temporarily disturbed locally and the mixture undergoes a series of geochemical reactions to reach a new equilibrium with calcite. While this process can be complicated by the participation of other co-existing minerals such as magnesite ( $\text{MgCO}_{3}$ ) or gypsum ( $\text{CaSO}_{4}\cdot \text{H}_{2}\text{O}$ ), the dominant effect can be captured with the following four reactions (Sanford & Konikow Reference Sanford and Konikow1989; De Simoni et al. Reference De Simoni, Sánchez-Vila, Carrera and Saaltink2007):

(2.4) $$\begin{eqnarray}\displaystyle & \text{HCO}_{3}^{-}\rightleftharpoons \text{CO}_{2}+\text{H}_{2}\text{O}-\text{H}^{+} & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \text{CO}_{3}^{2-}\rightleftharpoons \text{CO}_{2}+\text{H}_{2}\text{O}-2\text{H}^{+} & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \text{Ca}^{2+}\rightleftharpoons \text{CaCO}_{3}-\text{CO}_{2}-\text{H}_{2}\text{O}+2\text{H}^{+} & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \text{OH}^{-}\rightleftharpoons \text{H}_{2}\text{O}-\text{H}^{+}. & \displaystyle\end{eqnarray}$$

Eight chemical species participate in these reactions, namely calcium carbonate ( $\text{CaCO}_{3}$ ), water ( $\text{H}_{2}\text{O}$ ), carbon dioxide ( $\text{CO}_{2}$ ), hydrogen ion ( $\text{H}^{+}$ ), bicarbonate ( $\text{HCO}_{3}^{-}$ ), carbonate ( $\text{CO}_{3}^{2-}$ ), calcium ion ( $\text{Ca}^{2+}$ ) and hydroxide ( $\text{OH}^{-}$ ). Adopting the approach of De Simoni et al. (Reference De Simoni, Sánchez-Vila, Carrera and Saaltink2007) discussed in § 1.1, we track the mixing ratio  ${\it\alpha}$ , which is a conserved quantity in the system equivalent to ${\it\tau}$ in (1.3), using the advection–diffusion equation:

(2.8) $$\begin{eqnarray}{\it\phi}\frac{\partial {\it\alpha}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{u}{\it\alpha})-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\phi}D\boldsymbol{{\rm\nabla}}{\it\alpha})=0,\end{eqnarray}$$

where $D$ is the diffusion/dispersion coefficient, taken here as constant. We assume here that $D$ is the same for all chemical species in the system and is also the same as that of the fluid mixture. Using (1.4), we can then explicitly calculate the rate of calcite dissolution, which according to (2.6) is also the production rate for $\text{Ca}^{2+}$ :

(2.9) $$\begin{eqnarray}r=2{\it\phi}F({\it\alpha})\left(\boldsymbol{{\rm\nabla}}{\it\alpha}\boldsymbol{\cdot }D\boldsymbol{{\rm\nabla}}{\it\alpha}\right),\quad F({\it\alpha})=\frac{\partial ^{2}[\text{Ca}^{2+}]}{\partial {\it\alpha}^{2}},\end{eqnarray}$$

where $F({\it\alpha})$ is the speciation term associated with $\text{Ca}^{2+}$ that we obtain a priori with a geochemical code such as PHREEQC (see appendix B). Equation (2.9) states that the local rate of porosity change is controlled by two mechanisms: (i) chemical speciation of calcium ion, $F({\it\alpha})$ , which measures the vigor of the rock dissolution reaction in (2.6); and (ii) strength of fluid mixing, which is measured in terms of scalar dissipation rate of ${\it\alpha}$ . Finally, we update porosity locally using the calcite dissolution rate $r$ :

(2.10) $$\begin{eqnarray}c_{s}\frac{\partial {\it\phi}}{\partial t}=r{\it\Theta}({\it\phi}_{1}-{\it\phi}),\end{eqnarray}$$

where $c_{s}$ is the molar density of calcite, and ${\it\Theta}$ is the Heaviside step function that limits the increase in porosity to a maximum value ${\it\phi}_{1}$ , which accounts for the presence of residual minerals that are non-reactive in the system. Once ${\it\phi}$ reaches ${\it\phi}_{1}$ locally, the assumption of chemical equilibrium fails since there is no more calcite to react with the unequilibriated mixtures; nevertheless, we assume equilibrium for solutions that pass through the inert area. We keep this special case in mind and address its limitations when we interpret the results in § 3.2.

We choose as characteristic velocity the speed at which a fluid parcel sinks in the porous medium: $U_{c}=k_{0}{\rm\Delta}{\it\rho}g/{\it\mu}$ . The natural length scale is the length over which diffusion and advection are balanced: $L={\it\phi}D/U_{c}$ . We choose to rescale the problem with this intrinsic length scale $L$ as it is more relevant than domain height ( $H$ ) for characterizing local quantity changes. We set the other characteristic quantities as

(2.11a−f ) $$\begin{eqnarray}{\it\rho}_{c}={\rm\Delta}{\it\rho},\quad P_{c}={\rm\Delta}{\it\rho}gL,\quad T_{c}=\frac{{\it\phi}_{0}L}{U_{c}},\quad {\it\phi}_{c}={\it\phi}_{0},\quad F_{c}=\max (|F({\it\alpha})|),\quad r_{c}=\frac{2{\it\phi}_{0}F_{c}D}{L^{2}}.\end{eqnarray}$$

Applying these scales, we obtain three dimensionless parameters: the dimensionless domain height $\text{H}=H/L$ , the maximum relative porosity increase $R_{{\it\phi}}={\it\phi}_{1}/{\it\phi}_{0}$ , and the dissolution Damkhöler number:

(2.12) $$\begin{eqnarray}\mathit{Da}=\frac{r_{c}T_{c}}{c_{s}{\it\phi}_{0}}=\frac{2F_{c}}{c_{s}}.\end{eqnarray}$$

Since we use the intrinsic length scale here, the Rayleigh number ( $\mathit{Ra}=U_{c}L/D$ ) is identically equal to one, and $\text{H}$ plays the role of the traditional $\mathit{Ra}$ (Hidalgo & Carrera Reference Hidalgo and Carrera2009). The dissolution Damkhöler number measures the competition between rock dissolution rate and solute transport velocity. In dimensionless form, the governing equations read

(2.13) $$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad \boldsymbol{u}=-{\it\phi}^{3}(\boldsymbol{{\rm\nabla}}P^{\prime }-{\it\alpha}\boldsymbol{z}), & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\phi}\frac{\partial {\it\alpha}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{u}{\it\alpha})-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\phi}\boldsymbol{{\rm\nabla}}{\it\alpha})=0, & \displaystyle\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\phi}}{\partial t}=\mathit{Da}~r{\it\Theta}(R_{{\it\phi}}-{\it\phi}/{\it\phi}_{0}),\quad r={\it\phi}F({\it\alpha})(\boldsymbol{{\rm\nabla}}{\it\alpha}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\alpha}). & \displaystyle\end{eqnarray}$$

Here, $P^{\prime }$ is the dimensionless pressure with respect to a hydrostatic datum posed by ${\it\rho}_{2}$ . The boundary and initial conditions are

(2.16a,b ) $$\begin{eqnarray}\displaystyle {\it\alpha}(x,z=0)=1,\quad \left.\frac{\partial {\it\alpha}}{\partial z}\right|_{x,z=\text{H}}=0,\end{eqnarray}$$

(2.17a,b ) $$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}(x,z=0)=0,\quad \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}(x,z=\text{H})=0,\end{eqnarray}$$

(2.18a,b ) $$\begin{eqnarray}{\it\phi}(t=0)=1,\quad {\it\alpha}(t=0)=0.\end{eqnarray}$$

We adopt a stream function–vorticity formulation for (2.13) (Tan & Homsy Reference Tan and Homsy1988; Riaz & Meiburg Reference Riaz and Meiburg2003; Fu et al. Reference Fu, Cueto-Felgueroso and Juanes2013) and solve for the stream functions with an eighth-order finite difference scheme, implemented with a fast Poisson solver (Swarztrauber Reference Swarztrauber1977). To evaluate the derivatives used in the transport equation (2.14) and consumption rate (2.15b ), we use sixth-order compact finite differences in the vertical direction (Lele Reference Lele1992) and a Fourier discretization along the horizontal direction, which we assume to be periodic. We integrate (2.14) and (2.15a ) sequentially in time using a third-order Runge–Kutta scheme with dynamic time stepping (Ruith & Meiburg Reference Ruith and Meiburg2000). We accelerate the onset of the gravitational instability by perturbing the initial mixing ratios at the top boundary with small white noise (an uncorrelated Gaussian random function).

3. Dissolution regimes

In this section, we describe dissolution patterns as a result of the reaction with calcite, and we illustrate the various dissolution regimes. Specifically, we focus on $\langle {\rm\Delta}{\it\phi}\rangle$ , the domain-averaged porosity increase, and $\text{d}\langle {\it\alpha}\rangle /\text{d}t$ , the time derivative of the domain-averaged pore–volume ratio of solution 1 (the equilibrated $\text{CO}_{2}$ -rich brine), or cumulative mass of solution 1 in the domain. We solve the governing equations for Damkhöler numbers up to 20 and values of $\text{H}$ up to 6000 on a grid of $2500^{2}$ cells. Convergence tests show that the required grid resolution increases superlinearly with the product of $\text{H}$ and $\mathit{Da}$ since both parameters have a combined effect on generating small-scale details.

Figure 2. Snapshots of mixing ratio (a) and porosity percentage increase (b) for a simulation with $\text{H}=6000$ , $\mathit{Da}=20$ and $R_{{\it\phi}}=2$ at time $t=0.6\text{H}$ (top), 4H (middle) and 15H (bottom).

3.1. Mixing-controlled dissolution patterns

The diffusive layer of solution-1-rich mixture near the top boundary becomes gravitationally unstable after some initial onset time, resulting in finger-like structures (figure 2, top left) that migrate downwards, stripping away the dense fluid (figure 2, middle left). Fingering dynamics and scaling relations have been studied extensively (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Hassanzadeh et al. Reference Hassanzadeh, Pooladi-Darvish and Keith2007; Hidalgo & Carrera Reference Hidalgo and Carrera2009; Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012; Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2013; Slim et al. Reference Slim, Bandi, Miller and Mahadevan2013; Slim Reference Slim2014) and will not be the focus of this work. The density-driven flow leads to mixing of fluids and triggers rock dissolution reactions which alter the permeability field. Initially, these reactions chisel out patterns that follow closely the interface of the fluid fingers, leaving the interior of these fingering channels almost unaffected (figure 2, top and middle right). The ‘hollowness’ of the dissolution patterns can be understood as the direct consequence of a mixing-controlled dissolution rate (1.4). The rate of fluid mixing, as measured by the scalar dissipation rate, is strongest along finger interfaces where the concentration gradient is large (Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012) and is weak within the core of the fingers, where concentration gradients are small. Over time, the rock dissolution patterns evolve in two ways.

1. (a) Elongation towards the bottom: new fingers are continuously born at the top boundary; the cascade of fingers of different ages interact through a merging process where young fingers merge into old ones, small ones into bigger ones, allowing the surviving fingers to sink and coarsen laterally in size (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Hassanzadeh et al. Reference Hassanzadeh, Pooladi-Darvish and Keith2007; Pau et al. Reference Pau, Bell, Pruess, Almgren, Lijewski and Zhang2010; Slim Reference Slim2014). As a result, finger-like dissolution channels continue to extend downwards and expand laterally (figure 2, bottom), following the locations of finger interfaces where mixing is strong.

2. (b) Focusing near the top: while dissolution channels elongate towards the bottom due to downward migration of fingers, regions of high rock dissolution, or ‘dissolution hubs’ start to appear near the top boundary (figure 2, bottom). These hubs develop as a result of repeated events of newly born fingers carrying ${\it\alpha}\approx 1$ meeting the upwelling plume of ${\it\alpha}\approx 0$ , creating high concentration gradients that result in strong mixing and reaction.

Figure 3. Volume-averaged measures of (a) percentage increase in ${\it\phi}$ (b) chemical speciation and (c) scalar dissipation rate over time for a simulation with $\text{H}=6000$ . The transition in curve colour from light orange to dark indicates increase in $\mathit{Da}$ . The dark circles in (a) and (b) indicate the shutdown time $t_{sd}^{{\rm\Delta}{\it\phi}}$ and $t_{sd}^{F({\it\alpha})}$ , respectively, as defined in § 3.2. (d) Cumulative (volume-averaged) porosity increase plotted against $\mathit{Da}$ for different values of $\mathit{H}$ (square, triangle and circle). (e,f) Plots of $t_{sd}^{{\rm\Delta}{\it\phi}}$ and $t_{sd}^{F({\it\alpha})}$ against $\mathit{Da}$ for different values of $\mathit{H}$ .

3.2. Speciation-controlled reaction shutdown

The overall dissolution pattern is controlled by fluid mixing, as evidenced by the spatial correlation between the permeability and mixing fields. However, the magnitude of permeability increase, or how much reaction occurs when the two fluids mix, is set by the shape of the speciation term $F({\it\alpha})$ . As shown in figure 10 in appendix B, $F({\it\alpha})$ for our geochemical system is a highly nonlinear function, showing that reaction is favoured in the lower range of mixing ratios ( ${\it\alpha}<0.1$ ); when ${\it\alpha}>0.2$ , $F({\it\alpha})$ becomes so small that regardless of the strength of local mixing, almost no reaction will take place. The behaviour of this speciation curve is crucial for explaining the system’s evolution.

For a simulation with $\text{H}=6000$ , $\mathit{Da}=20$ , we track the domain-averaged porosity increase $\langle {\rm\Delta}{\it\phi}\rangle$ with time, and find that after a period of monotonic growth, $\langle {\rm\Delta}{\it\phi}\rangle$ reaches a plateau, indicating a shutdown of dissolution reactions globally (figure 3 a). We compute the time for dissolution shutdown, $t_{sd}^{{\rm\Delta}{\it\phi}}$ , as the time when $\langle {\rm\Delta}{\it\phi}\rangle$ reaches 95 % of the plateau value, denoted $\langle {\rm\Delta}{\it\phi}\rangle _{sd}$ ; the shutdown times are shown as dark circles in figure 3(a). We next track the domain-averaged speciation term over time and find that it also flattens out towards zero after some period of time. We compute the corresponding time for speciation shutdown, $t_{sd}^{F({\it\alpha})}$ , shown as dark circles in figure 3(b), as the time when $\langle F({\it\alpha})\rangle$ falls below 0.02. Shutdown time for dissolution corresponds closely to that of speciation (figure 3 a,b, dark circles) and this shutdown happens much earlier than the decaying of global degree of mixing (figure 3 c). In other words, while reaction-inducing fluid mixing is still active in the system, the fast decay of chemical speciation has caused the dissolution reactions to stop at a much earlier time. We refer to this as a speciation-controlled reaction shutdown. The concept of convective shutdown was coined in the context of the one-sided convection system to describe the effect of domain saturation on the mass flux (Hewitt et al. Reference Hewitt, Neufeld and Lister2013). Here, we find a shutdown in reaction caused by the decay of chemical speciation, not the loss in strength of the convective instability.

We extend the above analysis to a range of values of $\text{H}$ and $\mathit{Da}$ , and find that the plateau value $\langle {\rm\Delta}{\it\phi}\rangle _{sd}$ scales linearly with $\mathit{Da}$ (figure 3 d):

(3.1) $$\begin{eqnarray}\langle {\rm\Delta}{\it\phi}\rangle _{sd}\sim \mathit{Da},\end{eqnarray}$$

which is a consequence of the role that $\mathit{Da}$ plays in the porosity-update equation (2.15a ). The shutdown time for dissolution and speciation both decrease with increasing $\mathit{Da}$ (figure 3 e, f). This is because as $\mathit{Da}$ increases, the permeability near the top boundary increases faster, allowing fluid to enter and saturate the domain faster, and leading to an earlier speciation shutdown (see further discussion in § 4). We find that $\text{H}$ , on the other hand, has little impact on $\langle {\rm\Delta}{\it\phi}\rangle _{sd}$ and $t_{sd}^{{\rm\Delta}{\it\phi}}$ , an observation that is consistent with the fact that the rate of fluid mixing is independent of $\text{H}$ (Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012).

We analyse the time series of the porosity field and find that the maximum porosity ${\it\phi}_{1}$ is reached locally for simulations with large $\mathit{Da}$ and $\text{H}$ after some initial period of dissolution and this critical time arrives later for smaller $\mathit{Da}$ and $\text{H}$ . This indicates that the special case of non-equilibrium state as discussed in § 2 is indeed present in some of our simulations. For small enough $\mathit{Da}$ and $\text{H}$ (e.g.  $\text{H}=2000$ with any $\mathit{Da}$ ; $\text{H}=4000$ with $\mathit{Da}=1$ ; $\text{H}=6000$ with $\mathit{Da}=1$ ), however, ${\it\phi}_{1}$ is not reached anywhere. When ${\it\phi}_{1}$ is reached, we find that the spatial spread of the non-equilibrium state is focused near the top boundary and accounts for a small portion of the domain, between 0.17 and 0.69 %. This suggests that non-equilibrium effects are not critical in our examples. However, such non-equilibrium effect could be more important for simulations with large Damkhöler numbers as the maximum ${\it\phi}$ could be reached early on in the simulation and impact larger areas.

4. Impact on the macroscopic mass exchange rate

The macroscopic mass exchange rate measures how fast solution 1, the equilibrated $\text{CO}_{2}$ -rich brine at the top, enters the formation and mixes with the resident brine. It is a quantity of special interest in the context of geologic carbon sequestration as it determines the effectiveness of solubility trapping. Previous studies that do not account for geochemical reactions (Hassanzadeh et al. Reference Hassanzadeh, Pooladi-Darvish and Keith2007; Pau et al. Reference Pau, Bell, Pruess, Almgren, Lijewski and Zhang2010; Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012; Hewitt et al. Reference Hewitt, Neufeld and Lister2013; Slim et al. Reference Slim, Bandi, Miller and Mahadevan2013; Slim Reference Slim2014) have shown that the exchange rate, often referred to as the $\text{CO}_{2}$ flux, reaches a plateau after the onset of instability and remains constant until domain saturation starts to impact the top boundary, at which point the flux decays (figure 4, $\mathit{Da}=0$ ). Further, it has been shown that this constant flux is independent of the domain size $\text{H}$ (Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012; Slim Reference Slim2014). We address here how these observations change in the presence of rock dissolution reactions.

Figure 4. Volume-averaged measure of mass exchange rate of solution 1 for increasing Damkhöler numbers with $\text{H}=6000$ (a) and for three values of $\text{H}$ with $\mathit{Da}=20$ (b). The coarsest dashed curves in both figures correspond to the case with $\text{H}=6000$ and without reaction ( $\mathit{Da}=0$ ).

We define the mass exchange rate as the volumetric flux of solution 1, $\text{d}\langle {\it\alpha}\rangle /\text{d}t$ , which we calculate by taking numerical derivative of $\langle {\it\alpha}\rangle$ with respect to $t$ . For simulations with fixed $\text{H}$ (e.g.  $\text{H}=6000$ ) and $\mathit{Da}>1$ , we observe significant differences with respect to the non-reactive case. Instead of reaching a plateau after the onset of convection, the flux continues to grow with time as a result of increasing permeability from rock dissolution near the top (figure 4 a); the rate of flux increase grows with $\mathit{Da}$ . When $\mathit{Da}$ is sufficiently large (e.g.  $\mathit{Da}=20$ ), the increase in $\text{d}\langle {\it\alpha}\rangle /\text{d}t$ hits a ceiling early on (figure 4 a, $\mathit{Da}=20$ ); this is because the relative porosity increase $\langle {\rm\Delta}{\it\phi}\rangle$ at the top has reached $R_{{\it\phi}}$ and cannot increase further. In all cases, with or without reaction, the flux ultimately decays due to the effect of domain saturation and, as shown in figure 4, the flux decays at earlier times as $\mathit{Da}$ increases. This is, again, because an increase in permeability leads to a more convective system that can reach saturation faster. These observations are largely independent of the values of $\text{H}$ (figure 4 b), similar to what we observe in § 3.2.

5. Dissolution patterns in three dimensions

We simulate (2.13)–(2.15) with the same boundary and initial conditions in three dimensions, for a case with $\text{H}=2500$ , $\mathit{Da}=5$ and $R_{{\it\phi}}=2$ with a grid resolution of $368^{3}$ cells. The numerical scheme is the same as that for two dimensions, where we use Fourier discretization in the $x$ and $y$ directions and compact finite differences in the $z$ direction. Similar to the 2D results, the reaction shuts down at around $t/\text{H}=8$ . We present the final morphology (at $t/\text{H}=14$ ) of dissolved rock in figure 5. As the 3D columnar fingers move downwards (see more details in Fu et al. Reference Fu, Cueto-Felgueroso and Juanes2013) they chisel out columnar-shaped regions along the finger–brine interfaces (figure 5 a) as a result of the high local concentration gradients. Underneath the top boundary layer, the interiors of fingers are barely affected by reaction (figure 7 b–f, black–white frames). These empty columns are analogous to the hollow channels observed in two dimensions (see § 3.1). At the top boundary we observe the same dissolution hubs as seen in two dimensions (§ 3.1); here, they take the shape of upside-down dunes, scattered and hanging from the top surface (figure 5 b).

Figure 5. Dissolution patterns in three dimensions: for a simulation of $\text{H}=2500,~\mathit{Da}=5$ and $R_{{\it\phi}}=2$ , at $t/\text{H}=14.0$ , this figure shows surface contours for (a) ${\rm\Delta}{\it\phi}=4\,\%$ and (b) ${\rm\Delta}{\it\phi}=8\,\%$ .

While 3D visualization gives insight into the morphology of the rock dissolution patterns, it provides limited quantitative details. To obtain more information about the interior structure of the dissolution morphology, we plot horizontal slices of ${\rm\Delta}{\it\phi}/{\it\phi}_{0}$ at six depths: $z/\text{H}=0.01$ (top boundary layer), $z/\text{H}=0.14,027,0.54$ and $0.82$ (the bulk) and $z/\text{H}=1.0$ (bottom). As fingers start to arrive at a given depth, circular rings of dissolved rock emerge within the horizontal plane, reflecting the shape of fingering plumes at those depths. With the exception of the top boundary layer ( $z/\text{H}=0.01$ ), the dissolved rings continue to expand laterally to form a tessellation of polygons. The polygons are ‘hollow’ inside and separated from each other by a clearly defined polygonal network of undissolved edges (figure 7 b–f, black and white frames). Once the tessellation has evolved into its full form, usually within ${\rm\Delta}t=2500$ after the fingers invade that depth, the pattern ‘freezes’ in time due to a drastic reduction in dissolution activity.

We track the dissolution rate at each depth as a function of time, denoted here as $\text{d}\langle {\it\phi}\rangle _{z}/\text{d}t$ , where $\langle \cdot \rangle _{z}$ is the average across the $x{-}y$ plane, and observe that the dissolution rates peak at around the time that dissolution structures ‘freeze’ at that depth (figure 6). Similar to earlier observations in two dimensions (§ 3.2), the ‘freezing’ of rock dissolution structures is a consequence of a speciation-controlled reaction shutdown: once fingers arrive at a given layer, dissolution reactions take place for a short period of time; once the layer becomes weakly saturated with ${\it\alpha}$ (e.g.  $\langle {\it\alpha}\rangle _{z}>0.1$ ), the reactions stop. Figure 6 also clearly demonstrates the temporal dynamics of the shutdown process: instead of occurring all at once, the shutdown takes place at later times for deeper layers. Overall, this provides us with a more mechanistic description on how these tessellation patterns evolve: as fingers arrive at each layer, they imprint their ‘images’ through mineral dissolution onto the porosity field; once the reaction shuts down, the dissolution pattern becomes immune to further alterations by fingering dynamics.

Figure 6. Dissolution rate (laterally averaged) at various depths over time. Dots (cyan-coloured online) correspond to the ‘freezing’ events, which is around the time that the tessellation pattern at that depth adopts its final configuration.

5.1. Coarsening of dissolution patterns

Finger coarsening in non-reactive convective mixing has been observed and quantified in both experimental studies (Backhaus et al. Reference Backhaus, Turitsyn and Ecke2011; MacMinn & Juanes Reference MacMinn and Juanes2013) and simulations (Fu et al. Reference Fu, Cueto-Felgueroso and Juanes2013; Hewitt et al. Reference Hewitt, Neufeld and Lister2013; Slim et al. Reference Slim, Bandi, Miller and Mahadevan2013). The dynamics of coarsening are different in the boundary layer and within the bulk. Near the top boundary, finger roots initially coarsen through a series of merging events but then reaches a statistical steady state, where the finger root spacing stands at a quasi-steady value (Fu et al. Reference Fu, Cueto-Felgueroso and Juanes2013; MacMinn & Juanes Reference MacMinn and Juanes2013). On the other hand, finger coarsening in the bulk persists throughout the entire lifetime of the falling fingers (MacMinn & Juanes Reference MacMinn and Juanes2013). In this work, because dissolution structures are directly linked to the spatial arrangement of fingers, it is not surprising to find that the dissolution patterns also coarsen with depth (figure 7), in much the same way that fingers coarsen in the bulk. Here, we quantify this coarsening by counting the number dissolved polygonal rings ( $N_{polygon}$ ) as a function of depth ( $z$ ), as seen in figure 7 (black and white frames) at $t=14$ (figure 8). The analysis shows that the dissolution patterns coarsen with depth, and exhibit a robust power-law scaling:

(5.1) $$\begin{eqnarray}N_{polygon}\sim z^{-1}.\end{eqnarray}$$

Figure 7. Horizontal slices of the reaction rate $r={\it\phi}F({\it\alpha})\boldsymbol{{\rm\nabla}}{\it\alpha}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\alpha}$ (colour frames) and the relative porosity increase ${\rm\Delta}{\it\phi}/{\it\phi}_{0}$ (black and white frames) at different depths: (a) $z/\text{H}=0.01$ ; (b) $z/\text{H}=0.14$ ; (c) $z/\text{H}=0.27$ ; (d) $z/\text{H}=0.54$ ; (e) $z/\text{H}=0.82$ ; and ( f) $z/\text{H}=1.0$ . The colour frames for $r$ are at different times: (a) $t/\text{H}=6.0$ ; (b) $t/\text{H}=1.2$ ;(c) $t/\text{H}=2.7$ ; (d) $t/\text{H}=5.0$ ; (e) $t/\text{H}=7.8$ ; and ( f) $t/\text{H}=10.1$ (subtracted from a background image corresponding to $t/\text{H}=8.9$ for clarity). Red indicates high value (colourmap not shown). The black and white frames for ${\rm\Delta}{\it\phi}/{\it\phi}_{0}$ are all at time $t/\text{H}=14.0$ . The colourmap range varies with depth to best reflect the dissolution structure.

Figure 8. Number of dissolved polygonal rings as a function of depth at $t=14$ .

The above scaling was suggested by an earlier study in two dimensions (MacMinn & Juanes Reference MacMinn and Juanes2013), where the number of fingers in the bulk is found to scale with convective time as: $N_{fingers}\sim t_{conv}^{-1}$ . This is indeed consistent with (5.1) if we consider a constant finger falling speed so that $z\sim t_{conv}$ . However, understanding the mechanism that leads to such scaling remains a future task in this study.

6. Summary

In this paper, we study rock dissolution as a result of geochemical reactions during convective mixing in porous media, in the context of $\text{CO}_{2}$ sequestration. To couple geochemistry with flow and transport, we assume a mixing-limited reactive transport system where reactions reach chemical equilibrium instantaneously compared to the transport time scale. This allows us to adopt a formulation that completely decouples transport from reaction, and describes the local reaction rate as a function of scalar dissipation and chemical speciation (De Simoni et al. Reference De Simoni, Carrera, Sánchez-Vila and Guadagnini2005). Using high-resolution simulations, we investigate the interplay between flow and reaction as a result of local permeability changes from rock dissolution. Because the reactive system we study here is mixing-limited, we find that the rate of rock dissolution is initially high in regions of high fluid mixing, which leads to dissolution patterns that closely follow the structure of the mixing field. However, geochemical reactions shut down much earlier than convective mixing shuts down; a result of the highly nonlinear behaviour of chemical speciation. This feature of a speciation-controlled shutdown highlight the important role that the details of the geochemical equilibrium play in this hydrodynamics–reaction coupled process.

In both 2D and 3D simulations, we find that rock dissolution focuses on the top boundary, leaving the rock in the rest of the domain almost undissolved. As a result of the porosity increase at the top, we observe a significant increase in the rate at which the denser fluid enters the domain. This increase in solubility rate could enhance the effectiveness of $\text{CO}_{2}$ trapping in the context of a migrating plume (MacMinn et al. Reference MacMinn, Szulczewski and Juanes2011). In the bulk, we see weak coupling from the rock dissolution patterns to the fingering dynamics, since the permeability increase is small below the boundary layer. We investigate the dissolution morphology in three dimensions and observe that the pattern at each depth exhibits a polygonal tessellation structure. This structure corresponds closely to the spatial arrangement of columnar fingers, and it also coarsens with depth in a similar way that convective fingers coarsen with time in the bulk of the domain.