CO2 Dissolution Trapping Rates in Heterogeneous Porous Media

The rate of carbon dioxide (CO2) dissolution in saline aquifers is the least well‐constrained of the secondary trapping mechanisms enhancing the long‐term security of geological carbon storage. CO2 injected into a heterogeneous saline reservoir will preferentially travel along high permeability layers increasing the CO2‐water interfacial area which increases dissolution rates. We provide a conservative, first‐principles analysis of the quantity of CO2 dissolved and the rate at which free‐phase CO2 propagates in layered reservoirs. At early times, advection dominates the propagation of CO2. This transitions to diffusion dominated propagation as the interfacial area increases and diffusive loss slows propagation. As surrounding water‐filled layers become CO2 saturated, propagation becomes advection dominated. For reservoirs with finely bedded strata, ∼10% of the injected CO2 can dissolve in a year. The maximum fraction of CO2 that dissolves is determined by the volumetric ratio of water in low permeability layers and CO2 in high permeability layers.


Introduction
Worldwide carbon dioxide (CO 2 ) emission targets are unlikely to be met without large scale geological CO 2 storage (Intergovernmental Panel on Climate Change, 2018). Assessments of global CO 2 storage capacity (e.g., Michael et al., 2010) suggest that saline aquifers could account for around 90% of the total potential storage volume. The dissolution of injected CO 2 into the ambient brine within a saline aquifer is a key mechanism for increasing the security of long-term storage. At typical storage reservoir conditions, CO 2 is in the supercritical phase and is buoyant with respect to the surrounding reservoir fluid, presenting the risk of migration to the surface. As CO 2 dissolves into water the density of the water increases (Teng & Yamasaki, 1998), eliminating the buoyancy of free-phase CO 2 and reducing the risk of leakage. Quantifying total dissolution rates post injection is therefore important for assessing the contribution of CO 2 dissolution to the long-term security of stored CO 2 .
The rate of CO 2 dissolution in formation waters is controlled by the diffusive transport of dissolved CO 2 away from the CO 2 -water contact and the area of the contact. Because diffusive fluxes into a static system decrease as the square-root of time and the CO 2 diffusion coefficient is small (∼2 × 10 −9 m 2 s −1 (Cadogan et al., 2014)), CO 2 -enriched boundary layers in water in contact with free-phase CO 2 (Lindeberg & Wessel-Berg, 1997)

10.1029/2020GL087001
Key Points: • Injection of carbon dioxide into finely bedded reservoirs leads to enhanced contact area with water and hence enhanced dissolution rates • Propagation rates of free-phase carbon dioxide transition from advection to diffusion dominated as the contact area with water increases • For injection into the Salt Creek reservoir, Wyoming, nearly 10% of the injected carbon dioxide is predicted to dissolve in one year will grow to ∼10 cm thick in 1 year or ∼1 m in 100 years. The relative movement of water and CO 2 will therefore exert an important control on CO 2 dissolution rates. During injection and flow, the lower viscosity of the CO 2 will result in fingering (Saffman & Taylor, 1958) which will be strongly enhanced by reservoir heterogeneities. It is the resulting complexities in the geometry of the CO 2 -water interface and flow of CO 2 that makes CO 2 dissolution difficult to quantify.
There are few constraints on CO 2 dissolution during CO 2 injection. Measurements of CO 2 / 3 He ratios show that some natural CO 2 accumulations have lost more than 90% of their original CO 2 by dissolution over hundreds of thousands to millions of years (Gilfillan et al., 2009). At Green River, Utah, where natural CO 2 has been migrating up a fault system for several hundred thousand years, Bickle and Kampman (2013) estimated that less than 1% of the CO 2 escaped to the surface, the rest being dissolved in permeable horizons intersected by the fault system. Most modelling of CO 2 dissolution has concentrated on the impact of convective circulation of the brine beneath CO 2 accumulations driven by the density increase as brine saturates with CO 2 (e.g., Ennis-King & Paterson, 2005;Neufeld et al., 2010). However, the marked anisotropy of permeabilities in most reservoirs substantially reduces the convective circulation (Green & Ennis-King, 2014). Reservoir simulations using numerical models typically use grid sizes of more than 10 m which are unable to model flow heterogeneities on length scales of ∼1 m or less over which diffusion dominates. There are limited constraints on dissolution rates from measurements on small-scale injection experiments. In the Frio experiment, Texas, Freifeld et al. (2005) noted that the arrival times of the tracer krypton lagged the arrival of the tracers sulfur hexafluoride and perfluorocarbon and attributed this to the higher solubility of krypton in brine. Likewise Lu et al. (2012) observed a similar lag between sulfur hexafluoride and kyrpton tracers in the Cranfield, Mississippi CO 2 injection experiment. However, attempts to quantify such observations have had limited success (e.g., LaForce et al., 2014). In a CO 2 injection phase at the Salt Creek, Wyoming enhanced oil recovery site, Bickle et al. (2017) observed that dissolution of CO 2 in formation brines drove significant reactions with silicate minerals, but again, the difficulty in modelling the complex flows in a heterogeneous reservoir has so far precluded quantitative estimates.
The CO 2 -brine interactions which determine CO 2 dissolution will be controlled by the reservoir heterogeneities on all scales, and these are difficult to model properly, both because it is not possible to determine the reservoir structure at the submeter scales which matter for the diffusive processes, and because numerical models of CO 2 and brine flows in reservoirs are not currently capable of running at such resolutions.
In this paper, we consider dissolution during CO 2 injection into a simple representation of a layered reservoir and quantify the increased dissolution rates due to an increase in interfacial area between the CO 2 and the reservoir fluid. This provides a base case, given that the additional complexities are likely to substantially increase dissolution rates. The model is then evaluated using parameters appropriate to large-scale CO 2 injection such as reservoir bedding thickness, porosity, saturation, and injection flux.
We model a horizontally layered reservoir comprising alternating higher and lower permeability layers. Free-phase, low viscosity CO 2 flow will preferentially be confined to the higher permeability layers. This will be enhanced by capillary entry pressures which impede CO 2 entering the lower permeability layers (c.f. Sathaye et al., 2014). The model assumes horizontal strata within the reservoir with CO 2 flow confined to the high permeability layers and ignores buoyancy of the supercritical CO 2 . It is assumed that there is no flow of CO 2 between layers, but there is diffusive exchange of CO 2 across the static, water-filled low permeability layers. The modelling provides a minimum estimate for CO 2 dissolution against which the effect of additional processes or field observations may be assessed. Mixing of CO 2 and water along formation boundaries, fingering of low viscosity CO 2 into formation waters and the consequent dissolution of CO 2 ahead of the CO 2 finger, and the much more complex permeability structures in most sedimentary rocks would all be expected to enhance dissolution rates, most probably by an order-of-magnitude or more. It will likely only be possible to estimate the impact of these processes by experiments in field settings.

Single Layer Model
We first consider CO 2 propagating along a single, high permeability layer of width w, porosity , and permeability k in a low permeability, water saturated aquifer of porosity a , and permeability k a (Figure 1). CO 2 is injected into the high permeability layer at constant volumetric rate Q. For simplicity, we assume that a finite capillary entry pressure confines the flow of free-phase CO 2 to the high permeability layer. The CO 2 Figure 1. Schematic diagram of a two-dimensional finger of free-phase CO 2 propagating along a high permeability porous layer that is initially saturated with water surrounded by a low permeability porous aquifer also saturated with water.
finger has a total length L(t) the front moves with a flux V F (t). The length of the CO 2 finger is much greater than its width so dissolution across the CO 2 -water interface at the finger front is neglected, and the interface is modelled as planar for simplicity.
The diffusive CO 2 profile away from the CO 2 -water interface (in the z direction) for a given value of x is given by the solution for diffusion into a semiinfinite layer (Carslaw & Jaeger, 1959, p. 59), where erfc is the complimentary error function, c 0 is the maximum solubility of CO 2 in water, t is the time since injection commenced, and t 0 is the time at which the front passes position x = L(t 0 ). The effective diffusion coefficient of CO 2 in water is given by where D m is the molecular diffusion coefficient of CO 2 in water and a and are the porosity and tortuosity of the low permeability layer (Pismen, 1974).
The vertical CO 2 concentration gradient in the water is therefore and hence the vertical diffusive flux of CO 2 out of the high permeability layer at time t is, with a introduced as CO 2 only diffuses into water within the pores. As a nonwetting phase, CO 2 only partially displaces water in the high permeability layer. This reduces the fraction of the porosity occupied by CO 2 (s nw ), and some CO 2 dissolves in the water occupying the remaining pore space.
The velocity of the CO 2 front at x = L(t) is a function of the input flux and diffusive losses given by lateral diffusion from the finger and complete saturation of the residual water within the CO 2 finger, Here, v is the interstitial velocity of CO 2 at the front and v 0 is the CO 2 interstitial velocity at x = 0. By introducing the nondimensional variables Equation 5 may be rewritten in the generic form where = c 0 (1 − s nw )/s nw is a measure of how much CO 2 dissolves into residual water within the CO 2 finger. Equation 7 gives the dimensionless front velocity as a function of dimensionless time. For notational convenience, we drop the "∼" from all subsequent quantities. We solve Equation 7 numerically using a sequential iteration approach. The CO 2 concentration gradient at the interface is calculated every timestep allowing the total diffusive flux to be deducted from the input flux giving the CO 2 velocity, v, as a function of time. The new front velocity allows the position of the CO 2 front, L(t), to be calculated for the next timestep.

Figure 2a
illustrates the length of the CO 2 finger as a function of time. The length of the finger is governed by the input flux and the amount of dissolution. At early times (t ≪ 1), when the dissolution area to input flux ratio of the CO 2 finger is small, the amount of dissolution is negligible, and so the length evolves as L ∼ t. The diffusive CO 2 profile away from the CO 2 -water interface scales as t −1/2 ; hence, the total diffusive flux scales as t 1/2 . The total dissolution into the low permeability layers is the sum of the total flux over time and so scales as t 3/2 (see Figures S1 and S2 in supporting information for graphs of dissolution scaling).
At late times (t ≫ 1), the finger length evolves as L ∼ t 1/2 as diffusive loss dominates. This means that the total diffusive flux tends towards a constant, and hence the total dissolution evolves proportional to t. The transition between these two regimes happens when L and t are of the order one. Figure 2b shows the fraction of injected CO 2 lost by dissolution as a function of time. This fraction increases with time as the increase in the surface area increases the ratio of diffusive loss to input flux and tends to one as t → ∞. Figure 2 has been plotted for three values of . When = 0, there is no dissolution into residual water in the CO 2 finger. At low values of , the fraction of CO 2 dissolved in the residual water within the CO 2 finger is small compared with lateral loss to the surrounding water. For a typical reservoir, = 0.014 if s nw = 0.8 (Krevor et al., 2015) and c 0 = 5.5wt% (Dubacq et al., 2013).

Multilayered Model
The saline aquifers suitable for geological storage are characteristically sandstones bedded on 10 −2 to 1 m scales with permeabilities that vary by an order of magnitude or greater. Injection of CO 2 will primarily occupy the high permeability layers, and the diffusive fringes about the CO 2 -filled layers will overlap. We illustrate this behavior with a periodically layered reservoir with high permeability layers of width w, porosity and permeability k interbedded with low permeability layers of width 2h, porosity a , and permeability . Schematic diagram of periodically repeating high permeability porous layers of width w, separated by low permeability porous layers of width, 2h. The reservoir is initially saturated with water. CO 2 is injected into the high permeability layers. z = 0 at the high/low permeability interface and z = h at the midpoint between the high permeability layers. k a . CO 2 flows into each of the high permeability layers at volumetric rate Q. We assume that a finite capillary entry pressure confines advective flow of free-phase CO 2 to the high permeability layers (see Figure 3).
The diffusive profile between layers is given by the solution for diffusion into a layer bounded by two parallel planes (Carslaw & Jaeger, 1959, p. 100), (8) The velocity of the CO 2 front at x = L is a function of the input flux, lateral diffusive loss, and saturation of the residual water, Introducing the nondimensional variables Equation 9 can be rewritten as where = c 0 (1 − s nw )/s nw and = h a c 0 w s nw . For notational convenience, we drop the "∼" from all subsequent quantities. A similar iterative solution to the single finger case is used to solve Equation 11 to give the CO 2 front position and velocity, and total dissolution of CO 2 as a function of time for the multilayered model. Figure 4a illustrates the length of the CO 2 fingers as a function of time for = 5 and three values of .

Results
The length of the fingers evolves in three stages. At early times (t ≪ 1) the length of the fingers evolves as L ∼ t, at intermediate times they evolve as L ∼ t 1/2 , while at late times (t ≫ 1) they grow linearly L ∼ t. The total lateral dissolution of CO 2 in the early advection dominated regime scales as t 3/2 . The transition between the early time regime and the intermediate, diffusion dominated regime occurs as the increase in the surface area drives increasing dissolution. The transition from the intermediate diffusion dominated regime to the late-time advection dominated regime is due to CO 2 saturation of water in the low permeability layers which dampens diffusion over the more proximal parts of the CO 2 layers. The system tends to a steady state with a constant length zone at the front of the CO 2 finger in which dissolution of CO 2 is significant. The total dissolution scales proportionally with time in this late time regime (see Figures S1 and S2 for graph of dissolution scaling). Figure 4b shows the fractional loss of CO 2 by dissolution as a function of time. At late times, the fractional loss of CO 2 tends to a constant value which is less than one. The value of (= h a c 0 /w s nw ) determines the transition time between the three regimes, where is a ratio between the volume available for CO 2 to dissolve into the low permeability layers (h a c 0 ) and the volume of CO 2 in the high permeability layers (w s nw ). Larger values of allow more CO 2 diffusion leading to later transition to the late-time advection dominated regime ( Figure 5a). Importantly, also determines the maximum fraction of CO 2 dissolved at long times, reflecting the ratio of the volume of water available for saturation with CO 2 and the volume of CO 2 in the high permeability layers (Figure 5b).

Discussion
Evaluating the model using parameters appropriate to field settings establishes a practical sense of the permeability structures, lengths, and timescales for which significant dissolution of CO 2 will occur. The Salt Creek Oil Field in Wyoming has been the site of CO 2 injection for enhanced oil recovery since 2004. In 2010, there was a monitored injection of CO 2 into a 20 m interval of the second Wall Creek sandstone unit. This is a highly heterogeneous deltaic sequence made up of mudstones, siltstones, and sandstones in coarsening up sequences (Lee et al., 2005). Bickle et al. (2017) estimated the permeability profile of the injection interval using porosity measurements calculated from gamma ray density logs. Order of magnitude permeability variations were found on ∼0.5 m lengthscales. However, the resolution of the permeability distribution was limited by the resolution of the gamma ray density logs, which was ∼ 0.35 m, and it is probable that large variations in permeability on smaller length scales exist.  The periodically repeating layered model is evaluated using parameters from the CO 2 injection into the second Wall Creek sandstone unit. The parameters used in the calculation are effective diffusivity D = 2 × 10 −11 m 2 s −1 using D m = 2 × 10 −9 m 2 s −1 (Cadogan et al., 2014), a = 0.12 and = 5, CO 2 input velocity v 0 = 4 × 10 −5 ms −1 and porosity = 0.2 (Bickle et al., 2017), fraction of the porosity occupied by CO 2 s nw = 0.8 (Krevor et al., 2012), maximum saturation concentration c 0 = 5.5 wt% calculated at 15 Mpa, 50 • C and a salinity of 0.05 mol NaCl/kg(H 2 O) (Dubacq et al., 2013), and a high permeability layer spacing to thickness ratio h/w = 1.5 (Bickle et al., 2017). Three different widths for the high permeability layer have been plotted. These are w = 0.5 m, as calculated by (Bickle et al., 2017), as well as w = 0.1 m and w = 0.05 m, accounting for the limited resolution of the permeability distribution. Figure 6a shows the length of the CO 2 finger from the injection point as a function of time. A separate curve is plotted for each value of w, with the ratio between high permeability layer spacing to layer width held constant. The length of the finger if no diffusive loss occurs is also plotted (black dashed line). As the width of the high permeability layer decreases, that is, thinner and more finely spaced bedding, the distance the CO 2 propagates into the reservoir in a given time decreases becoming more pronounced at later times. When w = 0.05 m, the CO 2 travels around 10% less far than if no dissolution had occurred.
It is also useful to know how much of the injected CO 2 dissolves into the surrounding water. Figure 6b shows the fraction of the total injected CO 2 that has dissolved into the surrounding water as a function of time. The high permeability layers of width 0.1 and 0.05 m show total dissolution of around 5% and 9% of the total injection volume respectively within the first two years of injection. Less CO 2 dissolves if bedded layers are thicker.

Conclusion
Injecting CO 2 into saline reservoirs with interbedded high and low permeability layers substantially enhances dissolution rates. As the fluid travels further into the reservoir, the increase in surface area between the CO 2 and surrounding water causes increased diffusive loss. The velocity at which CO 2 travels in the reservoir is dominated by the advective input flux at early times and transitions to an intermediate diffusion dominated regime as diffusive loss increases. At late times, the water in the low permeability layers reaches CO 2 saturation, dampening diffusion and resulting in a return to an advection dominated regime. The transition times between these regimes is governed by the ratio between the volume available for CO 2 dissolution in the low permeability layers and the volume of CO 2 within the high permeability layers. This ratio also governs the maximum fraction of injected CO 2 dissolved at late times. In reservoirs with characteristic bedding thicknesses of ∼0.1 m, the modelling implies that a significant fraction of the CO 2 will dissolve in water within a few years of injection. The tendency of low viscosity supercritical CO 2 to finger and the much more complex flow paths in real reservoirs will likely increase CO 2 dissolution rates above the minimum estimates from this model.

Data Availability Statement
The data from numerical simulations in this study are in the repository (https://doi.org/10.17863/CAM. 47888).