Reactive Transport Modeling of the Enhancement of Density-Driven CO2 Convective Mixing in Carbonate Aquifers and its Potential Implication on Geological Carbon Sequestration

We study the convection and mixing of CO2 in a brine aquifer, where the spread of dissolved CO2 is enhanced because of geochemical reactions with the host formations (calcite and dolomite), in addition to the extensively studied, buoyancy-driven mixing. The nonlinear convection is investigated under the assumptions of instantaneous chemical equilibrium, and that the dissipation of carbonate rocks solely depends on flow and transport and chemical speciation depends only on the equilibrium thermodynamics of the chemical system. The extent of convection is quantified in term of the CO2 saturation volume of the storage formation. Our results suggest that the density increase of resident species causes significant enhancement in CO2 dissolution, although no significant porosity and permeability alterations are observed. Early saturation of the reservoir can have negative impact on CO2 sequestration.

Several studies including a few laboratory tests have been performed regarding convective mixing from perturbed saturated boundary [11][12][13][14][15][16][17][18][19] . However, the geochemical effect to this end is modeled to a lesser extent. Ennis-King and Paterson 20 investigated analytical and numerical methods and observed that the overall dissolution rate depends on the balance between the effects of permeability alterations and the consumption of dissolved CO 2 . They showed that the influence of ion concentrations (e.g., Ca +2 , Mg +2 ) on the fluid density alters the plume structure and favors faster dissolved plume development. Ghesmat et al. 21 used linear stability theory and direct numerical simulation to show that geochemical reactions stabilize the unstable diffusion boundary layer because of the consumption of dissolved CO 2 . Their results implied that more CO 2 can be trapped through mineral interactions. Zhang et al. classified dissolution-diffusion-convection process into four stages: (1) dissolution-dominated period; (2) diffusion-dominated period; (3) early convection-dominated period; and (4) late convection-dominated period 22 . Islam et al. 23 presented simulation results of convective mixing with reactions adding heterogeneity and geothermal effects. They concluded that at a fixed Damkohler number (Da), reaction orders make substantial difference of mixing over longer period of times. Geothermal gradient exhibited negligible impact. Ward et al. 24 studied reactive flow in the limit of high Ra in which the domain considered was deep, shallow or of intermediate depth, and for which the Da was respectively large, small or of order unity. For large Da the rapid reaction rate limits the plume depth and the boundary layer restricts the rate of solute transfer to the bulk volume, whereas for small Da the average solute transfer rate is ultimately limited by the domain depth and the convection is correspondingly weaker. All these previous researches were very generalized, in terms of Ra and Da.
Fu et al. 25 presented formation of rock-dissolution patterns that arose from a series of calcite dissolution reactions during convection. They used high-resolution simulations to examine the interplay between the density-driven hydrodynamic instability and the rock dissolution reactions and analyzed the impact on the macroscopic mass exchange rate. Their conclusion was the geochemical reactions terminate significantly earlier than the time when convective mixing stops. However, it should be clearly noted that geochemical effects may not always accelerate advection. A precipitation reaction such as that between the acidic brine and a rock formation rich in calcium feldspar promotes the deposition of solid calcite and kaolinite, removing CO 2 from the liquid phase. Such a reaction may actually attenuate convection motion 26 . Dai and co-authors 27 developed an integrated Monte Carlo method for simulating CO 2 and brine leakage from carbon sequestration and subsequent geochemical interactions in shallow aquifers. Their results showed shallow groundwater resources may degrade locally by reduced pH and increased total dissolved solids. Ennis-King and Paterson 20 and Fu et al. 25 showed results of porosity and permeability changes from reactions and subsequent impacts on convection mechanism. However, in this work we show, by performing reactive transport modeling with speciation of a series of calcite and dolomite dissolution reactions, that dissolution enhancement causes local   density increase and effects from porosity and permeability change are almost negligible even when the reservoir reaches near 100% dissolved CO 2 saturation.

Results and Discussions
We first have conducted sanity test of chemical speciation. CO 2 leak problem is solved considering chemical reactions of Eqs 3-8 by following the semianalytical approach in Yang et al. 28 . It is noteworthy that coupled reactive transport calculations strongly depend on accuracy of geochemical database. Different geochemical databases and uncertainty of thermodynamics data can severely affect the results. For the transport part analytical solution for the scenario with leaky well was used. Figure 1 shows that CO 2 stream in host carbonates after equilibrium dissolution can be increased by more than 1.5 folds. In our simulation the maximum equilibrium dissolved CO 2 concentration of 0.97 mol/l is used in saturated CO 2 -brine interface at the top. In order to quantify saturation area of the aquifer, average concentration formulation, , is used. The direct numerical simulations have been carried out for Ra numbers of 1000 and 10,000, where the associated mean permeability values are 10 and 100 mD, respectively. Figure 2 show time lapse of dissolved CO 2 fronts from 20 years period to the time the reservoir domain becomes completely saturated. In the case of no reactions, because Ra is low, convection is diffusion dominated, resulting in very slow plume advancement. Even after 900 years CO 2 propagation is still limited to upper 10% of the reservoir. However, reaction activities with calcite and dolomite make results completely different. By that time the entire porous formation becomes saturated with CO 2 . Whenever CO 2 arrives at any particular domain after equilibrium reactions, CO 2 stream in the form of − HCO 3 is dissolved from resident rocks, increasing local density of the brine. The increased density drives more instability which, in turn, causes plume boundaries to advance further downward. For this specific reservoir and thermophysical conditions in 900 years the pore volume is chock-full with dissolved CO 2 . The process adds significant feedback, however, in a negative sense because of early shutdown of both solubility and mineral trapping processes. Reactions constantly accelerate motion of the fronts. Because convective mixing shows diffusion dominance and the medium is homogeneous the cells do not form wormholes. Therefore, reaction fronts are planar throughout the process. No bifurcation occurs.   exhibit results of same Ra, however for heterogeneous permeability formation. Dykstra-Parson coefficient of 0.55 was applied in generating the distributions shown in Fig. 4. Initially, results do not vary much from its homogeneous counterpart other than the pattern of plume evolution. As convection proceeds in addition to heterogeneity effects concurrent geochemical reactions make noticeable difference in dissolution process. Because of local permeability variations and nonlinear flow dynamics from the beginning of plume development hopf-bifuraction of the cells occur and CO 2 stream added from carbonates help spread initially laterally and then finally dense phase sinks down. Thus, compared to no permeability contrasts, average concentration of CO 2 in the aquifer differs significantly with time. ~100% saturation is reached 50 years earlier. On the other hand, < 20% saturation only by reverse buoyant flow conveys clear message of reaction effects on CO 2 convection (see Fig. 5). We have also tested results of layered anisotropic heterogeneity = .
( ) Mean permeability is too small to add any effect resulting in almost same results as homogeneous reservoir.
Figures 6-8 plot concentration contours for the case of relatively high mean permeability 100 mD (Ra = 10,000). As expected, high permeability allows dissolving more CO 2 and consequently reaction effects become more pronounced. Among three cases, the most intense enhancement is observed in heterogeneous reservoir where CO 2 saturates completely in 500 years. Ra is high enough to initiate convection in early times and for heterogeneity the convection cells get more dynamicity because of availability of local high permeable favorable paths. Thus CO 2 dissolution swells, transporting more CO 2 saturated brine in contact with carbonates. The concurrent convection and reaction processes drive fast spreading of dense dissolved phase in the reservoir both laterally and downward. Figure 9 shows concentration profiles of Ra = 10,000. In this case layered heterogeneity adds little more reaction effect over the homogeneous distribution. For no reactions case maximum saturation obtained is ~40%. The early saturation due to geochemical effects has a negative effect on CCS in the sense that much injected CO 2 may remain as completely undissolved. As a result, potential of leakage through high permeable zones or abandoned wells can be greater. To our best knowledge, no experimental results of convection enhancement due to geochemical reactions are reported. This strongly suggests a potential future research. Figure 10 displays porosity changes, based on volume fraction deviations, for the case of Ra = 10,000 and heterogeneous reservoir. Though dissolution observed has significantly been boosted by the reactions, maximum porosity alteration until the reservoir becomes CO 2 saturated is only 0.02%. This slight increase is almost negligible and cannot affect hydrodynamics of the system by any means. By using porosity-permeability relation 29

Methods (reservoir model and governing equations)
Following the model development in Islam et al. 30 the continuity and transport equations read Here c, t, u, k, ψ represent concentration, time, velocity, permeability, and the stream function. x and z are coordinates in lateral and vertical (positive downward) directions. Ra is the primary parameter that defines buoyancy driven flow. The bottom, lateral boundaries are impervious, and the top dense CO 2 -brine interface is perturbed with saturated dissolved CO 2 . Numerical solutions of Eqs. 1 and 2 are explained in details in Islam et al. 31 .
As dissolved CO 2 comes into contact with rocks, it lowers the pH and triggers redistribution of carbonate aqueous species via the following reactions:

Conclusions
Convection from CO 2 -laden boundary is examined by coupling with reactive modeling of resident carbonate (calcite and dolomite) formations. Reactions involved were in locally equilibrium. We tested both homogeneous and heterogeneous cases with low and high mean permeability in term of Ra number. The sanity test delivered initial results of maximum equilibrium concentration of dissolved species. This showed dissolved CO 2 stream from reactions can increase local density even more than one order of magnitude. Geochemistry based on reaction intensity can cause serious alteration of flow dynamics. Our numerical results provide clear message that upon affecting local concentrations density-driven convection is enhanced substantially. For instance, Ra = 10,000 and heterogeneous distributions show that reservoir reaches ~100% CO 2 saturation in 500 years while only convective flow covers 40% upper area. Negligible porosity and permeability increase do not affect the hydrodynamics at all. The strong enhancement, in turn, adds negative impact on CCS in the sense that more free phase CO 2 may exist in the storage formations, increasing the potential for leakage and prolonging the period for long-term monitoring.