Geological Carbon Sequestration by Reactive Infiltration Instability

We study carbon capture and sequestration (CCS) over time scales of 2000 years by implementing a numerical model of reactive infiltration instability caused by reactive porous flow. Our model focuses on the mineralization of CO2 dissolved in the pore water—the geological carbon sequestration phase of a CCS operation—starting 10–100 years after the injection of CO2 in the subsurface. We test the influence of three parameters: porosity, mass fraction of the Ca-rich feldspar mineral anorthite in the solid, and the chemical reaction rate, on the mode of fluid flow and efficiency of CaCO3 precipitation during geological carbon sequestration. We demonstrate that the mode of porous flow switches from propagation of a planar front at low porosities to propagation of channels at porosities exceeding 10%. The channels develop earlier for more porous aquifers. Both high anorthite mass fraction in the solid phase and high reaction rates aid greater amounts of carbonate precipitation, with the reaction rate exerting the stronger influence of the two. Our calculations indicate that an aquifer with dimensions 500 m × 2 km × 2 km can sequester over 350 Mt solid CaCO3 after 2000 years. To precipitate 50 Mt CaCO3 after 2000 years in this aquifer, we suggest selecting a target aquifer with more than 10 wt% of reactive minerals. We recommend that the aquifer porosity, abundance of reactive aluminosilicate minerals such as anorthite, and reaction rates are taken into consideration while selecting future CCS sites.


INTRODUCTION
Carbon capture and permanent sequestration in the subsurface is becoming an important mechanism in carbon neutralization (Hosa et al., 2011;Odenberger et al., 2013;Celia, 2017). Current anthropogenic CO 2 emissions (42 Gt/yr, Le Quéré et al., 2018;Friedlingstein et al., 2019) are three orders of magnitude higher than the current combined global CCS injection rate (35-40 Mt/yr, Celia, 2017). Recognizing the importance of CCS in mitigating global carbon emissions, a number of different incentives are offered. For example, the 2018 US Bipartisan Budget Act (26 U.S.C. § 45Q) provided a tax credit of $20 per ton of CO 2 sequestered permanently. To encourage CCS using existing infrastructure, a number of locations-Sleipner site in the Norwegian North Sea (1 Mt/yr, Hosa et al., 2011), Kimberlina site in the western US (0.25 Mt over 4 years, Doughty, 2010), and Goldeneye reservoir in the North Sea (10 Mt CO 2 over 10 years, Spence et al., 2014)-are either targeted for future operations or are being used as permanent CO 2 storage sites. Hosa et al. (2011) provided an overview of the storage capacity, reservoir characteristics, and injection rates of 20 active CCS sites over the world. This article identifies a number of controls on the efficiency of CCS in potential future sites and testing facilities. Kingdon et al. (2019) provided details of such planned experimental facilities at the United Kingdom Geoenergy Observatory Cheshire Energy Research Facility Site, which will be used for in-situ testing of these controls.
Carbon capture and storage can broadly be envisioned as a two-step process (Bachu, 2000;Friedmann, 2007;Jiang, 2011;Aminu et al., 2017). During the early phase, supercritical CO 2 is injected into the subsurface (Doughty, 2010;Bacci et al., 2011;Emami-Meybodi et al., 2015). Following subsurface injection, the supercritical CO 2 remains a separate, migrating phase before slowly reacting with the connate water in the saline aquifer. This slow phase of reaction with saline groundwater leads to the formation of carbonic acid by the reaction.
We assume that the target aquifer has minimal buffering capacity, which leads to lower pH values. Therefore, given the time frame of the model used in this paper, H 2 CO 3 will be dominant over bicarbonate species. Once the solution of CO 2 is completed, the process of geological carbon sequestration (GCS), or mineralization of carbon, ensues. During this phase, dense, carbonic acid-rich groundwater descends into the porous host formation, reacting with the aluminosilicate and silicate minerals to form solid carbon deposits in the form of carbonate minerals, which takes place over time scales of 1,000 years. In this article, we focus on the mechanisms involved in this stage.
The process of carbonate mineralization by assimilation of dissolved H 2 CO 3 takes place through a combination of densitydriven porous flow, chemical reaction, and diffusion. The combination of these processes leads to a mechanism known as the reactive infiltration instability (RII) (Chadam et al., 1986), which has been associated with subsurface cavern or "wormhole" formation (Hinch and Bhatt, 1990;Szymczak and Ladd, 2013;Szymczak and Ladd, 2011) and chemical reaction between magmatic melts and mantle rocks (Aharonov et al., 1995;Spiegelman et al., 2001;Takei and Hier-Majumder, 2009;Sun et al., 2020a). An important manifestation of RII is the propagation of the reactive fluids into the aquifer through transient "channels" or "fingers", which merge and interact with each other during the penetration of the reactive fluid into the saline aquifer (Riaz et al., 2006;Soltanian et al., 2016). Laboratory experiments on reactive porous flow of basaltic melt into olivine (Pec et al., 2015) and dense brine injection into a porous medium also support this observation (Kneafsey and Pruess, 2010;Vosper et al., 2014).
In addition to the flow regime, the physical and chemical characteristics of the saline aquifer or subsurface reservoir also exert an influence on the efficiency of GCS. Previous models of porous flow and laboratory experiments studied the efficiency of mass transfer as a function of two important dimensionless numbers: the Damköhler number (Da), representing the ratio between the rates of reaction and advection, and the Péclet number (Pe), representing the ratio between rates of advection and diffusion (Chadam et al., 1986;Riaz et al., 2006;Ghesmat et al., 2011;Szymczak and Ladd, 2013). These studies, however, ignored the crucial role of parameters such as the lithology and porosity of the aquifer on the efficiency of GCS. These parameters can be potentially important in identifying CCS sites and quantifying the amount of long term carbon sequestration in currently active sites. They can also vary significantly across different sites. For example, Hosa et al. (2011) showed that currently active sites of CCS can be characterized by a low porosity of 3% and a permeability of 0.1 mD (MRCSP R.E. Burger project, Ohio, USA) to a high porosity of 37% and a permeability of 5 D (Sleipner site, Norwegian N. Sea). The role of porosity and permeability variations is generally incorporated into models of supercritical CO 2 injection over periods of time spanning over a few years (Class et al., 2009;Doughty, 2010;Soltanian et al., 2016). Systematic analyses of the efficiency of GCS as a function of porosity, formation lithology, and fundamental dimensionless numbers still remain scarce.
In this work, we bridge this gap by studying the influence of formation porosity and the mineral composition on the efficiency of GCS in a porous saline aquifer. Our model begins after the supercritical CO 2 has been injected into the aquifer and has reacted with the brine in the pore space to create a dense, H 2 CO 3 rich pore fluid. In our numerical simulations, this dense pore fluid sinks into the unaffected formation by reactive porous flow. To model the deposition of solid carbonate resulting from the reaction between the H 2 CO 3 rich fluid and the solid matrix, we carry out a series of 2D finite element simulations using a massively parallel open source software MuPoPP 1.2 (Hier-Majumder, 2020) running on the Oracle cloud computing platform. In these simulations, we vary the dimensionless Damköhler (Da) and Péclet (Pe) numbers, aquifer porosity, and the initial concentration of the Ca bearing Da feldspar mineral phase anorthite (An) in the aquifer. As we demonstrate, these parameters influence (a) the mode of fluid transport by a transition from propagation of a planar front to channelized flow and (b) the amount of carbonate deposited in the solid phase after 2000 years of fluid flow.

Governing Equations and Boundary Conditions
We model the porous flow of a reactive, H 2 CO 3 rich fluid as a reactive infiltration instability problem. As discussed above, this model focuses on the long term GCS process and begins after the supercritical CO 2 has been pumped into the subsurface and equilibrated with the groundwater to create a dense, H 2 CO 3 rich pore fluid. The density of the pore fluid, ρ, depends on the concentration of the dissolved H 2 CO 3 in the pore water, c 0 , and is related to the concentration by ρ ρ 0 (1 + cc 0 ), where ρ 0 is the density of H 2 CO 3 free water, and c is a dimensionless constant. Driven by the concentration dependent density, pore fluid from the aquifer percolates through the porous medium reacting with the aluminosilicate minerals to precipitate CaCO 3 . In our model, we use the chemical reaction between anorthite (An CaAl 2 Si 2 O 8 ) in the solid with H 2 CO 3 in the pore fluid, leading to the precipitation of CaCO 3 and kaolinite (Ka Al 2 Si 2 O 5 (OH) 4 ), given by the reaction.
While a number of chemical reactions can take place between H 2 CO 3 and subsurface minerals, the most common set of reactions involves aluminosilicates in the rock (Ghesmat et al., 2011;Zhang and Song, 2014). In this article, we focus on the influence on GCS driven by this chemical reaction between the solid and the pore fluid. We discuss the assumptions made in our model in detail in Section 2.2. We provide the full set of dimensional equations, nondimensionalization schemes, and details of the numerical methods in the Supplementary Material. For brevity, here we present the final, dimensionless equations. In this set-up, the velocity, u, of the pore fluid is governed by Darcy flow, where ϕ is the constant porosity of the rock, k (ϕ/ϕ 0 ) 3 is the permeability, ϕ 0 0.1 is a reference porosity, p is the fluid pressure, and z is a unit vector in the vertically upward direction (Sun et al., 2020a). The last term arises from the concentration dependent density of the fluid, where we set c 1. In addition to Darcy flow, the continuity equation requires that the velocity is nondivergent.
The rate of transport of dissolved H 2 CO 3 in the pore fluid is controlled by advection via porous flow, chemical diffusion, and chemical reaction, where X 0 is the stoichiometric mass fraction of H 2 CO 3 in Eq.
(2), and the dimensionless quantities Pe and Da are defined below in Eq. 8. We consider a second order chemical reaction driven by the rate of the forward reaction (Ghesmat et al., 2011). As a result, the rate of consumption and production of reactants and products depend on the concentrations of H 2 CO 3 in the pore fluid, c 0 , and mass fraction of anorthite in the system, c 1 , as evidenced by the second term on the right hand side. Finally, the rate of consumption of anorthite in the rock and the rate of production of CaCO 3 are given by the two equations where X 1 and X 2 are the mass fractions of anorthite and CaCO 3 in the chemical reaction (Eq. 2), and c 2 is the mass fraction of CaCO 3 produced by GCS.
The dimensionless quantities are defined by the relations where u 0 is a characteristic fluid velocity, H is the height of the domain (thickness of the aquifer), D is the chemical diffusivity, and Γ 0 is a characteristic rate of chemical reaction. In this work, we chose to carry out a number of simulations over a range of plausible values of these dimensionless quantities, which allows us to quantify the influence of not only the dimensionless numbers, but also the aquifer characteristics on the efficiency of GCS. We discuss the rationale behind selecting these ranges in the following subsection. We solved the five governing Eq. 3-7 for the five unknowns-u, p, c 0 , c 1 , and c 2 -subject to a set of Dirichlet and Neumann boundary conditions (BCs) and initial conditions (ICs). These BCs and ICs are illustrated in the schematic diagram in Figure 1. On the top boundary, zΩ t , we assume a Dirichlet boundary condition, namely a uniformed constant downward velocity w 0 z. On the vertical boundaries, zΩ v , the fluid velocity and concentrations are defined to be periodic. On the top and bottom boundaries, a Neumann boundary condition is applied on zΩ N .
In addition, we impose the initial conditions, where [An] is a constant value of initial anorthite mass fraction in the aquifer. In order to explore the influence of aquifer composition on the efficiency of the GCS, we carried out a number of simulations with different values of [An]. The unknown c 2 (x, z, t) gives the concentration of CaCO 3 at each point within the domain at a given time. It is useful to also have an estimate for the mass of CaCO 3 that can be sequestered by the process. To calculate the cumulative CaCO 3 mass M over time, we use the definition of production rate of CaCO 3 in Eq. 7. Using the right hand side of the equations from the known quantities at each time step and integrating, we get where M(t) is the transient mass of CaCO 3 deposited over the time period 0 − t within the domain Ω (the aquifer), of height H and depth L. Since our simulations are 2D, we assume that the length of the domain in the third dimension is also L. Finally, the density of sandstone in the aquifer is given by ρ sst .
The dimensional values of these constants are listed in Table 1.

Model Assumptions
In this model, we focus on the influence of three parameters: aquifer porosity, initial abundance of Anorthite in the aquifer, and the chemical reaction rate on the structure of the fluid percolation and the amount of solid carbonate precipitation.
To address the complex, nonlinear interplay among these factors, we made a few assumptions in our model. First, we only consider the chemical reaction given in Eq. 2 and test the outputs of our models over a range of Da, a method followed by previous numerical studies (Aharonov et al., 1995;Ghesmat et al., 2011). In addition to the reaction between H 2 CO 3 rich fluid and anorthite, reactions between the fluid and CaCO 3 and, although much slower, reactions between the fluid and quartz are possible (Espinoza et al., 2011). By ignoring these possible reactions, especially the dissolution of CaCO 3 , we assume that the pH level of the pore fluid is buffered such that the fluid is saturated with CaCO 3 . While we do not explicitly include such buffering reactions in our model, Snaebjörnsdóttir et al. (2017) reported that samples of injected fluid reached CaCO 3 saturation within months of injection at the CarbFix site in Iceland. Given the substantially longer timescale of our simulations, this observation of Snaebjörnsdóttir et al. (2017) suggests that our assumption of CaCO 3 buffering in the pore fluid is likely reasonable.
Second, we ignore the influence of calcite precipitates on reducing the porosity and permeability of the aquifer. Within the timescale of our simulations, this effect is likely to play an insignificant role. As we show in Figure 6, the highest concentration of precipitated CaCO 3 in our simulations is less than 0.5%. Recent microtomographic studies on permeability reduction by cementation show that porosity and permeability reduction by precipitates is negligible for cement volume fractions less than 3-5% ( Thomson et al., 2019;Thomson et al., 2020). The assumption of constant porosity, therefore, is reasonable for the range of time considered in these simulations. With these caveats, we next discuss the range of parameters explored in this work.

Parameter Ranges
In this study, we focus on the influence of a number of factors on the modes and efficiency of GCS by reactive infiltration instability. These parameters are highlighted on Table 1. As discussed below, it is often difficult to assign a particular value of a characteristic dimensionless number such as Da and Pe. In addition, aquifer features such as porosity and the anorthite content can vary between sites. It is, therefore, useful to have the information on the influence of these factors on the efficiency of GCS, which can be used for selection and planning of sites such as those associated  with the UKGEOS initiative. In this section, we discuss the range of values of these parameters discussed in literature and the rationale behind the ranges selected in this study. Despite a number of reports on the second order reaction rate for dissolution of silicate and aluminosilicate minerals (Oelkers and Schott, 1995;Palandri and Kharaka, 2004;Espinoza et al., 2011), direct estimates of Da from experimental measurements remain scarce. Here, we present a range of most likely natural values of Da relevant to the CCS problem. Following previous works on reactive infiltration instability (Steefel and Lasaga, 1994;Aharonov et al., 1995), we notice that the rate of mass fraction Γ 0 can be expressed as a function of dissolution rate, R, specific surface area A, and the density of anorthite ρ ′ 2750 kg/m 3 , as where M An 0.28 kg/mol is the mass of anorthite per mole. Using the dissolution rate of anorthite from Espinoza et al. (2011), we calculate a lower and an upper limit of R 4.3 × 10 − 12 and R 1.3 × 10 − 10 mol/m 2 /s, respectively. Following the argument outlined by Steefel and Lasaga (1994), we take the specific surface area as the area of contact between the pore fluid and grains per unit volume. We use the microstructural model of Hier-Majumder and Abbott (2010) to obtain this specific surface area (unit of m 2 /m 3 ). Next, we combine the value of A with reaction rate R in Eq. 14, yielding values of Γ 0 ranging between 1.0 × 10 − 11 and 3.1 × 10 − 10 s −1 . Using the values of characteristic velocities and length scales listed in Table 1, these lead to values of Da ranging between 0.3 and 9.7. We ran the simulations for a more conservative range of Da ranging between 0.1 and 5. It is possible that the value of Da is higher than our estimates, in this case the mass of deposited CaCO 3 will be higher than the amount reported here. For the range of Pe, we used a range of values between 0, completely diffusion dominated transport, to 10 4 , convection dominated transport (Chadam et al., 1986;Hinch and Bhatt, 1990). Of these two regimes of flow, the latter, convection dominated mode is more dominant in natural reactive infiltration instability problems (Chadam et al., 1986;Hinch and Bhatt, 1990;Riaz et al., 2006;Ghesmat et al., 2011;Soltanian et al., 2016). Consequently, we show the results from a run of Pe 500 in this article. Using the definition of Pe from Eq. 8 and using the values of the dimensional constants in the table, we get Pe 3963. Thus, the results shown here depict a more conservative estimate of the flow rate than natural conditions. Finally, we use a range of values of aquifer porosity and anorthite concentration in the solid, as shown in Table 1.

RESULT
Our simulations identify the modes of penetration of dense H 2 CO 3 rich fluid into the porous aquifer as a function of various parameters controlling the process. Some of these factors, such as porosity and aquifer composition, are potentially critical in GCS site selection. In Figures 2-4 we show sets of three snapshots of the concentration map of H 2 CO 3 dissolved in the pore water (c 0 (x, z, t)). To isolate the effect of one controlling parameter in each figure, we kept the other three constant and show the concentration map for four different values of the controlling parameter. A number of previous studies on reactive infiltration instability focused on the role of Pe and Da (Aharonov et al., 1995;Riaz et al., 2006;Ghesmat et al., 2011;Szymczak and Ladd, 2013;Sun et al., 2020a), as we show here, the aquifer properties such as porosity and composition exert equally important controls on the mode and efficiency of GCS by reactive infiltration instablity. The source code for the simulation, MuPoPP1.2.0 (Multiphase Porous flow and Physical Properties) is publicly available . Simulation data is also available for download through the Royal Holloway Figshare repository (Sun et al., 2020b).

The Effect of Porosity
Rock porosity exerts a strong influence on the mode of invasion of H 2 CO 3 rich pore fluid into the aquifer, by controlling the shape of the reaction front. As outlined in Figure 2, for small values of porosity (ϕ 0.01 or 0.05), the H 2 CO 3 rich pore fluid descends slowly, as a planar front. In contrast, simulations with higher values of porosity (ϕ 0.10 or 0.15) depict the propagation of H 2 CO 3 rich fluid as growing channels, characteristic of reactive infiltration instability. We also notice that the channels for ϕ 0.15 begin forming earlier, grow faster, and are characterized by the thinnest top boundary layer for all four sets of simulations reported here. The rate of infiltration of the H 2 CO 3 rich fluid into the aquifer is also much faster with the presence of channels. The snapshots in Figure 2 illustrate that while the H 2 CO 3 rich channels penetrate the aquifer (500 m) by 2000 years for ϕ 0.15, the planar interface traverses only half of the aquifer after the same time for ϕ 0.01. We also notice that for ϕ 0.15, the onset of channelization is marked by a large number of small wavelength channels, which later merge into fewer, larger wavelength channels. This behavior of reactive infiltration instability arises from the fundamentally nonlinear nature of the process and interaction between channels (Aharonov et al., 1995;Riaz et al., 2006) and indicates that-at least within the parameter space explored in this work-GCS is dominated by processes far from equilibrium. Vosper et al. (2014) tested the formation of channels in a nonreactive Hele-Shaw cell experiment containing close packed glass beads. The porosity of close-packed systems, typically exceeding 25% (Wimert and Hier-Majumder, 2012), is much higher than the range tested here. Qualitatively, their observation of channel growth and low spacing between the channels are similar to the trend we report at higher porosities.
While previous studies of reaction infiltration instability noticed the formation of channels (Hinch and Bhatt, 1990;Aharonov et al., 1995;Ghesmat et al., 2011;Sun et al., 2020a) and evolution of the channels (Riaz et al., 2006;Soltanian et al., 2016), none of these studies reported the transition of the flow from a slow planar front to a highly nonliner channeling instability with an increase in the porosity of the aquifer. The transition of the flow model from a planar front to channels indicates that GCS models need to incorporate the nonlinear effects arising from reactive Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 533588 5 FIGURE 3 | The effect of Anorthite concentrations for simulations with Da 5.0, Pe 500, and ϕ 0.1. Each of the four columns is marked by a different initial anorthite content in the solid, annotated at the top. The values of dimensional time at each step are annotated on the right. The snapshots show that a higher mass fraction of anorthite in sandstone weakens the propagation of H 2 CO 3 .

FIGURE 4 |
The effect of reaction rates for simulations with ϕ 0.1, Pe 500, and [An] 0.1. The snapshots show that efficient reactions (higher Da) will constrain the propagation of dissolved CO 2 , namely H 2 CO 3 in groundwater. Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 533588 6 infiltration instability. The spatial distribution of fluid flow rate is apparently heterogeneous, especially with the formation of channels. Therefore, if we use a parameterized, one-dimensional fluid flow rate when calculating CaCO 3 precipitation, the efficiency of the GCS can be significantly under or over estimated.

The Effect of Anorthite Concentration and Reaction Rate
The extent of reaction between H 2 CO 3 rich pore fluid and the aquifer matrix depends on both the availability of reactants (in this case, An) and the rate of chemical reactions. These two factors, while controlled by two different parameters, have a similar effect on the channel propagation.
A higher anorthite concentration [An] in the solid leads to more rapid consumption of the dissolved H 2 CO 3 in the pore fluid, as observed in the series of simulations with increasing anorthite content in Figure 3. With an increasing consumption of the H 2 CO 3 in the pore fluid, formation of the channeling instabilities at the reaction front are more subdued, as evidenced during each time step of the comparative snapshots. We observe the strongest channel formation after 2 ka in the simulations with [An] 0.01. In contrast, the simulations with [An] 0.15 (15 wt% anorthite in the initial solid composition) display weak channels containing lower concentrations of the dissolved H 2 CO 3 as evidenced by the colormap. While these channels penetrate into the aquifer, the dissolved H 2 CO 3 in the channels are quickly consumed by the high concentration of anorthite in the surrounding, unreacted rock. As a result, the amount of CaCO 3 precipitation is much higher in these cases, an issue we discuss in Section 4.1.
High reaction rates, represented by high values of Da, also reduce the concentration and the propagation rate of H 2 CO 3 rich pore fluid. The images in Figure 4 display that in the simulations with higher Da, the concentration of H 2 CO 3 in the channels is lower compared to simulations with a lower Da. Similar to high initial anorthite abundance, a high reaction rate also aids rapid consumption of the acidic pore fluid and CaCO 3 precipitation. Although arising from two separate mechanisms, the influence of initial anorthite abundance in the solid and the rate of chemical reaction both exert a similar kind of control on the propagation of channels during GCS.

DISCUSSION
The results presented in Section 3 demonstrate the influence of the control parameters on the structure of the RII. In this section, we summarize the influence of porosity, initial abundance of anorthite, and the rate of reaction on the mass of CaCO 3 precipitated during GCS. For an aquifer with a given porosity, the reaction rate exerts a stronger influence than the initial anorthite concentration, while an increase in any of the three parameters leads to an increase in the amount of carbonate mineralization. Based on the structure of the porous flow and magnitude of carbonate mineralization, we present our recommendations for selecting future sites of CCS.

Capacity of Deep Saline Aquifers as GCS Sites
An important outcome of numerical modeling of GCS is an estimate for the amount of carbon permanently sequestered in solid mineral phases. The formulation in Eq. 13 provides the dimensional form of this time-dependent mass. As discussed earlier, we integrate the CaCO 3 deposition rate over a volume of 500 m thick and 2 km wide aquifer over 2000 years. The resultant values of permanently sequestered CaCO 3 mass are strongly Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 533588 dependent on the parameters studied in this work, as illustrated by the plots in Figures 5A,B. Aquifer properties characterized by initial anorthite mass fraction and porosity strongly influence the sequestered CaCO 3 mass. The plot in Figure 5A shows that this mass increases with the initial mass fraction of anorthite in the solid. Moreover, the trend of the growth is greatly modulated by the value of the Da, with a more strongly nonlinear trend evident at higher values of Da. Of the four cases shown in Figure 5A, the flattest curve corresponding to Da 0.1 shows that the CaCO 3 mass varies approximately between 1 and 15 Mt, while for a value of Da 5, the same increase in An concentration in the solid leads to a nonlinear increase of mass from 30 to 390 Mt. The comparison between these two influencing factors-chemical reactivity (Da) and abundance of reactant [An]-demonstrates that the chemical reactivity exerts a stronger influence on the amount of carbon permanently sequestered in minerals.
Porosity of the aquifer exerts a similar, but slightly more moderate influence on the mass of precipitated CaCO 3 . The plot in Figure 5B demonstrates this influence for values of Da exceeding 0.1. For Da 1, the mass increases from 35 to 60 Mt by increasing the porosity from 0.01 to 0.2. For a value of Da 5, this effect is more pronounced, as the mass increases from 65 to 260 Mt for the same increase in aquifer porosity. Porosity influences the flow by controlling the magnitude of permeability. A more porous aquifer permits faster transport of H 2 CO 3 rich fluid, which, subject to the reaction rate, can increase the precipitation of CaCO 3 .
The cumulative mass fraction of CaCO 3 sequestered by GCS varies with time. The plot in Figure 6 shows the evolution of the volume-averaged CaCO 3 mass fraction with time for four different aquifers. The top curves, shaded in gray, correspond to an aquifer containing 20% anorthite in the bulk, and the upper and lower limits are fixed by porosity values of 15 and 5% respectively. The bottom set of curves correspond to an aquifer containing 10% anorthite in the bulk. While following the conclusions drawn from the previous results, this plot also outlines that the process of GCS starts to produce appreciable CaCO 3 deposition at least after 500 years, but gains pace with time. This increase in rate is related to the penetration of the channels to the bottom of the aquifer and reacting with a larger amount of unreacted anorthite.

Geological Criteria for Future CCS Site Selection
As discussed by Hosa et al. (2011) and Celia (2017), current CCS sites contain rocks with a wide range of porosity values. To meet the challenge of carbon neutralization, and increase the current global injection rate of ∼ 40 Mt/yr (Celia, 2017), more CCS sites need to be identified. In identifying these sites, a number of factors such as the existence of a cap rock, cost of CO 2 transport to the site, and public perception need to be weighed against the geological criteria. Keeping these mitigating factors in mind, we use the results from our study to suggest parameter ranges that can ensure efficient GCS by CaCO 3 precipitation.
The amount of immobilized CO 2 , namely precipitated CaCO 3 , is influenced by initial anorthite amount and the reaction rate, as shown in Figure 5A. For a moderate reaction rate (Da 1), in order to sequester at least 50 Mt CaCO 3 after 2000 years, we suggest choosing an aquifer with more than 10 wt% of reactive minerals, such as anorthite. However, when the reaction rate is high (Da 3 − 5), an amount of 5 wt% anorthite in the target formation is sufficient to sequester 50 Mt CaCO 3 after 2000 years.
We observe that in deep saline aquifers containing > 10% porosity, the mode of fluid propagation shows the transition from a planar front to channels. Our simulations indicate that the channels are faster and more efficient in transferring H 2 CO 3 rich pore fluid into the aquifer, and as the curves in Figure 6 demonstrate, channels also lead to larger amounts of CaCO 3 precipitation. Thus, to ensure GCS by channel formation, we recommend selection of future sites containing > 10% porosity. Additionally, CaCO 3 precipitation increases with an increase in anorthite concentration in the solid. Comparison of experimentally determined reaction rates indicates that anorthite or anorthite-rich members of the feldspar solid solution group are more reactive than the albite end member (Palandri and Kharaka, 2004). If the composition of the aquifer rock is known, we recommend selecting an aquifer containing anorthite or other reactive aluminosilicate minerals (Oelkers and Schott, 1995;Espinoza et al., 2011) over one containing primarily unreactive quartz.

CONCLUSION
Our simulations demonstrate the presence of two distinct regimes of porous flow: planar front propagation at Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 533588 8 porosities less than 10% and channel formation at higher porosities. When comparing percolation rates into the aquifer, H 2 CO 3 rich pore fluids percolate faster along channels than planar reaction fronts. The amount of CaCO 3 precipitation by GCS grows with an increase in the reaction rate or the initial anorthite concentration of the aquifer. Of these two, the former exerts a stronger influence on the mass of CaCO 3 deposited after 2000 years. Under the most reactive conditions considered in this study, up to 390 Mt CaCO 3 is sequestered in the aquifer over a period of 2000 years. Also in order to sequester at least 50 Mt CaCO 3 after 2000 years, the amount of reactive minerals should exceed 10 wt%. While an increase in porosity also leads to an increase in the amount of CaCO 3 deposition, the effect of the reaction rate is stronger than the influence of porosity on the cumulative mass of CaCO 3 . Based on the results of the numerical models, we suggest the selection of future GCS sites containing > 10% porosity and characterized by the presence of anorthite or similarly reactive aluminosilicate minerals in the solid.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author. The source code for the simulation, MuPoPP1.2.0 (Multiphase Porous flow and Physical Properties) is publicly available . Simulation data is also available for download through the Royal Holloway Figshare repository (Sun et al., 2020b).

AUTHOR CONTRIBUTIONS
YS, RP, and SH-M developed the source code. YS carried out the simulations. AK provided information about the parameters around the planned UKGEOS observatory. All authors contributed toward writing the manuscript.

FUNDING
The computational resources for this work was provided by a computational grant from Oracle Research. RP was funded by the NERC London Doctoral Training Partnership, grant number NE/ L002485/1 and matching case funding grant from the BGS.