Dynamics of CO2 Density‐Driven Flow in Carbonate Aquifers: Effects of Dispersion and Geochemistry

The dissolution of carbon‐dioxide (CO2) in deep saline aquifers is an important trapping mechanism in carbon storage. This process is triggered by unstable high‐density CO2 front, which later promotes density‐driven mixing, hydrodynamic dispersion of CO2, and favors the long‐term sequestration. In many former studies, effects of hydrodynamic dispersion and multispecies geochemical reactions have been ignored. This work elaborates the impacts of these simplifications on the dynamics of convective mixing by numerical simulations. Geochemical effects were studied by the implementation of rock‐fluid and fluid‐fluid interactions for a typical carbonate aquifer. Results show that accounting for the hydrodynamic dispersion decreases the convection onset time and increases the CO2 dissolution flux, which is more significant in larger dispersivities and Rayleigh numbers. Results indicate that carbonate geochemical reactions intense the long‐term overall efficiency of the process, while decrease the total amount of sequestered carbon in the diffusion‐dominated period. Results also reinforce the importance of realistic geochemical representation and importance of spatial and temporal dependence of the reactions pathway, subsequent to the finger development for more detailed simulation of the CO2 storage process.


Introduction
Net-zero carbon emission is a long-term strategy to save the planet from eminent danger of catastrophic "global warming" (Rogelj et al., 2015). Carbon capture and storage (CCS) is one of the key technologies to be used to target the long-term goal of decarbonization and net-zero carbon emission in the energy sector (Pye et al., 2017). Among available options for carbon storage (CS), deep saline aquifers have the highest capacity (Celia, 2017;Godec et al., 2011). In this scheme, CO 2 is separated from emission sources before being released to the atmosphere, and then injected into saline aquifers, where it is supposed to be trapped for a long time. CO 2 is usually injected in supercritical or even liquid state, which occurs in aquifers located roughly at depths of 800 m or deeper (Bachu, 2003;Bachu & Adams, 2003). In such a scenario, possible trapping mechanisms include structural trapping, residual trapping, solubility trapping, and trapping through mineralization (Bakhshi et al., 2018;De Silva et al., 2015;Emami-Meybodi et al., 2015;Soltanian et al., 2019).
In a typical carbon storage injection scenario, supercritical CO 2 density is roughly 260-760 kg/m 3 with the viscosity ranging between 0.02 and 0.06 cP. So, due to the smaller density and viscosity when the supercritical CO 2 is injected into a deep aquifer, it will eventually rise upward (i.e., gravity segregation) until it reaches the impermeable caprock (Ide et al., 2007;Juanes et al., 2008). Then, CO 2 will dissolve into the aqueous phase, resulting in a denser (roughly 1%) phase being on top of the native aquifer fluid. This instability promotes CO 2 mixing due to convection, which is beneficial for the whole process. Beside CS, convective mixing is relevant to a wide range of other natural and industrial applications, such as enhanced oil recovery (EOR) (Kahrobaei et al., 2012;Sabet et al., 2018;Sheng, 2011), thermal and compositional mixing in the sea and oceanic water and karstification (Anderson et al., 1979;Ronen et al., 1988;Smith, 2004;Class et al., 2020), industrial chemical production (Wylock et al., 2014(Wylock et al., , 2017 and enhanced heat transfer using nanofluids (Heris et al., 2007;Li et al., 2003).
Abstract The dissolution of carbon-dioxide (CO 2 ) in deep saline aquifers is an important trapping mechanism in carbon storage. This process is triggered by unstable high-density CO 2 front, which later promotes density-driven mixing, hydrodynamic dispersion of CO 2 , and favors the long-term sequestration. In many former studies, effects of hydrodynamic dispersion and multispecies geochemical reactions have been ignored. This work elaborates the impacts of these simplifications on the dynamics of convective mixing by numerical simulations. Geochemical effects were studied by the implementation of rockfluid and fluid-fluid interactions for a typical carbonate aquifer. Results show that accounting for the hydrodynamic dispersion decreases the convection onset time and increases the CO 2 dissolution flux, which is more significant in larger dispersivities and Rayleigh numbers. Results indicate that carbonate geochemical reactions intense the long-term overall efficiency of the process, while decrease the total amount of sequestered carbon in the diffusion-dominated period. Results also reinforce the importance of realistic geochemical representation and importance of spatial and temporal dependence of the reactions pathway, subsequent to the finger development for more detailed simulation of the CO 2 storage process.
The convective mixing of CO 2 is usually simulated by convection and diffusion solute transport mechanisms. In natural porous media, the local fluid velocities distribute around the mean velocity due to rock heterogeneity. The effect of such a phenomenon on solute transport is modeled by hydrodynamic dispersion (Dentz et al., 2017), which is usually neglected in modeling the convective mixing process (Emami-Meybodi et al., 2015). Using numerical simulations Hidalgo and Carrera (2009) showed that dispersion can reduce the onset of instability by two order of magnitudes. Ghesmat et al. (2011a) showed that earlier convection onset and higher dissolution flux are expected as the medium dispersivity is increased. The effects of the transversal to longitudinal dispersivity ratios were found to be minimal on the dissolution flux while having an impact on the fingering pattern. Wen et al. (2018) showed that the increase of mechanical dispersion coarsens the fingering pattern by increasing the plume spacing, which was shown experimentally by Wang et al. (2016) and Liang et al. (2018) as well. Recently, Michel-Meyer et al. (2020 experimentally studied the effect of aquifer background flow on CO 2 dissolution and found that the background flow suppresses the development of fingers, while the dissolution rate is not sensitive to the background flow. Rock-fluid and fluid-fluid interactions are inseparable parts of CS, and can increase or decrease the efficiency of the process under different circumstances (Audigane et al., 2007;Cardoso & Andres, 2014;Liu et al., 2011;Sainz-Garcia et al., 2017). Although the reactive CO 2 convective mixing has been widely studied in the literature, the effects of complex geochemical interactions on the dynamics of the process from the early to the late (convection shutdown, which takes place as the domain starts saturating with CO 2 ) stage have not been studied. The conducted research on the effects of geochemical interactions on CO 2 convective mixing can be categorized into three main groups, dealing with sandstone, carbonate and generic (in the form a sink/source term or generic A + B → C reaction) interactions.

Generic Studies
In these studies, the reactions are usually coupled with flow equations by introduction of a sink/source term in the form of Damköhler (Da) number, defined as the ratio of the diffusion time scale to the reaction time scale. Lack of interpretation of the physical-chemical processes that are relevant for geological carbon storage, ignoring different paths of the geochemical reactions (by reaction path we aim determining if the reaction is proceeding in the forward or backward direction, i.e., dissolution/precipitation in the case of rock-fluid interaction) in time and space as well as complex multimineral interactions can be mentioned as limitations of these studies. However, they provide a clear insight into the effects of reaction on the process by simplifying the chemistry part. Andres and Cardoso (2011) studied the effect of a first-order reaction on the onset of density-driven convection. They concluded that the stability of the system can be determined by the dimensionless group of Da/ 2  , where  stands for Rayleigh number, which is a measure of the natural convection in porous media as a ratio of the driving buoyancy forces to the diffusive forces. They showed that above the critical value of Da/    2 3 2 10  the geochemical reactions stabilize the diffusive boundary such that no convective mixing is possible, and all dissolution occurs by diffusion. A similar conclusion was drawn by Kim and Choi (2014). Loodts et al. (2017) and Jotkar et al. (2019) studied the effect of a simple A + B → C reaction on the dynamics of the process. They showed that the effects of the reaction is highly dependent on the Rayleigh ( )  number, as well as the reaction direction (forward and/backward paths). Ghoshal et al. (2017) showed that in such systems, the impact of reaction on the process is highly determined by the ratio of density change coefficients for the product and the reactant (β C /β A ) in the aqueous phase. They showed that as β C /β A increases, the large density contribution from the product elongates fingers, resulting in a deeper chemically active ERFANI ET AL. 10.1029/2020WR027829 2 of 26 domain. As the implemented reaction is simplistic and cannot mimic the complex nature of the geochemical interactions of CO 2 -rock system, it is difficult to upscale the obtained results and provide field-scale implications and insights for CS. The same discussion applies to studies by Ward et al. (2014) and Hill and Morad (2014). They just considered the consumption of the solute in the simulation domain by neglecting complex multimineral interactions, which can result in spatial and temporal nonuniform solute consumption/production. To this end, Lei and Luo (2019) concluded from reactive lattice Boltzmann simulations that the reaction path determines if the chemical reaction enhances or damps the density instabilities. Babaei and Islam (2018) simulated reactive convective dissolution in a double-porosity domain with mass transfer between mobile and immobile regions. By simplifying geochemical interactions into a single second-order reaction, they showed the importance of the interplay between reaction rate and mass transfer coefficient intensities as the determining factor for the effect of geochemistry on the CO 2 storage. They concluded that due to a lower intensity of the geochemical reactions in the mobile region, intermediate values of mass transfer rate serve as a threshold below which the mass transfer coefficient does not affect the overall CO 2 storage; however, for high mass transfer intensities, overall CO 2 storage increases.

Sandstones
Sandstones are usually represented by silicate-rich (such as quartz and feldspar minerals) mineralogies, which are mostly mixed with clay (such as kaolinite and mixed-layer minerals) and some carbonates (calcite and/dolomite), considered as cementing agents (Tallman, 1949). Due to relatively slow reactions, the kinetic approach is usually implemented to study the effects of geochemistry on CO 2 convective mixing in sandstones (Erfani et al., 2020). Ennis-King and Paterson (2007) concluded that the density increase due to ion concentrations of Ca 2+ , Mg 2+ can be of the same order of magnitude as the density increase due to dissolved CO 2 . So, due to the strong impact of density changes on the CO 2 convective mixing, it is crucial to account for the effect of all species on the density. Zhang et al. (2009), also concluded that the geochemical interactions increase the amount of dissolved CO 2 in the aquifer. Ghesmat et al. (2011b) lumped all possible rock-fluid and fluid-fluid interactions into a single term and considered scenarios with carbon precipitation reaction, mostly relevant to sandstone aquifers. They showed that the geochemistry usually stabilizes the diffusive boundary. However, different reaction paths at different time scales and different locations of the aquifer (subsequent to plume propagation) were neglected. Such considerations were later shown to be crucial in sandstone aquifers by Erfani et al. (2020). Cardoso and Andres (2014) experimentally studied the effect of an irreversible precipitation reaction on convective mixing. They concluded that such form of interactions can curtail the solute spreading and convective onset in the system. Erfani et al. (2020) showed the scale-dependency and time-dependency of geochemical interactions in a typical sandstone aquifer. They showed that in the field-scale, the geochemical interactions decrease the solubility trapping, while increasing the total amount of sequestered carbon through mineralization (same conclusion as Ghesmat et al., 2011b andAndres, 2014). However, at early time the rock dissolution was shown to be the dominant mechanism, which stimulated the instabilities (contradictory to Ghesmat et al. (2011b)). They reinforced the importance of considering the time and spatial dependence of the reaction path (dissolution/precipitation). Moreover, important precipitating carbonate minerals in sandstone aquifers were shown to be dawsonite, calcite/dolomite, ankerite, and magnesite in different circumstances (Erfani et al., 2020;Moore et al., 2005;White et al., 2005;Xu et al., 2005).

Carbonates
Calcite and dolomite are the most abundant minerals found in carbonate rocks. Gypsum is also found in sulfur-rich aquifers. The reaction rates of these minerals in circumstances relevant to carbon storage are 3-4 orders of magnitude higher than major sandstone minerals (Palandri & Kharaka, 2004). Hence, the local equilibrium approach (LEA) is widely used for coupling reactions and flow (especially for convective mixing which is relatively slow) in carbonate aquifers (Emami-Meybodi et al., 2015). Fu et al. (2015) and Hidalgo et al. (2015) studied the reactive density-driven mixing in carbonate aquifers through 2D/3D simulations and concluded that the rock dissolution mostly takes place at the upper boundary, they have also shown the importance of chemical speciation on the results. They assumed that all species share the same ERFANI ET AL.

10.1029/2020WR027829
3 of 26 diffusion coefficient and simplified the geochemical interactions by introducing a mixing ratio between boundary and resident fluids (this approach was first introduced by De Simoni et al. (2007) for equilibrium reactions). In these studies, PHREEQC was used to get the chemical speciation at different mixing ratios between the boundary and resident fluids. Sainz-Garcia et al. (2017) showed that the chemical interactions enhance the overall process in carbonate aquifers through rock dissolution at the early stages of the process. They also emphasized the importance of noncarbonate mineral (i.e., gypsum) interactions on the results. Islam et al. (2016) coupled geochemical interactions of carbonate rock with convection-dissolution transport and found that the porosity and permeability changes (0.06%) are negligible and do not affect the hydrodynamics at all (a similar conclusion was drawn by Erfani et al. (2020), Fu et al. (2015), and Sainz-Garcia et al. (2017)).
From the above considerations, it is clear that a specific study accounting for the effect of complex representative geochemical interactions on the entire CO 2 convective mixing dynamics is still missing, and this study tries to fill this gap. For this purpose, a complementary set of geochemical (rock-fluid and fluid-fluid) reactions relevant to carbonate aquifers (including calcite, dolomite and gypsum minerals) is coupled with high-resolution 2D/3D multicomponent convection-diffusion numerical solvers. The effects of different ions concentrations on density are considered, which delineate the effect of precipitation and dissolution reactions on the density. Then, extended simulations are performed over a wide range of -numbers to study the effects of geochemical interactions on the dynamics of convective mixing, relevant to CS in carbonate aquifers. Moreover, former studies mostly assume diffusion and convection transport of specie. However, we have demonstrated the impact of hydrodynamic dispersion on the process. Lastly, we study the convective mixing process in a heterogeneous aquifer to demonstrate the impact of geochemistry and dispersion in nonidealized circumstances.
The manuscript is organized as follows: in Section 2 the model setup, as well as the formulation and postprocessing schemes, will be presented. In Section 3, the model properties and numerical methods will be provided. The numerical results and discussion will be provided in Section 4, then the paper will be closed with concluding remarks.

Problem Description and Formulation
Two-dimensional and three-dimensional rectangular/cubic domains were considered to represent a confined isothermal aquifer, saturated with brine with length and thickness of L (m) and height H (m), respectively. Initially, the aquifer brine was assumed to be at equilibrium with the rock constituents (primary minerals) at the aquifer temperature and pressure. Such initial condition is based on the assumption that deep saline aquifers were intact for a long time (Hellevang et al., 2013;Sainz-Garcia et al., 2017). All domain boundaries were assumed to be closed (no flow boundary) except the top boundary, which was considered to be at a constant concentration, a schematic of the simulation domain is presented in Erfani et al. (2020). Additional assumptions are as follow: • The top CO 2 -brine interface is assumed to be sharp and remain at a constant concentration (constant concentration boundary condition) • The top boundary is assumed as an infinite source of CO 2 , so the pressure change due to dissolution is neglected (Ennis-King & Paterson, 2005; Islam et al., 2013) • Gaseous CO 2 has no interaction with the rock (Balashov et al., 2013;Gaus, 2010) • Homogeneous and isotropic permeability K (m 2 ) and porosity (ϕ) fields are considered. The aquifer is assumed to have the aspect ratio of unity with height (H) of 100 m • Initial rock composition/mineralogy is assumed to be calcite, dolomite, and gypsum • The Oberbeck-Boussinesq approximation (Chandrasekhar, 2013;Spiegel & Veronis, 1960;Valori et al., 2017) and Darcy's law are assumed valid. Convective flow does not generate large flow velocities and the flow is expected to be laminar (Landman & Schotting, 2007;Sheremet et al., 2015) The governing equations of flow and species transport are: is the gravitational acceleration, and  z n is the unit vector pointing upward.
where ϕ is porosity, c (mol 1 kg w ) is the species concentration and i  denotes the sink/source term of the species i due to the geochemical interactions. D i can represent diffusion or dispersion coefficient, depending on the assumptions made. The use of diffusion coefficient in simulation of convective mixing is conventional, while some researchers (Chevalier et al., 2015;Ghesmat et al., 2011a;Hidalgo & Carrera, 2009;Xie et al., 2011;Wen et al., 2018) tried to relax this assumption. In Equation 2, D (m 2 s −1 ) is the hydrodynamic dispersion tensor and is defined as (Bear & Bachmat, 2012): where α L (m) and α T (m) denote longitudinal and transversal dispersivities, respectively. δ ij is Kronecker delta,  operator is the Euclidean norm of a vector, is the molecular diffusion coefficient and I is the identity matrix.
(c) Density change due to species concentration is represented by where β i is defined as the coefficient of density change for species i as

Dimensionless Form of the Governing Equations
To take into account the aqueous phase density change and simplifying the governing equations, we lump different species with the same component into pseudocomponents (i.e.,  ).
Then, we can write Equation 4 in a differential form for each pseudocomponent and take the carbon pseudocomponent coefficient of density changes as the reference value (β r ).
where n ps denotes the total number of pseudocomponents. The -number is also defined based on the diffusion coefficient of the carbon pseudocomponent (i.e., reference diffusion coefficient ) and the density difference at the top boundary as: ERFANI ET AL.
10.1029/2020WR027829 5 of 26 Horne (1979) has shown that introducing a vector potential of the form u = ∇ × ψ into the formulation results in faster and more stable numerical solutions for the problem. It has also been shown (Aziz & Hellums, 1967;Hewitt et al., 2014;Hirasaki & Hellums, 1968) that the velocity and stream function are solenoidal (incompressible flow), so ∇ ⋅ ψ = 0. To derive the dimensionless form of the equations for a domain with an aspect ratio of unity, dimensionless variables are defined as: where * superscript denotes dimensionless variable,   c is the carbon pseudocomponent coefficient of density change, c r is the reference and c 0 is the initial concentration at the upper boundary and the aquifer, respectively. Then, one can derive the following equations.
Moreover, the reactive species transport equation for different species reads as: .

Geochemistry
In this work, geochemistry was modeled by coupling PhreeqcRM (Charlton & Parkhurst, 2011;Parkhurst & Wissmeier, 2015) with an in-house stream function based convective-dispersive transport simulator. Geochemical reactions were coupled with the flow equations by sequential noniterative approach (SNIA) (Carrayrou et al., 2004;Charlton & Parkhurst, 2011;Erfani et al., 2019;Farajzadeh et al., 2012;Kanney et al., 2003;An et al., 2021). As the carbonate mineral reactions are fast and the convective flow is relatively slow, LEA was followed for all geochemical reactions including solid and aqueous phase reactions. The primary minerals considered were calcite, dolomite, and gypsum, which are typical minerals in carbonate aquifers. Moreover, a complementary set of aqueous phase interactions between different species were considered for a more nature representative modeling of rock-fluid and fluid-fluid interactions (Erfani et al., 2019;Fu et al., 2013;Sainz-Garcia et al., 2017). Table 1 summarizes all considered solid and aqueous phase reactions. Notably, the reaction equilibrium constants were taken from the literature and temperature dependency was modeled by the analytical five-term formulation or Arrhenius equation, in case, the analytical coefficients were not available (Blanc et al., 2012;Laidler, 1984;Parkhurst & Appelo, 1999). Based on the provided reactions, 20 species, as well as six pseudocomponents and 13 reactions, were defined in the model.
where τ t and τ r are transport and reactive residence times,  denotes the average transport velocity, V p is the pore volume, c eq is CO 2 equilibrium concentration, r i is ith mineral reaction rate and i  is the corresponding reactive surface of the mineral. Damköhler number compares the relative importance of transport and reaction. Even for the slowest reaction in the system (i.e., dolomite, Palandri & Kharaka, 2004) and for a relatively high -number in convection-dominated region ( , which was obtained by the displacement velocity of the concentration front in convection-dominated regime) the Damköhler number will be much higher than 1 (Da I ≫ 1), which shows the validity of LEA assumption (Ennis-King & Paterson, 2007;Erfani et al., 2019;Fu et al., 2015;Nasralla et al., 2015).

Convective Dissolution Dynamic Measures
To study the dynamic behavior of convective transport and characterizing the impact of geochemical interactions on the process, several quantitative measures are defined. In each time step, the amount of stored carbon due to mineral trapping from the start of the process, in mole is calculated as:  (Blanc et al., 2012;Parkhurst & Appelo, 2013) total amount of precipitated/dissolved mineral at each time step, we can update the porosity distribution in the domain as: where V m,n (m 3 mol −1 ) denotes the molar volume of the nth mineral and n m is the total number of minerals (i.e., calcite, dolomite, and gypsum in this case) and  t i is the porosity of the ith grid block at time t. Moreover, the amount of sequestered carbon due to solubility trapping mechanism is calculated by multiplying aquifer pore volume by average carbon pseudocomponent concentration      c C : The total amount of stored carbon can be obtained by summation of solubility-based and mineral-based carbon storage   : where A (m 2 ) denotes the domain area normal to dominant flow direction (i.e., z-direction). In the diffusion-dominated regime for nonreactive case the dissolution flux decays as Paoli et al., 2017;Loodts et al., 2017). The next global characteristic measure is the dispersion width (σ z (t)), which describes the average width of a particular spatial distribution in the dominant direction of the flow (Amooie et al., 2018;Sudicky, 1986). The second central spatial moment (σ z ) is usually analyzed to quantify spreading or macrodispersion during solute transport.
where c denotes the concentration of the main component, which is defined based on the dimensionless carbon pseudocomponent concentration. z c is the vertical position of plume center of mass (i.e.,     / c z Cz C ) and  is used as the domain averaging operator.

Model Description
To simulate the convective-diffusive transport process, at each time step Equation 7 should be solved over the domain to get the stream functions   * * ( , ) x y by which we can calculate the velocity map * * * ( , , ) x y z u u u . Next, Equation 8 (excluding the sink/source term) should be solved for all species separately to get the dimensionless concentration map of all species * ( ) i c . Afterward, the concentrations are transferred into the dimensional domain and imported into PhreeqcRM module to be updated based on the geochemical reactions mentioned in Table 1.
We solved Equation 7 with the discrete Fourier transform (DFT) (Cooley et al., 1970;Schumann & Sweet, 1988) method while Equation 8 was solved by the forward-time central-space (FTCS) (Kurtz et al., 1978) scheme. FTCS is a computationally efficient scheme while it has instability issues at large time steps. To avoid this, we chose a relatively small time step, which improves the accuracy of flow-geochemistry coupling in SNIA (Carrayrou et al., 2004).
To study the effects of dispersion and geochemical interactions on the dynamics of density-driven mixing in carbonate aquifers, the domains under investigation were field-scale 2D/3D domains with dimensions of 100 m with the aspect ratio of unity. The dimensionless time step was set to 10 −6 , for low to moderate ERFANI ET AL.

10.1029/2020WR027829
8 of 26 -numbers. Smaller numbers were chosen for higher -numbers to guarantee the convergence of the numerical methods. In order to develop finger-like front propagation, the CO 2 -water interface concentration was perturbed by a random noise of 0.5% (Riaz et al., 2006). The coefficient of density change for different pseudocomponents were calculated by fitting a line on the plot of density versus component concentration, obtained from PHREEQC output (Table S2). Diffusion coefficients of different species for different temperatures were taken from the literature (Cussler, 2009;Lasaga, 2014) and extrapolated to the simulation temperature, which are mentioned in Table S3.
Initial water composition was taken from Sainz-Garcia et al. (2017) and the composition of equilibrated water with rock minerals is shown in Table S4. Simulation temperature and pressure were set to 80 °C and 100 atm, respectively. Initial domain porosity (ϕ) was set to 0.15 and the -number was controlled by changing the permeability, which varied between 6 and 400 mD.

Results and Discussion
In this section, we investigate the density-driven CO 2 mixing in carbonate aquifers. 2D and 3D simulation results are presented for reactive and nonreactive cases over extended simulation times. Moreover, the effect of domain dispersivity was studied on the dynamics of the process.

Convective Mixing in Nonreactive Cases Without Dispersion
For the nonreactive cases, the simulator produces consistent results and hydrodynamics, compared to the existing previous modeling codes of subsurface CO 2 convective mixing (Amooie et al., 2018;Cheng et al., 2012;Islam et al., 2013;Sainz-Garcia et al., 2017;Slim, 2014). The model was also validated quantitatively and qualitatively against Cheng et al. (2012) in Section S1 as well as Slim (2014) shown in Figure S3. Figure 1 shows the time evolution of carbon pseudocomponent concentration profile in 2D domain for two different -numbers (namely  3000  and  12000  ), excluding the effect of dispersion (α L = 0). The sensitivity of the simulation results on the model resolution in 2D was also studied in Figure S3 and we chose the 300 × 300 grid blocks for the simulations. As the process starts, CO 2 enters the domain by diffusion and it forms a diffuse layer on top of the simulation domain. This regime is called diffuse regime, in which the dissolution flux ( ) decays proportional to t −0.5 , where t is time Slim (2014). After some times, sufficient dense fluid accumulates beneath the upper boundary, which amplifies the perturbations and triggers the convective instabilities. This is called linear-growth regime. This regime ends with the onset of convective mixing. Estimation of the onset time has been the focus of many studies in the literature (Amooie et al., 2018;Emami-Meybodi et al., 2015;Javaheri et al., 2009;Rees et al., 2008).
After the convective mixing onset, the dissolution flux deviates from t −0.5 and starts growing into a local maximum, so this is called flux-growth regime. The convective fingers become macroscopic and they mostly move downward. Once the fingers are sufficiently long, they start lateral movements and interact with each other. This is called merging regime, because the fingers start merging from the roots (i.e., root zipping mechanism) as a result of lateral movements. Root zipping is a phenomenon which affects the macroscopic fingers and results in merging from the top part (i.e., close to top boundary). The merged fingers form megaplumes (which are bigger secondary plumes) and move downward. Formation of megaplumes make lateral gravity counter-current around them, that hinders the downward movement of the small neighbor plumes (i.e., wake mechanism). As a result, later plumes are hindered from downward movement and eventually merge with existing megaplumes by drafting (i.e., the lateral movement of the plumes due to existing upward velocity around bigger adjacent plumes). During this period, the dissolution flux stabilizes, which is recognized as constant-flux regime (De Paoli et al., 2017;Green & Ennis-King, 2018;Slim, 2014).
Finally, the fingers reach the lower boundary of the domain and the flow enters its final dynamic stage, usually called shutdown regime. After the fingers hit the bottom wall, the CO 2 rich fluid starts to move upwards with the returning flow. Once this dense fluid reaches the top wall, the driving force for convection is ERFANI ET AL.  decreased, the flux declines rapidly, and eventually the convection is shut down. In this period, the domain starts saturating with dense fluid, which results in a drastic decline of the dissolution flux. Figure 2 shows the development of the carbon pseudocomponent concentration in the 3D domain with the permeability of 100 mD, in which the same phenomena can be observed. Figure 3 shows the carbon pseudocomponent concentration distribution in a horizontal plane close to the upper boundary of the 3D simulation domain at different times. The snapshots demonstrate the coarsening pattern during the convection process. In early times ( Figure 3a) the islands can be seen near the top boundary due to amplification of the imposed perturbations in the boundary concentration. Figure 3b shows the configuration of maze-like structures, which later convert to cellular network (Figure 3c). The evolution of such patterns and coarsening dynamics have been discussed earlier by Fu et al. (2013)   The described flow regimes can be elaborated by the development of dissolution flux ( ) versus time. Figure 4a shows the dissolution flux for four different cases (   12,000, 6,000, 750, and 375) in the 2D domain versus time. The flow regimes discussed earlier are marked on the main panel based on  6,000  case, which corresponds to the domain permeability of 100 mD. Side panels (Figures 4b  and 4c) show the time behavior of the concentration field in a horizontal plane close to the top boundary (z D = 0.1). The development of initial fingers and their interactions are obvious in these subfigures. An interesting point is that, after the formation of megaplumes, the protoplumes (small plumes generated in middle-to-late-time intervals inside or next to megaplumes) generated later will join them as it has been discussed earlier. Figure 5 shows the average concentration of the carbon pseudocomponent in the domain as well as the carbon dissolution rate for the 2D and 3D simulation domains for different -numbers. Noticeably, the 3D simulations show a higher dissolution rate and correspondingly higher average concentration in the domain, which was also reported by other researchers (Amooie et al., 2018;Hewitt et al., 2014;Pau et al., 2010) for comparable simulations. As it is shown, the 2D simulations show more fluctuations in the dissolution ERFANI ET AL.  flux, compared to 3D cases, which is expected due to averaging over a smaller domain. The density-driven mixing characteristics such as convection onset and maximum dissolution flux are fairly preserved and ≈ 16% difference between the stabilized dissolution flux can be observed between 3D and 2D simulations for a domain of 100 m in the case of  6,000  (permeability of 100 mD). Pau et al. (2010) reported 25% difference in the stabilized mass fluxes between 3D and 2D simulations for a domain of 32 m height. Moreover, Hewitt et al. (2014) reported 40% difference between 3D and 2D simulations. Notably, both have seen higher dissolution flux in the case of 3D compared to 2D simulations. ERFANI ET AL.

Effects of Dispersion on CO 2 Convective Mixing
In this section, we investigate the effects of dispersion on density-driven mixing at different -numbers. Figure 6 shows the evolution of carbon dissolution flux for different -numbers ( = 1,500, 6,000, and 12,000, which correspond to different domain permeabilities) and different values of longitudinal dispersivity (α L ). In all simulations, the value of transversal dispersivity was considered as α T = 0.1α L (Gelhar et al., 1992;Hidalgo & Carrera, 2009) and the model was run by the resolution of 300×300 grid blocks. This dispersivity is equivalent to the macroscopic dispersion in Liang et al. (2018). As it is shown, the effect of dispersion is considerably more significant in higher -numbers. Hydrodynamic dispersion (instead of only diffusion) results in earlier onset time because the dependence of dispersion on local velocities results in more nonuniformity of the diffusive front in early times. It also amplifies early instabilities due to velocity distribution in the system and hence, an earlier onset of convection.
For the simulation cases with α L /H = 10 −3 , there is no significant difference compared to diffusive condition (α L = 0), while the difference is more visible in the case of α L /H = 10 −2 , in which the dispersion effect is more pronounced in higher -numbers because of larger velocities indeed. The dynamics of the process considerably change for the case of α L /H = 10 −1 , (longitudinal dispersivity might seem artificially high in this case but this will amplify the effect and hence makes it easier to conclude the effect of dispersion) and especially in Rayleigh numbers of 6,000 and 12,000.
As seen in Figure 6, the onset time is lower and the dissolution flux stabilizes on a higher value compared to the diffusive case due to dispersion. Moreover, in the diffusive regime, the dissolution flux is independent of longitudinal dispersivity as well as -number because the velocities in the domain are still not significant to affect the dissolution rate and dispersion, so the system remains diffusion controlled (Note: in Figure 3 of Hidalgo and Carrera (2009)), the dissolution rates show a significant difference in the diffusive regime for different values of α L and interestingly the higher the longitudinal dispersivity the lower the dissolution rate. It can be due to different definition of dimensionless time and dependence of the dimensionless time to the dispersivity. Also, the dimensionless time scale is different for the case with α L = 0 in Figures 3 and Figure 2 which is unexpected.).
Carbon concentration distribution shows a different fingering pattern for the simulation cases with α L /H = 10 −1 compared to other cases. In these cases, fewer plumes are generated and the dissolution flux from the finger stem is higher as a result of higher local velocity in these parts. The formed fingers are also thicker, more laterally spread with larger spacing between them compared to cases with a lower α L . ERFANI ET AL.   Noticeably, the diffuse layer on top of the domain is also thicker compared to other cases. Similar results were reported in Liang et al. (2018).

Effects of Geochemistry on Solutal CO 2 Convection
Mixing of CO 2 -rich acidic fluid into the aquifer disrupts the equilibrium condition and triggers the geochemical interactions. The acidic water will dissolve calcite and dolomite from the aquifer rock, which increases the concentration of calcium and carbon in the swept area. Such interactions increase the local density difference while decreasing the local CO 2 concentration difference, so in diffusion-dominated regime the carbon dissolution flux ( ) is lower for the reactive case. In contrast, the geochemical interactions promote enhanced convective mixing and result in earlier onset and higher flux in the convection-dominated regime. As the calcium concentration increases due to the rock dissolution, gypsum precipitation is stimulated, which works as a calcium sink, resulting in more calcite and dolomite dissolution. Figure 7 shows the concentration distribution of different components and pH in convection-dominated regime at a snapshot. the flux in the diffusion-dominated regime is lower for the reactive cases corresponding to the mechanisms discussed. Figure 8b depicts the carbon dissolution flux from the aquifer rock for different simulation cases. At early times the dissolution flux is controlled by diffusion of CO 2 from the boundary into the domain. As the diffusion starts the water acidity increases and the generated protons in the aquifer are consumed by calcite and dolomite dissolution (hence the effect of geochemical interactions in this regime is not significant). After the convection onset, the dissolution rate considerably increases which enhances the total flux ( t  ) and mixing in the aquifer as well. After this period, the CO 2 and rock dissolution flux balance and rock dissolution rates show a constant rate regime. Afterward, the domain starts saturating with CO 2 which decreases the rock dissolution as well as the total flux significantly (i.e., shutdown regime).  Figure 9. Distribution of change in (a) calcite, (b) dolomite, and (c) gypsum content for the simulation case with  = 6,000 at t D = 10 −3 (corresponding to ≈100 years for a domain of 100 m). The changes in mineral content are illustrated in (log mol 1 kg w ). Notably, calcite and dolomite minerals dissolve from the aquifer rock into the aqueous phase, while gypsum precipitates due to CO 2 storage in the aquifer. are the dominant mechanisms. Such rock-fluid interactions enhance the CO 2 sequestration flux and spreading in the aquifer. At this snapshot, the calcite and dolomite dissolution are dominant on the edge of the plumes. Since at these regions, the CO 2 concentration, as well as pH, is low, while in the center of the plumes the CO 2 concentration is close to the saturation state, which hinders the rock dissolution and local precipitation takes place (see Figure S4). The opposite takes place for gypsum. These results reinforce the importance of considering different paths for different reactions subsequent to plume propagation in the domain for a complex representation of the geochemical effects, representative of the nature of the process (Islam et al., 2016;Sainz-Garcia et al., 2017). Figure S5 shows the distribution of porosity improvement in the domain for the simulation case of  = 6,000 at t ≈ 100 years (corresponding to t D = 10 −3 ). As shown, the permeability increase happens mostly at the top boundary of the domain and the porosity increase is not significant to change the hydrodynamics. Notably, in the current geochemical setup, the porosity increase is moderated by the gypsum precipitation. Also, as the molar volume of gypsum (74.69 cc/mol) is almost double that of calcite (36.93 cc/mol), the gypsum precipitation has a higher effect on the porosity change compared to calcite dissolution. Figure 10 represents the effects of geochemical interactions and dispersion on the global dynamic measures of convective mixing. Figure 10a shows that both geochemical interactions and hydrodynamic dispersion increase the average concentration of the carbon pseudocomponent in the domain (i.e., enhanced flux), especially in the convection-dominant regime. The effect of geochemical interactions on total stored carbon ( tot  ) is nonuniform and different for different -numbers. Before the convection onset (t c ), rock-fluid interactions decrease the total amount of stored carbon through decreasing the concentration difference, which is the main driving force in the diffusion-dominant regime. However, after the onset time, the total amount of stored carbon is higher for the reactive cases as the mixing and CO 2 flux are higher in these cases (see Figures 8 and 10b).
The geochemical interactions increase the mixing length (z m ) in the simulation domain by rock dissolution and enhanced mixing at finger-tips through dissolution/precipitation (see Figure S4). Conversely, the dispersive simulations with α L /H = 10 −2 show smaller mixing length, while the average concentration of carbon is higher in the domain. Such a phenomenon is a result of smaller vertical velocity component at the finger-tips compared to the finger stems. The dispersion results in enhanced flux in the finger-stem areas compared to finger-tips especially for the small to the average longitudinal dispersivities (see Figure 6). In high dispersivities, the mixing length increases as the effect of dispersion as well as the onset time are considerably enhanced. Figure 10d illustrates the dispersion width versus dimensionless time for different simulation cases. In early times, spreading and mixing are driven by diffusion and the diffusive boundary expands, so σ z increases with a classical Fickian scaling (  0.5 z D t ) for all simulation cases, demonstrating the minimal impact of geochemical interactions and hydrodynamic dispersion in this period. In early convection-dominated regime (t > t c ) the dispersion width deviated from the Fickian behavior and at the late convection period, the linear growth of the global dispersion width was observed for all cases, with ballistic ∼t D scaling (Amooie et al., 2018). This shows that hydrodynamic dispersion and geochemical interactions in the numerical model do not change the convective mixing fundamentally, while weakens or strengthens the existing mechanisms.

Field-Scale Implications for Heterogeneous Aquifers
In this section, we study the convective mixing in a heterogeneous media in the field-scale. For this purpose, we have generated a correlated permeability map of 300×300 mesh with a minimum, maximum and average permeability of 30, 550, and 300 mD, respectively, and with a normal distribution. The permeability map and distribution are shown in Figure S6. Figure 11 presents the results obtained for the base case (nonreactive and nondispersive) as well as reactive and dispersive (α L = 1 m) simulation case. The permeability heterogeneity results in an earlier onset of convection while the discussed regimes in Section 4.1 are still visible but with a higher fluctuation due to domain permeability heterogeneity. The reactive simulation showed an earlier onset of convection and a higher domain average concentration in agreement with the mechanisms discussed.

10.1029/2020WR027829
18 of 26 The carbon concentration distributions in Figure 11 show that in the heterogeneous simulation domain the fingering pattern is mostly imposed by local high permeability regions. Moreover, the dissolution flux shows higher fluctuations compared to homogeneous simulations. The obtained results for the heterogeneous aquifer show that the phenomenological study of the effect of different parameters on the process as well as flow regimes in the homogeneous simulation domain can be observed in the heterogeneous domains as well.

Concluding Remarks
This study presents numerical simulations of multicomponent reactive density-driven CO 2 flow in deep saline aquifers coupled with hydrodynamic dispersion and a complex set of carbonate geochemical reactions. The geochemical reactions make the overall process more complicated due to interconnection between local geochemical reactions and hydrodynamics of the flow and transport. The provided results enhance ERFANI ET AL. our understanding of the reactive density-driven convection of CO 2 through coupling the multicomponent solute mixing with PhreeqcRM to include detailed geochemical modeling.
First, the dynamics of the process was reviewed in 2D and 3D domains as well as the coarsening patterns near the top boundary of the 3D simulations. It was shown that 2D and 3D simulations provide consistent ERFANI ET AL.  results with ≈16% difference between the stabilized dissolution flux of 3D and 2D simulations for a domain of 100 m in the case of  6,000  (permeability of 100 mD). In such cases, the 3D simulations showed a higher stabilized flux compared to 2D simulations which were in line with findings of other researchers.
The presence of hydrodynamic dispersion was shown to increase the amount of dissolved carbon in all mixing regimes. Such a phenomenon destabilizes the diffusive front and results in earlier dominance of convection in the process. The effects of dispersion were shown to be significantly different for different values of dimensionless dispersivity (α L /H). The effects of hydrodynamic dispersion were shown to be minimal for O(α L /H) ∼−3. Such effects strengthen and become important as the dispersion increased to O(α L /H) ∼−2. For dispersivities in the range of O(α L /H) ∼−1, the onset time and the value of stabilized dissolution flux were significantly different compared to convection-diffusion cases especially for high -numbers. The pattern and the shape of plume propagation, as well as the thickness of the diffusive layer on top of the domain, were also shown to be different for the case of O(α L /H) ∼−1, as presented in Section 4.2.
The rock dissolution was shown to be the dominant mechanism in reactive CO 2 convection in carbonate aquifers. While such a mechanism decreased the amount of stored carbon in the diffusion-dominant regime, it stimulated the convection onset and the amount of the total flux after the onset time. Such effects may not be captured unless by considering the effects of different ions on the density as well as taking advantage of realistic geochemical modeling. These effects were captured in this study by the modified stream function formulation, which takes into account the effect of all ions on the aqueous phase density. As it was shown, the carbonate (calcite and dolomite) precipitation and dissolution were dominant at the stem and the edge of the plumes, respectively. Notably, the opposite took place for gypsum. Moreover, the simulation results showed that the maximum porosity increase was in the order of 10 −2 after 100 years for a domain of 100 m, which mostly took place near the top boundary of the domain.
The presented simulation results for heterogeneous aquifers show that the permeability heterogeneity leads to an earlier convection onset time. Moreover, in this case, the pattern of heterogeneity imposes the fingering pattern and plume shapes. The fingers tend to spread mostly in the high permeability regions which result in a higher level of fluctuations in the dissolution flux. Notably, the phenomena observed in the homogeneous aquifers are still visible in the heterogeneous case which shows the applicability of the discussions provided.