Multi-Layered Systems for Permanent Geologic Storage of CO 2 at the Gigatonne Scale

The effectiveness of Carbon Capture and Storage (CCS) as an imperative decarbonization technology relies on the sealing capacity of a fine-grained caprock to permanently store CO 2 deep underground. Uncertainties in assessing the caprock sealing capacity increase with the spatial and temporal scales and may delay CCS deployment at the gigatonne scale. We have developed a computationally efficient transport model to capture the dynamics of basin-wide upward CO 2 migration in a multi-layered setting over geological time scales. We find that massive capillary breakthrough and viscous flow of CO 2 , even through pervasively fractured caprocks, are unlikely to occur and compromise the storage security. Potential leakage from the injection reservoir is hampered by repetitive layering of overlying caprocks. This finding agrees with geologic intuition and should be understandable by the public, contributing to the development of climate policies around this technology with increased confidence that CO

The density contrast of the injected CO 2 relative to resident brine promotes buoyant upward flow of CO 2 (Bachu, 2003).The CO 2 plume can rise several hundreds of meters through permeable sediments in the reservoir during and shortly after the injection, leaving CO 2 bubbles trapped in rock pores by capillary forces (capillary or residual trapping) (Juanes et al., 2010;Krevor et al., 2011).The remaining free-phase CO 2 eventually becomes permanently trapped by dissolution into the resident brine (solubility trapping) and mineral precipitation (mineral trapping) at relatively slow rates (Benson & Cook, 2005;Gunter et al., 2004).The buoyant free-phase CO 2 and the associated leakage hazard may exist for hundreds of thousands of years, necessitating the presence of tight continuous caprock(s) to ensure storage integrity.Competent sealing caprocks, commonly comprising clay-rich formations (referred to here as shales), are characterized in the laboratory by ultra-low intrinsic permeability in the nanodarcy range, small porosity, and high capillary entry pressure of several MPa (Boulin et al., 2013;Espinoza & Santamarina, 2017;Heath et al., 2012;Hildenbrand et al., 2002;Makhnenko et al., 2017).These factors narrow down the sources of leakage to improperly sealed wellbores and conductive faults (Alcalde et al., 2018;Song & Zhang, 2013).
Nevertheless, the transport properties of shales are reported to be strongly scale-dependent: the permeability may increase by several orders of magnitude from the lab to the field and regional scales (Hart et al., 2006;Neuzil, 2019).The caprock sealing capacity may also gradually degrade by a number of factors: (a) overpressure creating (micro)fractures or reactivating pre-existing faults (Rutqvist, 2012;Vilarrasa et al., 2014), (b) CO 2 dissolution into water that acidifies the pore fluid and likely promotes mineral dissolution and/or transformation (Gaus et al., 2005;Rezaee et al., 2017), (c) CO 2 adsorption on clay platelets and generation of swelling pressure (Busch et al., 2016), and (d) CO 2 entry in free phase and resulting dehydration of the pore or clay interlayer giving rise to desiccation cracking (Bhuiyan et al., 2020;Espinoza & Santamarina, 2012).These mechanisms have the potential to create pervasive pore-scale heterogeneities that may act as paths of least resistance for CO 2 leakage.An example is the Sleipner storage project, where the CO 2 plume has ascended through eight thin fractured intra-formational claystone barriers (total thickness of 20 m) over the period of 3 years (Cavanagh & Haszeldine, 2014).Furthermore, even an effective capillary seal remains amenable to transport by aqueous CO 2 molecular diffusion.Indeed, diffusive transport, although inherently slow, is ubiquitous and contributes to growing uncertainties in long-term leakage evaluations (Busch et al., 2008;Kivi et al., 2022).
In this study, we address two fundamental questions: how far will the leaking CO 2 front migrate toward the surface over geological time scales, and how is the subsurface CO 2 containment affected by the caprock sealing properties?To address these questions, we substantially advance existing numerical modeling efforts that capture the dynamics of CO 2 leakage but are restricted to decadal to multi-century time scales due to high computational costs (Cavanagh & Haszeldine, 2014;Hou et al., 2012;Zahasky & Benson, 2016).To extend the time and length scales of analysis, we propose here a one-dimensional upscaling model that is arguably capable of dealing with the underpinning physics of leakage at the basin scale and over long-time periods relevant to gigatonne-scale CCS (thousands of years).We estimate the temporal evolution of CO 2 mass distribution across geological strata, finding that the possibility of CO 2 leakage through a multi-layered system of aquifers and caprocks is significantly lower than when a single, even much thicker caprock exists.Such a multi-layered geological setting builds confidence in permanent geologic storage of CO 2 at the gigatonne scale.

Description of the CO 2 Leakage Model
To simulate long-term CO 2 storage, we consider continuous CO 2 injection for 30 years at the regional scale through a dense grid of 10-km-spaced wellbores into a deep homogeneous saline aquifer (Figure 1).The planar symmetry of the problem in the x-and y-directions imposes quasi no-flow boundaries across injection sites, controlling a laterally extensive gravity override CO 2 plume and brine pressure evolutions (Bandilla & Celia, 2017;De Simone et al., 2019;Juanes et al., 2010;Nordbotten et al., 2005;Vilarrasa et al., 2010).Meanwhile, the interfering pressure perturbations of neighboring sites build up a large pressurized region touching the basin margins (Zhou et al., 2010).After the injection ceases, the plume gradually redistributes mainly due to buoyancy and the lateral expansion tapers off, consistent with short-term post-injection observations at the resolution of seismic data at the demonstration projects in Ketzin, Germany (Lüth et al., 2017) and at the CO2CRC Otway site in Australia (Dance et al., 2019).The CO 2 plume shape and pressure distributions in the storage aquifer are expected to approach quasi-equilibrium states a few tens up to hundreds of years after the injection stops (Pawar et al., 2020;Zapata et al., 2020).Injection-induced and post-injection hydraulic perturbations can be considered as instantaneous in geologic time, and one may assume that the long-term CO 2 migration will occur almost uniformly in the z-direction.This physically justifiable simplification allows use of a one-dimensional transport model to estimate an upper limit on the CO 2 leakage possibility as the plume can migrate in the post-injection period under natural groundwater flow and the reservoir slope effects (MacMinn et al., 2011).The model geometry comprises an idealized multi-layered system of aquifers and caprocks overlying the low-permeability crystalline basement (Figure 1b).

Model Parameters
The transport properties of aquifers and caprocks are given in Table S1 in Supporting Information S1.For simplicity, we assign to all aquifers the same set of properties equivalent to those of Berea sandstone as a typical reservoir formation (Vilarrasa et al., 2016).We stipulate two simulation scenarios where all caprocks share either a high or low sealing capacity representing the best-case and worst-case leakage scenarios, respectively.The former scenario adopts the sealing properties of an intact clay-rich shale caprock evaluated in the lab with a typically high-capillary entry pressure p 0 of 2.5 MPa and low intrinsic permeability k of 10 −20 m 2 (corresponding to 10 nanodarcy) (Hildenbrand et al., 2002;Makhnenko et al., 2017).The latter case considers k = 10 −16 m 2 (corresponding to 0.1 millidarcy), representative of a degraded regional-scale shale permeability (Hart et al., 2006;Neuzil, 2019), and p 0 = 0.1 MPa, comparable to that inferred from the Sleipner CO 2 plume migration (Cavanagh & Haszeldine, 2014).Another uncertain parameter is the relative permeability to CO 2 (Kivi et al., 2021), which is hypothesized here to be a power-law function of saturation with an exponent n = 6 for the high-quality seal and n = 3 for the degraded one.We assume that the diffusive transport of aqueous CO 2 through aquifers is established at the same rate as in bulk water with a coefficient D of 2 × 10 −9 m 2 /s (Tewes & Boury, 2005), whereas it slows down through the tortuous pore network of caprocks with D = 2 × 10 −10 m 2 /s, reaching the upper limit of experimentally determined values (Busch et al., 2008;Rezaeyan et al., 2022).

Numerical Simulations
We aim at numerically simulating the fate of injected CO 2 in the long term.Therefore, we implicitly account for transient injection-induced effects during the initialization of the post-injection period.Pressure distribution is initially hydrostatic with a uniform dissolved CO 2 concentration of 2.54 × 10 −4 kg/kg of low-salinity water.The model is isothermal with a temperature of 40°C for all layers.We inject CO 2 for 30 years at a constant overpressure of 6 MPa to the upper lateral segments of the reservoir while the bottom side permeates water, resembling its displacement by the spreading CO 2 plume or potential leakage to the crystalline basement or at basin margins.We make several trials by varying the water leakage parameter to achieve a CO 2 saturation column and pressure distributions representative of real practices after attaining the aforementioned CO 2 plume stabilization.The long-term CO 2 migration analysis is subjected to constant atmospheric pressure and no-flow conditions at the top and bottom boundaries, respectively.We solve the one-dimensional transport problem for 1 million years (Vilarrasa, 2022) using the finite element code CODE_BRIGHT (Olivella et al., 1996), which was extended to model the CO 2 injection (Vilarrasa et al., 2010, see Section S1 in Supporting Information S1 for the model formulation).We assume that the long-term CO 2 migration is primarily controlled by buoyancy and dissolution into the brine, neglecting potential residual and mineral trapping contributions.The mesh consists of 7,100 quadrilateral elements with a constant thickness of 5 cm for the lowermost 200-m-thick sediments while gradually enlarging toward the surface.The solver discretizes the 10 6 -year simulation period to near 14 million unequal time steps.Sensitivity analyses show that the solution in terms of the spatial and temporal evolutions of the pore pressure and CO 2 saturation and concentration does not change with further refinement of the mesh and time steps.The total CPU time is around 140 hr on a Xeon CPU of speed 2.5 GHz with 32 GB memory.

Best-Case Assessment of the Caprock Sealing Capacity
The combined effects of injection overpressure and buoyancy arising from the CO 2 column height increase the capillary pressure at the base of the caprock to a maximum value of near 2 MPa (Figure 2a).Thus, CO 2 in free phase enters the caprock, but cannot percolate through it.CO 2 remains immobile due to negligible changes in its already very small relative permeability (Figure 2b) and becomes confined to the lowermost portion of the caprock (10 cm).On the contrary, the high relative permeability of the caprock to brine facilitates its flow ahead of CO 2 .CO 2 dissolution into brine leads to a gradual decrease in the excess pressure and, thus, the capillary pressure, initiating an imbibition path after the injection ceases.Brine is imbibed into the pores from the reservoir below until recovering the initial state of saturation.Due to these phenomena, the possibility of a rapid advective leakage of free-phase CO 2 is thoroughly ruled out, and CO 2 transport is exclusively dominated by molecular diffusion.Gravity establishes a hydrostatic brine pressure distribution after 10 5 years that significantly prolongs further CO 2 pressure drop in the reservoir and boosts the diffusive leakage.The dissolved CO 2 is able to propagate through 200 m of overlying formations after 1 million years (Figure 3a).Accordingly, the amount of free-phase CO 2 within the reservoir progressively declines, whereas the storage witnesses a growing contribution of dissolved CO 2 in overlying and underlying formations (Figure 4a).

Worst-Case Assessment of the Caprock Sealing Capacity
The low-capillary entry pressure and relatively high permeability of the caprocks permit CO 2 to pass through the whole primary caprock and penetrate the secondary one in the course of the injection.CO 2 displaces a considerable portion of the mobile wetting phase out of the pore space and makes a transition to a two-phase flow (Figures 2c and 2e).The relative permeability to CO 2 increases to values as high as 0.6, significantly promoting its advective flux (Figures 2d and 2f).In the early post-injection period, CO 2 redistributes on the imbibition path across the already invaded pore network, including the whole primary caprock and the base of the secondary one, while it continues penetrating the overlying rock layers.Interestingly, the pore fluid of the overlying aquifers and caprocks becomes saturated with dissolved CO 2 (the region with a nearly constant concentration gradient proportional to pressure changes in Figure 3b), leading to the formation of two new CO 2 accumulations.The advective CO 2 migration lasting for more than 200 years strongly dominates over diffusive transport.Nevertheless, the third caprock in the sequence cuts the percolating path, and molecular diffusion begins to govern CO 2 leakage.The diffusive CO 2 front rise is limited to 300 m after 1 million years.
A significant portion of the injected CO 2 initially exists in free phase (>50%) or dissolved in brine (>10%) outside the target aquifer (Figure 4b).The contribution of secondary traps to storage slightly decreases in the long term as more CO 2 diffuses upward.Remarkably, the CO 2 accumulation below the third caprock completely dissolves into brine after several hundred thousand years, leading to a drop in CO 2 concentration (Figure 3b).Furthermore, almost one-third of the mobile CO 2 within the reservoir becomes trapped by dissolution into underlying formations.

Discussion and Conclusion
The one-dimensional transport model developed in this study provides a high-resolution, computationally cost-effective, yet realistic assessment of CO 2 leakage from a gigatonne-scale storage scenario over geological time scales.The relatively fast advective flow may dominate CO 2 migration during and shortly after injection for caprocks with a low entry pressure and relatively high permeability.However, CO 2 does not migrate far up across the multi-layered system even if the caprocks are pervasively fractured.As the CO 2 pressure gradually drops, imbibition of brine into caprocks further hinders the CO 2 leakage.The CO 2 plume that migrates out of the target storage formation accumulates below the overlying caprocks forming secondary traps.Our observations agree well with field examinations of naturally occurring CO 2 in stacked reservoirs that have rarely been accompanied  b) CO 2 relative permeability curves with four scanning times t 1 = 30 years, t 2 = 100 years, t 3 = 1,000 years, and t 4 = 10,000 years (corresponding to colors gradually shifting from dark to light blue) for a point at the base of the primary caprock in the best-case leakage assessment scenario.Similar curves are plotted for the worst-case leakage assessment scenario at the base of the primary caprock (c and d, respectively) and the base of the secondary caprock (e and f, respectively).The observation points at the basal interfaces of the primary and secondary caprocks are shown by circles and squares, respectively, in the right-hand side of the figure .by leakage to the surface over geological time scales (Miocic et al., 2016).In the long term, the diffusive leakage will not jeopardize storage security because the aqueous CO 2 can barely diffuse to a distance of 1 m over thousands of years.Our estimates of the diffusive migration rate pose an upper bound limit to previous predictions by reactive-diffusive transport modeling of a few millimeters to a few centimeters on a 1,000-year time scale (Gherardi et al., 2007;Kampman et al., 2016).
A single caprock would overall behave in a similar manner.CO 2 may create a long percolating path into the caprock as it does not anymore displace large brine volumes of alternating aquifers when forming secondary CO 2 accumulations.If CO 2 breaks through the caprock, a trade-off will likely arise between CO 2 dissolution in resident brine and buoyant rise (Sharma & Van Gent, 2018).CO 2 density decreases as it ascends, featuring its significant expansion at depths of around 800 m, where CO 2 undergoes a phase transition to the gaseous state, favoring strong buoyant rise (Pruess, 2004).Thus, the possibility of CO 2 leakage back to the surface will be governed by the amount of leaked CO 2 from the caprock available for migration through the overlying high-permeability sediments.Taking the Span and Wagner (1996) equation of state to calculate CO 2 density and Henry's law of solubility with coefficients extracted for CO 2 from experimental data at a constant temperature of 40°C (Spycher et al., 2003), we argue that a CO 2 accumulation height of 35 m at a depth of 1,500 m is required to statically saturate the pore fluid of rock formations up to the surface (assuming a uniform porosity profile).At the same time, the self-limiting hydrodynamic mechanisms could appear to drastically attenuate leakage (Oldenburg & Unger, 2003;Pruess, 2004).Nonetheless, if injection in this type of stratigraphic setting is needed to unlock more storage capacity, a caprock of even low-sealing capacity but enough thickness should exist to ensure secure retention of CO 2 underground (see Section S2 in Supporting Information S1 for more details).In the meantime, adopting appropriate techniques to monitor the evolution of the free-phase CO 2 front in the more critical injection and early post-injection periods would be of paramount importance (Romanak et al., 2012;Vilarrasa et al., 2019;Zheng et al., 2022).The isothermal conditions considered in our analysis overlook the temperature decline corresponding to the geothermal gradient as CO 2 migrates toward the surface.The temperature drop may be more pronounced due to the Joule-Thomson cooling effect, resulting from CO 2 expansion and depressurization following its phase transition from supercritical to gas at shallow depths (Vilarrasa & Rutqvist, 2017).However, the isothermal assumption would be conservative in CO 2 leakage assessment as the temperature decrease drives an increase in CO 2 viscosity (Altunin & Sakhabetdinov, 1972) and solubility in the resident brine (Spycher et al., 2003), as well as a reduction in the diffusion coefficient of CO 2 in water (Cadogan et al., 2014), all retarding CO 2 flow and transport through overlying formations.Yet, shallow sediments accommodate more CO 2 in free phase as a result of higher CO 2 density at lower temperatures (Span & Wagner, 1996).Thermally induced changes in CO 2 density and viscosity may affect its flow and, thus, alter pressure profiles, although such thermal effects are of secondary importance to CO 2 flow owing to the intrinsically slow heat transfer compared to hydraulic diffusivity of fluids (longer characteristic time of thermal vs. hydraulic processes, Taron et al., 2009).
By neglecting residual and mineral trapping mechanisms in the aquifer alongside transport-limited chemical reactions in the caprock, our evaluations of the vertical CO 2 migration are on the conservative side from the viewpoint of leakage potential.In contrast, the presence of a preferential flow conduit like a high permeability fault or fracture zone may violate our assumption of laterally uniform CO 2 leakage.Localized flow imposes minor changes in the reservoir pressure that may sustain a long-lived and vertically extended leakage (Oldenburg & Rinaldi, 2011).However, the knowledge of flow properties of faults in clay-rich formations is yet to be developed.Insights from recent experiments at the Mont Terri underground rock laboratory (Switzerland) point to significant shale fault heterogeneity with permeability variations, depending on the fault architecture, from 10 −18 to 10 −20 m 2 (Kneuker et al., 2017;Zappone et al., 2021).These values are comparable to those of Eau Claire shale or Opalinus Clay that are considered as caprock analogs (Kim & Makhnenko, 2020).The presence of flow heterogeneities and the ductile nature of clay-rich rocks, enabling fault reactivation with slight shear-induced dilation and permeability creation, are argued to considerably retard localized CO 2 migration up along the fault (Ikari et al., 2009;Rutqvist et al., 2016;Vilarrasa & Carrera, 2015).Indeed, significant permeability enhancement is anticipated after tensile fracturing occurring at notably high overpressure (Guglielmi et al., 2021), which is commonly avoided in the course of CO 2 injection.Besides, CO 2 flow through the fault may promote mineral precipitation (Nooraiepour et al., 2018) and dissolution-enhanced fine mobilization (Noiriel et al., 2007), resulting in self-sealing.Conclusively, although key research challenges concerning focused leakage (im)possibility through fault-like features have yet to be addressed, clay-rich caprocks most likely remain prominent seals for permanent gigatonne-scale CO 2 storage.

Figure 1 .
Figure 1.Conceptual representation of basin-wide CO 2 injection (not to scale).(a) A bird's eye perspective of CO 2 injection into the reservoir through a dense pattern of wellbores (filled black circles) with a spacing of 10 km.Under buoyancy, the CO 2 plume migrates along the top of the aquifer (blue circular regions; the same color is used for illustration on the x-z plane).Assumed symmetries in x-and y-directions, regardless of the interactions between individual CO 2 plumes, simplify the post-injection leakage phenomenon to a 1D (z-direction) transport problem.(b) The model comprises a sequence of aquifers and caprocks with representative thicknesses and depths.We extend the succession to the surface, allowing us to impose atmospheric boundary conditions at the top, while the lower boundary immediately above the non-porous crystalline basement is subjected to a no-flow boundary condition.

Figure 2 .
Figure2.Two-phase flow dynamics of CO 2 leakage.(a) Capillary pressure and (b) CO 2 relative permeability curves with four scanning times t 1 = 30 years, t 2 = 100 years, t 3 = 1,000 years, and t 4 = 10,000 years (corresponding to colors gradually shifting from dark to light blue) for a point at the base of the primary caprock in the best-case leakage assessment scenario.Similar curves are plotted for the worst-case leakage assessment scenario at the base of the primary caprock (c and d, respectively) and the base of the secondary caprock (e and f, respectively).The observation points at the basal interfaces of the primary and secondary caprocks are shown by circles and squares, respectively, in the right-hand side of the figure.

Figure 3 .
Figure 3. Dissolved CO 2 concentration profiles.Temporal evolution of the CO 2 concentration along a sequence of aquifers and caprocks for the (a) best-case and (b) worst-case leakage assessment scenarios.

Figure 4 .
Figure 4.The long-term fate of the injected CO 2 underground.Temporal evolution of the CO 2 mass fractions for (a) the best-case and (b) the worst-case leakage assessment scenarios.
19448007, 2022, 24, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022GL100443by Csic Organización Central Om (Oficialia Mayor) (Urici), Wiley Online Library on [10/01/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License I.R.K. and V.V. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Program through the Starting Grant GEoREST (www.georest.eu)under Grant agreement No. 801809.I.R.K. also acknowledges support by the PCI2021-122077-2B project (www.easygeocarbon.com)funded by MCIN/ AEI/10.13039/501100011033 and the European Union NextGenerationEU/ PRTR.IDAEA-CSIC is a Centre of Excellence Severo Ochoa (Spanish Ministry of Science and Innovation, Grant CEX2018-000794-S funded by MCIN/ AEI/10.13039/501100011033). R.Y.M. acknowledges the support from US DOE through Carbon SAFE Illinois Corridor Project DE-FE0031892.Additional funding for completing this manuscript was provided to C.M.O. and J.R. by the U.S. Department of Energy under contract No. DE-AC0205CH11231 to the Lawrence Berkeley National Laboratory.