Two-Phase Flow Mechanisms Controlling CO2 Intrusion into Shaly Caprock

Geologic carbon storage in deep saline aquifers has emerged as a promising technique to mitigate climate change. CO2 is buoyant at the storage conditions and tends to float over the resident brine jeopardizing long-term containment goals. Therefore, the caprock sealing capacity is of great importance and requires detailed assessment. We perform supercritical CO2 injection experiments on shaly caprock samples (intact caprock and fault zone) under representative subsurface conditions. We numerically simulate the experiments, satisfactorily reproducing the observed evolution trends. Simulation results highlight the dynamics of CO2 flow through the specimens with implications to CO2 leakage risk assessment in field practices. The large injection-induced overpressure drives CO2 in free phase into the caprock specimens. However, the relative permeability increase following the drainage path is insufficient to provoke an effective advancement of the free-phase CO2. As a result, the bulk CO2 front becomes almost immobile. This implies that the caprock sealing capacity is unlikely to be compromised by a rapid capillary breakthrough and the injected CO2 does not penetrate deep into the caprock. In the long term, the intrinsically slow molecular diffusion appears to dominate the migration of CO2 dissolved into brine. Nonetheless, the inherently tortuous nature of shaly caprock further holds back the diffusive flow, favoring safe underground storage of CO2 over geological time scales.


Introduction
The emission of huge amounts of anthropogenic CO 2 into the atmosphere has given rise to global warming and subsequently climate change. It is suggested that large-scale carbon capture and storage (CCS) in deep geological media should be part of the solution to diminish CO 2 emissions and the associated harmful consequences (IPCC 2005;Iglauer et al. 2015;Celia 2017). Injection into reservoirs deep enough to bring CO 2 to high pressure (> 7.38 MPa) and temperature (> 31.04 °C) ensures that CO 2 will remain in its supercritical state (Bachu 2003). The advantages of doing that are the high storativity and injectivity owing to the liquid-like density and gas-like viscosity of supercritical CO 2 . Yet, CO 2 remains lighter than the reservoir brine and thus, is buoyant, tending to flow upward. Further migration of CO 2 out of the target storage formation is deemed to be restricted by a low permeability (typically below 10 −18 m 2 , or equivalently 10 −3 mD) and high CO 2 entry pressure (~ 1-10 MPa) caprock overlying the reservoir (IPCC 2005;Busch et al. 2008;Song and Zhang 2013;Espinoza and Santamarina 2017). Therefore, one of the crucial scientific and technical issues that have to be addressed for the successful deployment of CCS is assessing the sealing capacity of caprocks to provide a measure of CO 2 containment in place or, conversely, its leakage back to the surface over geological timescales (IEAGHG 2011).
CO 2 is partially miscible with water, dissolving between 1 and 4 mol% into the resident brine at storage conditions (Spycher et al. 2003). The resulting concentration gradient of CO 2 across the interface between brine-saturated caprock and the storage formation serves as a driving force for molecular diffusion into the caprock. This phenomenon is believed to provide a continuous transport of CO 2 as long as the caprock constitutes an interconnected system of brine-saturated pores (Schlömer and Krooss 1997;Busch et al. 2008). Likewise, the injection-induced overpressure drives the reservoir brine, which bears dissolved CO 2 into the caprock. Contrarily, the vertical flow of CO 2 as a free phase into the caprock is restricted by capillary forces induced by fluid-fluid and fluid-rock interfacial tensions acting across the pore system. A sketch of these simultaneous flow mechanisms in the subsurface is illustrated in Fig. 1. The lowest capillary resistance against the bulk intrusion of the non-wetting CO 2 phase is imposed at the largest pore throats and is characterized by the capillary entry pressure p 0 (Hildenbrand et al. 2002;Espinoza and Santamarina 2010;Makhnenko et al. 2017). Therefore, the caprock sealing potential may be compromised if the injection-and buoyancy-driven differential pressures between wetting and non-wetting fluids at the reservoir-caprock interface surpass the entry pressure. From this point on, the Fig. 1 Sketch of the CO 2 intrusion into low-permeability caprock leakage through the caprock may be accelerated by two-phase flow mechanisms (Gherardi et al. 2007; Amann- Hildenbrand et al. 2013). The entry pressure is hence a critical parameter diagnosing a possible transition from a diffusion-dominated to a mixed advective-diffusive CO 2 flow. This necessitates in the first step the employment of appropriate laboratoryscale techniques to characterize the intrusion into the caprock and the dynamics of CO 2 leakage across it under representative subsurface conditions.
Recently, some efforts have been devoted to the experimental assessment of the CO 2 penetration into typical sealing lithotypes, mostly including shales and mudrocks (e.g., Hildenbrand et al. 2002;Li et al. 2005;Egermann et al. 2006;Carles et al. 2010;Heath et al. 2012;Espinoza and Santamarina 2017;Makhnenko et al. 2017). The well-established mercury intrusion porosimetry (MIP) technique has been widely exploited for providing rough estimates of the CO 2 entry pressure after converting the gaseous-liquid mercury to CO 2 -brine system with appropriate interfacial properties (IEAGHG 2011). Meanwhile, direct experimental approaches involving CO 2 injection into initially brine-saturated rock specimens have also been developed. The implementation of direct methods is technically challenging and extremely time-consuming, from a few weeks to months (Li et al. 2005;Amann-Hildenbrand et al. 2013), which hinders an extensive use of direct methods. Nonetheless, these methods provide virtue of accounting for factors contributing to the sealing capacity of the caprock, including thermo-mechanical dependence of the pore network characteristics (Angeli et al. 2009;Vilarrasa and Rutqvist 2017), rock fabric anisotropy (Schlömer and Krooss 1997), thermal effects on fluid-fluid and fluid-rock interfacial properties (Iglauer et al. 2015;Chalbaud et al. 2009) and CO 2 -triggered geochemical reactions causing mineral dissolution and precipitation (Noiriel et al. 2007;Kim et al. 2018), among others. The direct method involves injection of CO 2 at excess pressure, which is slowly and gradually increased until continuous CO 2 flow is observed (Thomas et al. 1968). The breakthrough pressures obtained via this method are reproducible because further water imbibition closes the pathways created by the first breakthrough experiment in the case of poroelastic deformations. Hildenbrand et al. (2002) presented an experimental scheme, widely referred to as residual pressure technique, which has received particular attention in the assessment of caprock sealing capacity. The proposed method relies on initiating a pressure-driven CO 2 stream through a brine-saturated specimen and inspecting the evolution of upstream and downstream (also saturated with CO 2 ) pressures. Using analytical models, Hildenbrand et al. (2002) interpreted the typical response acquired for mudrock specimens to coincide with a drainage path proceeded by imbibition until the pressure difference between the two sides approached a constant threshold. The CO 2 outflow was found to stop due to the capillary obstruction of the largest interconnected pore throats extended across the specimen. The corresponding residual pressure was taken equal to the breakthrough pressure P brth , defined as the minimum capillary pressure required to drive CO 2 into the pore network and form a continuous flowpath through it. Adopting almost the same testing methodology, Egermann et al. (2006) attempted to modify interpretations by breaking down the differential pressure across the specimen into two parts: capillary pressure drop at the CO 2 -brine interface after penetration of CO 2 into the rock and pressure gradient across the brine column displacing brine out of the specimen. For the latter differential pressure component, single-phase Darcy flow of brine was assumed, while possible multi-phase flow effects and CO 2 dissolution into the brine were neglected. This approach led to the measurements of the breakthrough pressure on the drainage path and relatively shorter experimental times (on the order of hours to days). It has been shown that the application of direct and indirect methods can give very similar values for clay-rich caprock (Makhnenko et al. 2017), even 1 3 though some controversies between the two methods have also been reported (Boulin et al. 2013).
In spite of the achieved progress in testing procedures, the employed analytical models are unable to effectively capture the complexities of the CO 2 breakthrough and flow in initially brine-saturated caprock. Key issues are mainly associated with the complex interplay between different flow mechanisms through the caprock sample, and controlling parameters are yet to be investigated. Numerical simulations performed to reproduce experimental observations are perceived to be highly instructive in giving insight into the underlying processes governing the problem (Damians et al. 2020). Furthermore, numerical models can be employed to conduct mechanistic and predictive studies of CO 2 flow into/through low-permeability caprocks over longer temporal and broader spatial scales than what can be achieved in the laboratory.
This study aims at shedding light on the flow processes governing CO 2 leakage across caprock by establishing an integrated workflow that brings together experimental observations and numerical simulations. We employ the residual pressure approach under realistic subsurface conditions to evaluate CO 2 penetration and flow in two analogues clay-rich caprock samples: one representing intact caprock and the other a shear zone in the caprock. Obtained results provide a comparison between the sealing capacity of intact caprock and possible faults crossing it. Next, based on the experimental data and physically sound constitutive models, we parameterize a two-phase flow model in deformable porous media to numerically simulate the core-scale flow of CO 2 -brine in the rock specimens. Through numerical simulations, we present an improved understanding of the breakthrough process across the caprock and associated flow mechanisms. We distinguish between our direct model parameterization and calibration approaches in which numerical simulations are brought to a forced fit with experimental observations by adjusting the model parameters. We conduct sensitivity analyses on uncertainties pertinent to model parameters and the testing scheme with direct relevance to underground storage of CO 2 . Relying on the flow mechanisms inferred in the core scale, we discuss and draw some general conclusions on the potential risk of CO 2 leakage in geologic carbon storage.

Materials
Shaly facies of Opalinus Clay are chosen as a caprock representative for this study. The material is considered to be a potential host for high-level nuclear waste storage and target formation for studying the effect of in situ CO 2 injection on caprock properties (Bossart et al. 2017). Opalinus Clay cores are recovered from the Mont Terri Underground Rock Laboratory in Switzerland from a depth of 1-2 m inside a borehole drilled perpendicular to the apparent bedding planes. The cores were vacuum-sealed right after the extraction and prepared specimens had the initial degree of saturation of about 0.85-0.90 and no visible major cracks. The mineral composition of the rock is 50% clay (mainly illite, kaolinite, chlorite, and smectite), 24% calcite, 20% quartz, 5% organic matter and traces of other minerals. Containing about 50% of clay size particles and being fissile means that this material can be called "shale" in its geological sense (Weaver 1988). Ductility, low permeability ~ 10 −20 m 2 (10 −5 mD), and high CO 2 entry pressure on the order of a few MPa (Makhnenko et al. 2017) indicate that this shale satisfies the requirements for a potential caprock (IEAGHG 2011). Two types of Opalinus Clay specimens are considered: intact and remolded (i.e., reconstituted) rock. The latter is introduced to increase the permeability and compressibility of an original caprock and investigate the role of these properties on CO 2 flow.
Remolded shale specimens are prepared by crushing the intact material in a grinder and sieving the particles with a size smaller than 0.5 mm. Then, brine corresponding to approximately 1.5 times the liquid limit (or 60%) is added and mixed with the powder. The obtained material is initially consolidated for 72 h at the axial stress of 350 kPa. The obtained cylindrical specimens have a porosity of 0.33 and a degree of saturation of 0.90. Crushing of the clay preserves the flake-like structure of an intact material but significantly changes its permeability and compressibility at low effective mean stresses. Subsequently, the remolded shale is subjected to high effective mean stresses, presenting a slightly larger fraction of submicron scale pores, but possessing material properties close to the ones of the intact rock: porosity reduces to 0.15, permeability drops by two orders of magnitude (from 10 −18 to 10 −20 m 2 , or equivalently from 10 −3 to 10 −5 mD), and bulk modulus increases by two orders of magnitude from 20 MPa to 2 GPa (Makhnenko et al. 2017;Makhnenko and Podladchikov 2018).
The pore size distribution in the tested intact and remolded materials is measured by mercury intrusion porosimetry tests performed on specimens with linear dimensions not exceeding 10 mm. The specimens are freeze-dried prior to mercury injection to ensure that the in situ pore fluid would not influence the test results. As the applied mercury pressure increases (up to 400 MPa), pores with smaller throat diameters (down to 4 nm) are getting filled, and the mercury volume intruded in the specimen at each step is measured. The technique also provides accurate measurements of the interconnected pore volume and total volume of the specimens; hence, interconnected porosity can be calculated. The porosity of intact shale is measured to be 0.125 and the dominant pore throat diameter is 15 nm. The diameter of pore throats in remolded shale ranges from 10 nm to 20 microns, but 15 nm remains to be the dominant pore throat diameter because the grain crushing does not affect the pores at such a small scale. Consolidation at high effective mean stresses significantly reduces the size of large pores and moves the pore size distribution curve for remolded shale closer to the one of the intact material (Makhnenko et al. 2017). The saturation curves calculated for CO 2 intrusion into brine-saturated intact and remolded shale indicate entry pressures of approximately 4 and 2 MPa, respectively (Fig. 2). One should note that the apparent intrusion of the non-wetting fluid into shale samples at low pressures (< 0.01 MPa) is associated with the mercury volume filling surface irregularities and probably, artificial fissures in samples. The interfacial tension between supercritical CO 2 and brine is taken as 30 mN/m, and the contact angle for the shale-CO 2 -brine system is assumed to be 40° (Espinoza and Santamarina 2010).

Measurements of Mechanical and Flow Properties
The pore pressure diffusion time is proportional to the characteristic length squared and is equal to days-weeks for 100-mm-long shale specimens that can be tested in conventional triaxial compression or core flooding experiments . In order to simulate the in situ condition but keep the diffusion times significantly shorter (order of tens of minutes), 1D deformation (oedometric) tests are performed on thin specimens of remolded and intact shale. Cylindrical specimens are cored and perfectly ground from a 88 mm in diameter shale core to the height h = 10 mm and diameter d = 35 mm. The cylinders are squeezed into a thick-walled stainless steel ring and put inside the oedometric cell, which provides zero lateral strain conditions (Fig. 3). The axial stress is applied and measured by a 100-kN load cell. Assuming elastic deformation of the rock specimen, the lateral stress lat can be calculated from the generalized Hooke's law as lat = 1− ax − p + p from knowing the applied vertical stress ax , Poisson's ratio , and Biot coefficient . The axial deformation is measured by a set of vertical linear variable differential transformers (LVDTs) attached to the cell and considers the correction factor for system deformation (Makhnenko et al. 2017).
The cell is instrumented with two separate upstream pressure channels and one downstream channel. The upstream and downstream brine pressure/volume controllers (GDS Instruments, UK) are utilized to induce pore pressure inside the rock (with 32 MPa capacity and 0.01 MPa accuracy) and measure the volume of fluid that enters or leaves the specimen (with 1 mm 3 resolution). The injection of CO 2 can be realized in parallel with the brine via a separate upstream channel connected to a syringe pump (Teledyne ISCO, USA) that has a pressure capacity of 53 MPa with 0.01 MPa accuracy and allows the flow control at rates as low as 1 mm 3 /min. The upstream and downstream pressure lines are instrumented with valves that can prevent the flow and pore pressure transducers providing additional pore pressure measurements (70 MPa capacity and with ε ax is axial strain, σ ax is axial stress and σ lat lateral stress 0.07 MPa accuracy) right outside of the cell, hence very close to the specimen boundaries. Saturated filter paper, previously soaked in deaired brine, is placed on the bases of the cylindrical specimen surfaces, in contact with the upper and bottom porous stones. The porous stones (initial porosity = 32%, permeability = 10 −12 m 2 or 1000 mD, thickness = 6 mm) are kept submerged in the brine to allow full saturation. Once the specimen is placed inside the cell, the upper piston is sealed with a rubber O-ring. In this work, the imposed boundary conditions are the axial stress equal to 42 MPa, pore fluid pressure of 8 MPa, and temperature equal to 40 °C. All the equipment is placed inside a thermostat room providing constant and well-maintained (within 0.2 °C) temperature during testing and excluding undesirable CO 2 phase transition within the system.
After being installed in the cell, the shale specimens are fully saturated with the brine having the chemical composition of the formation pore fluid, containing 6.13 g/l sodium chloride (NaCl), 1.63 g/l sodium sulfate (Na 2 SO 4 ), 1.00 g/l of Ca 2 Cl 2 ·2H 2 O and MgCl 2 ·6H 2 O, and traces of other chemicals (Pearson 2002). This synthetic brine is supposed to have no chemical effect on the material if tested at in situ conditions. The observed specimen swelling at the beginning of the fluid saturation stage appears to be small compared to the deformation the specimens experienced due to the change in applied axial stress, pore pressure, or temperature. At the same time, the swelling of the specimens allows achieving perfect contact between the specimens and the walls of the cell that guarantees no fluid leakage along the sidewall. The bulk modulus of the brine is measured to be 2.0 GPa-slightly less than the one of water (Makhnenko et al. 2017). The back pressure saturation method is implemented for ten days with a graduate increase of the upstream and downstream pressures and promoting the brine flow through the specimen (ASTM D4767 2011). It allows achieving full saturation at brine pressures p brine > 8 MPa that is confirmed by measurements of identical Skempton's B coefficient values, while the effective mean stress is preserved to be constant and equal to 20 MPa. Recorded undrained pore pressure increments are corrected for the influence of "dead" volume and provide B = 0.85 for the remolded shale and B = 0.81 for the intact shale, though the latter one might be affected by the anisotropy of loading. The drained Poisson's ratio is reported to be 0.30 for the remolded shale and 0.33 for the intact material (perpendicular to the bedding planes), both measured in conventional triaxial tests. Similarly, the drained Young's moduli at the in situ conditions are measured to be 2.1 GPa for remolded shale and 2.5 GPa for the intact shale. The Biot coefficient is then calculated from the data on the drained bulk response of the shale and measurements of its unjacketed bulk modulus = 8.9 GPa: = 0.69 for the Opalinus Clay and 0.74 for the remolded shale (Makhnenko and Podladchikov 2018).
Intrinsic permeability k is determined from Darcy's law (see Eq. 2) after achieving steady-state flow conditions in the brine saturated specimens by applying the differential pressure of 1 MPa and measuring the fluid outflow rate. The brine density is measured to be 1015 kg/m 3 , and its viscosity at 40 ℃ is taken as 0.0011 Pa·s. Considering all the contributing factors, permeability is measured with an accuracy of 2.5% when the duration of steady-state flow exceeds 24 h . Consistency is ensured by obtaining constant permeability values with a standard deviation of 2.5% for several consecutive measurements at constant differential fluid pressure. Results of flow tests conducted in the oedometric cell are compared to those performed with the core flooding device at the same effective mean stress values and show good agreement for both materials: k = 3 × 10 −21 m 2 for the intact shale and k = 3 × 10 −20 m 2 for the remolded shale.

CO 2 Injection Experiments
After shale's permeability is measured using brine as the only pore fluid, supercritical CO 2 (temperature = 40 ℃ and pressure above 8 MPa) is injected into the specimen. It needs to be noted that for high-permeability materials, such as sandstones and limestones, CO 2 can flow as a non-wetting fluid without fully dissolving in the brine that initially fills the pores (see e.g., Kim et al. 2018;Tarokh et al. 2020). However, in low-permeability rock-it first takes a few days (depending on the specimen size) for injected CO 2 to fully dissolve in the brine and only after that a free-phase flow can be observed.
Residual pressure technique is adapted for the breakthrough pressure measurements (Hildenbrand et al. 2002;Makhnenko et al. 2017). The axial load is preserved to be constant at 42 MPa throughout the test. The downstream pressure/volume controller (Fig. 3) is set to hold the brine pressure at a constant value of 8 MPa. The valve to the upstream line is initially closed and the upstream CO 2 controller is maintaining a higher pressure than the downstream brine pressure plus the estimated value of supercritical CO 2 breakthrough pressure (above 15 MPa), such that a continuous CO 2 flow through the rock can be promoted. Then, the valve to the upstream controller is closed, but the upstream fluid line is being opened to the specimen and CO 2 pressure starts decreasing as the flow goes on. At this time, the downstream controller is collecting the fluid. As soon as the differential pressure reaches the breakthrough pressure, the reading of the upstream pressure transducer becomes constant and the downstream pressure controller indicates approximately constant volume confirming that the fluid flow is blocked by the capillary forces. The measured residual value of the differential pressure (CO 2 pressure after equilibration minus downstream brine pressure) is considered as the supercritical CO 2 breakthrough pressure and is reported with 5% accuracy. Following CO 2 breakthrough in the shale, its viscous flow through the material is assumed. Finally, this process is followed by the imbibition experiments that allow re-saturating specimens with brine and measuring the replaced CO 2 volume-hence, the non-wetting fluid saturation after each injection stage.

Two-Phase Flow in a Deformable Porous Medium
Injection of CO 2 into a brine-saturated specimen can be approached as a two-phase flow problem. We deal with this problem in a coupled hydromechanical context to avoid introducing uncertainties to our study stemming from rock deformation even of small magnitude. To explain the coupling, the mass and momentum conservation equations for each phase (wetting and non-wetting) or each component (brine and CO 2 ) and the momentum balance equation of the porous medium have to be solved simultaneously. A detailed derivation of the governing equations has been explained by Olivella et al. (1994), and only a concise description is presented in this section.
To account for the molecular diffusive flux of CO 2 in brine, the approach of balancing the mass of each component i , based on pertinent contributions in each phase, is adopted to write the fluid mass conservation equations as: and ML −3 are the saturation and density of -phase, respectively, and is either aqueous phase a or CO 2 phase c . i represents mass fraction of component i in each phase, f i ML −3 T −1 is an external mass supply of the respective component, and i is either CO 2 or brine component. In this formulation, the mass transport of components is expressed as a combination of their molecular diffusive fluxes i i ML −2 T −1 in each phase, and the phase advective fluxes q L 3 L −2 T −1 .
We assume that Darcy's law is valid as a macroscopic momentum equation for each fluid phase , and reads as in which k L 2 is the intrinsic permeability tensor, k r [−] is the -phase relative permeability, ML −1 T −1 is the -phase viscosity, g LT −2 is the gravity acceleration, and z[L] accounts for the vertical position.
The molecular diffusive flux of dissolved species relates to their concentration in each phase through Fick's law where D i L 2 T −1 is the molecular diffusion coefficient of component i in phase , is the pore structure tortuosity, and I[−] is the identity matrix. The product of pore tortuosity and phase diffusion coefficient is defined here as the effective diffusion coefficient D.
In the above formulations, densities and viscosities, specifically those of the CO 2 phase, are highly dependent on temperature, pressure and concentrations of the coexisting components in each phase (Spycher et al. 2003). Diffusion coefficient values are, in general, governed by the state of temperature and pressure (Bear 2018). Relative permeabilities of fluid phases vary with the respective degree of saturation, which by itself is a function of the capillary pressure. The interdependence of different variables imposes high nonlinearities in the flow model that have to be accurately dealt with during model parameterization and numerical solution.
The linear momentum conservation of the porous medium, in the absence of inertial effects, is reduced to the stress equilibrium where ML −1 T −2 is the stress tensor and b ML −2 T −2 is the vector of body forces applied per unit volume of the porous medium as a whole.
Assuming linear elastic and isotropic deformation of porous media (Cheng 2016), the constitutive relationship between the total stress, strain, and pore pressure is established as: where LL −1 is the strain tensor, v = 1 + 2 + 3 is the volume strain, p = max p a , p c is the fluid pressure, E ML −1 T −2 is the Young's modulus, [−] is Poisson's ratio, and [−] is Biot effective stress coefficient. Here, the sign convention of geomechanics, i.e., positive for compressional stresses and compactive strains, is adopted. In the formulation, described through Eqs. (1)-(5), displacements and pressure of each phase are the primary unknowns that are solved numerically for the given initial and boundary conditions.

Model Setup and Numerical Solution
To better understand and quantify the governing mechanisms of supercritical CO 2 intrusion and flow in low-permeability shaly caprocks, we simulate the breakthrough experiments carried out on two representative shale specimens, remolded and intact Opalinus Clay, in an oedometric cell. Figure 3 schematically represents the geometry of the testing setup and boundary conditions in a two-dimensional axisymmetric model with the symmetry axis parallel to the specimen axis. The numerical domain contains a homogeneous shale specimen surrounded by stiff porous stones, which, in turn, are overlain and underlain by loading platens that represent the dead fluid volumes of the experimental apparatus. To mimic the experiments, the model is initialized by exerting constant axial stress of 42 MPa at the top boundary, while other boundaries are subjected to the constraint of no perpendicular displacement. Brine initially occupies the whole pore system at a pressure of 8 MPa. The initial porosity 0 is set to 0.15 and 0.125 for remolded and intact shale specimens, respectively, according to the experimental measurements under equivalent effective mean stresses.
We simulate CO 2 injection at the bottom boundary in a pressure-controlled mode through four numerical steps: (1) CO 2 pressure increase to 8 MPa, (2) the pressure is preserved constant for a time period long enough to achieve an almost uniform diffusion-driven distribution of dissolved CO 2 concentration across the model, (3) CO 2 pressure enhancement to 8.5 MPa to drain the upstream compartment, and (4) increasing CO 2 pressure to the desired value to induce CO 2 flow across the specimen. During CO 2 injection upstream (steps 1-3), the hydraulic boundary conditions account for no-flow on lateral boundaries and constant brine pressure of 8 MPa on the top and bottom boundaries. As soon as injection stops (onset of step 4), the simulation proceeds by closing the upstream to any flow while monitoring the temporal and spatial evolutions of the pressure-perturbed domain. The prescribed mechanical boundary conditions, as well as a constant temperature of 40 °C over the whole model, are maintained throughout the simulation to make sure that CO 2 stays in supercritical conditions.
The relationship between the capillary pressure p cap and degree of saturation is assumed to follow the van Genuchten (1980) model where p 0 ML −1 T −2 and m[−] are the model parameters. The former parameter is commonly rated as the gas entry pressure at a specified pressure and temperature. S ea , S ra and S maxa are the effective, residual, and maximum saturation degrees of the aqueous phase, respectively. We fit the van Genuchten model to the experimental capillary pressure curves for the CO 2 -brine system (illustrated in Fig. 2). Assuming that the residual CO 2 saturation after imbibition approaches zero (or equivalently S maxa = 1 ), we use the variables S ra , p 0 and m as fitting parameters. Best-fit van Genuchten parameters for both shale samples are provided in Table 1. The fit yields S ra for the tested specimens to be around 0.2, which is substantially lower than the values reported for some shales (e.g., Bennion and Bachu 2008). This discrepancy is important because it may affect CO 2 mobility within the specimen and, thus, our simulation results. To check whether this low value controls our system, we perform the fitting of the capillary pressure properties by fixing S ra to 0.2, 0.3, 0.4 or 0.5 and adjusting p 0 and m to achieve the best fit. The fitting results show that the CO 2 entry pressure is quite insensitive to the residual brine saturation, having values around 2 MPa and 4 MPa for the remolded and intact shales, respectively, while m goes as high as 0.9 at Nonetheless, inferred from a sensitivity analysis, increased residual brine saturations do not produce a perceptible impact on replicating CO 2 injection experiments, nor do they alter dominant flow mechanisms. Therefore, we build our reference model based on capillary pressure data included in Table 1.
The intrinsic permeability of rock may change as a result of deformation (Rutqvist and Tsang 2002). To capture these indirect hydromechanical coupling effects, we model variations of the intrinsic permeability k as a power-law function of porosity where the porosity can change due to the applied effective mean stress (Eq. 5; see also Cheng 2016 for more details). k 0 is the intrinsic permeability at initial stress and pressure conditions (i.e., prior to CO 2 injection), which was evaluated to be 3 × 10 −20 m 2 (3 × 10 −5 mD) and 3 × 10 −21 m 2 (3 × 10 −6 mD) for the remolded and intact samples, respectively. The exponent n is taken equal to 16 for the intact Opalinus Clay and 13 for the remolded sample, as derived experimentally by Kim and Makhnenko (2020). Therefore, Eq. (7) provides a measure of the intrinsic permeability changes in the examined specimens due to the deformation induced by CO 2 injection changing the pore pressure and the effective mean stress.
Relative permeabilities to both CO 2 and aqueous phases are supposed to be governed by power-law functions of the corresponding effective saturations. A power of six is postulated to replicate the two-phase flow of brine and supercritical CO 2 in low-permeability caprocks (Bennion and Bachu 2008;Vilarrasa et al. 2014Vilarrasa et al. , 2015 and is picked out for the remolded shale. To improve the numerical fit on experimental data, we use a relative permeability exponent of 2 for the intact specimen. Nevertheless, the appropriateness of our choice of relative permeability functions for the specimens still remains uncertain. We further deliberate the associated uncertainties in the Discussion section. The porous stone is assumed to have a linear relationship between the relative permeability and the degree of saturation. On the other hand, neglecting the tortuosity of the pore network, an effective diffusion coefficient of 1.6 × 10 −9 m 2 s −1 is set to deliver the diffusive flow of dissolved CO 2 species. This value lies within the experimentally measured bounds of the CO 2 diffusion coefficient in pure water under representative pressures and temperature (Tewes and Boury 2005). The parameter values for Young's modulus, Poisson's ratio, and Biot effective stress coefficient are also extracted from experiments conducted on the same materials under similar effective mean stress conditions (Makhnenko et al. 2017;Makhnenko and Podladchikov 2018). Table 1 includes the hydromechanical parameters of the shale specimens and porous stones. Except for the relative permeability curves and the molecular diffusion coefficient, whose effects on numerical simulations will be discussed through sensitivity analyses, other properties have been directly measured in the laboratory tests reproduced in this study.
In addition to our two reference numerical models with the associated parameters presented in Table 1, we numerically inspect the sensitivity of CO 2 penetration and flow behavior in the remolded shale specimen to three parameters: sample length, entry pressure of porous disks, and effective diffusion coefficient. The subsurface representativeness of the experiments may be affected by scale dependency and capillary properties of the porous disks. The scale dependency is explored by increasing the specimen length by one order of magnitude (i.e., L = 100 mm) while keeping other parameters constant. Conducting the breakthrough experiment on such a long specimen is time-consuming. Similarly, to examine the sensitivity of the system response to the capillary entry pressure of porous disks, two cases are run with higher entry pressures of 0.1 and 0.2 MPa. Eventually, varying the diffusion coefficient provides an insight into the role of the diffusive flow in the system behavior. Taking the pore structure tortuosity equal to unity, an effective diffusion coefficient identical to the aqueous diffusion coefficient (i.e., D = 1.6 × 10 −9 m 2 /s) was adopted for simulating the reference case. However, the diffusive transport of CO 2 in low-permeability mudstones is likely to be disturbed by the tortuous character of the nanometer-sized pore throats, chemical CO 2 -brine-rock interactions, and physical sorption of CO 2 , mainly on clay minerals (Busch et al. 2008), whose combined effects are partly reflected in reduced effective diffusion coefficients. Therefore, we also simulate a case with a lower diffusion coefficient of 1.6 × 10 −10 m 2 /s, slowing down the molecular diffusion of CO 2 by one order of magnitude with respect to the reference case. The selected value is considered rational for tight rock under subsurface conditions (c.f., Song and Zhang 2013).
This hydromechanical problem is simulated using the fully coupled finite element code CODE_BRIGHT (Olivella et al. 1994(Olivella et al. , 1996, which solves non-isothermal two-phase flow problems in deformable porous media. Concerning fluid properties and mutual solubilities of CO 2 and water, Vilarrasa et al. (2010Vilarrasa et al. ( , 2013 adapted the code to simulate CO 2 injection in deep geological formations. The CO 2 phase density is computed using a modified Redlich-Kwong equation of state with parameters obtained by Spycher et al. (2003). The CO 2 viscosity is calculated from the empirical correlation proposed by Altunin and Sakhabetdinov (1972). CO 2 dissolution in the aqueous phase is linearly correlated with the CO 2 partial pressure according to Henry's law. However, experimental CO 2 solubility data are better-fitted by bilinear functions of pressure (Spycher et al. 2003). Therefore, we estimate Henry's law constants separately at the temperature and pressure ranges of interest in each experiment. The fitted Henry's constant (in the mole fraction scale) are 3.6 GPa and 3.3 GPa, respectively, for CO 2 injections into the remolded and intact shale specimens.
We discretize the whole model with a structured triangular mesh comprising 6800 elements and linear interpolation functions for both displacements and pressures. All elements share laterally the same edge length of 0.5 mm. Vertically, we discretize the upstream and downstream compartments with 0.35 mm thick elements. The element size progressively decreases in the specimen from 0.35 mm at the middle part to 0.1 mm at the contacts with porous stones. This mesh structure is used for all simulations conducted in this study. As a result, more than 20,000 elements are used to simulate CO 2 intrusion and flow across the 100-mm-long specimen. This fine mesh allows for resolving the two-phase flow and rock deformation to greater accuracy at the lower specimen boundary where CO 2 injection causes a high pressure gradient.  compartment causes an increase in brine pressure within the specimen. Yet, given the difference between the CO 2 injection pressure and the initial brine pressure, i.e., 6.5 MPa, the capillary pressure exceeds the entry pressure (recall Table 1), and CO 2 enters into the rock. As CO 2 penetrates into the pores, the upstream pressure continuously declines. Additionally, the increase in brine pressure in the upstream induces brine flow through the specimen that also appears at the downstream outlet after a short period of time. The outflow volume during the first 3 h is only constituted of brine. As the free-phase CO 2 reaches the downstream side, the outflow fluid is enriched by CO 2 . Comparing the experimental data (discrete points) with simulation results (continuous curves) discloses a general consistency in terms of the described evolution trends and observations even though a perfect quantitative match is missing (Fig. 4a). In particular, the simulation yields a slower pressure decline in the upstream side, which is directly linked to the slightly lower outflow at early times. More importantly, the numerical simulation of the outflow volume deviates from the experimental trend after 10 h. Experimental measurements may suggest that outflow stops, but the terminating fall-and-rise flow behavior is correlated with the observed temperature variations and slight changes in axial load (Fig. S1). The experiment was terminated after about 14 h. However, simulation results show that the long-term upstream pressure and the cumulative outflow volume evolve continuously (Fig. 4b).

CO 2 Intrusion and Flow in the Remolded Shale-Reference Case
To gain a fundamental insight into the governing flow mechanisms, we analyze the distributions of the capillary pressure and brine saturation (Fig. 5) and the individual components of CO 2 mass flux (Fig. 6), all derived numerically at several discrete times. CO 2 injection at the upstream compartment increases the brine pressure up to 9.8 MPa (equivalent to a brine pressure buildup of 1.8 MPa) at the specimen interface right after pressurizing CO 2 to the target value (t = 0). As the specimen and the downstream compartment are initially fully saturated with brine, its relative permeability is equal to unity, and fluid flow occurs exclusively by single-phase brine advection. The capillary pressure at the upstream,  (Fig. 5a), exceeds the entry pressure of the remolded shale, and thus, CO 2 penetrates into the specimen displacing the brine and making a transition to twophase flow conditions in the upstream portion. As brine initially flows toward the downstream in advance of CO 2 , it experiences a relatively higher pressure drop at the upstream that builds up capillary pressure and promotes the desaturation of the upstream portion of the specimen for a short period of 5 h. Thereafter, the capillary pressure at the upstream decreases and the pores are getting gradually refilled with brine. Meanwhile, CO 2 dissolution into the brine creates a concentration gradient across the specimen, establishing a diffusive transport of CO 2 through the brine-filled pore network. Diffusive flux curves show peaks, corresponding to the maximum local CO 2 concentration gradient within the sample. These diffusive peaks progressively penetrate into the specimen and attenuate with time as CO 2 concentration across the specimen becomes more uniform (Fig. 6a). An unambiguous distinction between flow mechanisms in shale is achieved by comparing the distribution patterns of CO 2 mass fluxes: the diffusive flow of CO 2 in brine (Fig. 6a), CO 2 advection (Fig. 6b), and CO 2 transported dissolved into the advected brine (Fig. 6c).
Simulation results show that the volumetric gas flow (advection) driven by pressure perturbations develops only in the lowermost two millimeters of the 10-mm-thick sample next to the upstream boundary. This observation can be attributed to the low intrinsic permeability of the rock and small relative permeability to CO 2 , making CO 2 almost immobile at brine saturations greater than 0.85 (based on the power-law model considered in this study). The advected brine conveys additional amounts of dissolved CO 2 across the specimen and the downstream compartment. However, the contribution of brine advection to CO 2 mass transport is small and attenuates as the injection-induced brine overpressure decays with time. The bulk advective flux of CO 2 also sharply decreases with distance from the upstream compartment, so that the molecular diffusion stands merely as the dominating flow process. As the CO 2 front advances with time, fluid-connected pores successively Numerically simulated profiles along the remolded shale specimen and platens of individual CO 2 flux components: a molecular diffusion, b advection, and c CO 2 carried dissolved into brine at several times of the breakthrough test experience a gradual increase in the CO 2 pressure (drainage path) followed by a pressure drop (spontaneous imbibition, see Fig. 5b) due to the decline in the upstream pressure. The resaturation of the pore pressure system results in a progressive decrease of the relative permeability to CO 2 . Consequently, the pressure-driven volumetric flux dissipates with time, and flow becomes almost purely diffusion-controlled.
The processes occurring at the downstream compartment and the outlet can also be inferred from Fig. 6. After a short period of time (< 2 h), free-phase CO 2 forms at the downstream side and begin to desaturate the porous disk once its capillary entry pressure is exceeded. It takes 3 h for the supercritical CO 2 to create a path through the downstream compartment to the outlet. At about 82 h, coinciding with the development of a linear distribution of the CO 2 pressure along the entire specimen, the downstream face also stops draining and starts uptaking the wetting phase back into the pores from the adjacent porous disk. From this moment on, the whole specimen undergoes imbibition, which is maintained until recovering the initial state of complete saturation with the aqueous phase. This transition can also be assessed from the corresponding reduction in the diffusive CO 2 mass flux at the downstream face, which converts into an advection-dominated flow in the porous disk owing to its high permeability (compare Fig. 6a, b at the interface between the specimen and the porous disk at the downstream). The downstream CO 2 discharge rate reaches peak value at the critical time of 82 h and then decreases. Note that CO 2 mass transport along the downstream compartment is also under the control of brine advection, whose contribution intensifies toward the outlet as the compartment drainage accompanies the brine stream.
CO 2 penetration into the specimen, which is accompanied by changes in fluid pressures, induces hydromechanical coupling effects and perturbations in the state of effective stress. Given that the rock deformation is fully constrained in the lateral direction, the rock expands in the axial direction, i.e., vertically, and experiences a porosity enhancement as the pore pressure increases. This porosity enhancement increases the intrinsic permeability by a factor of up to 1.2 at upstream (according to Eq. 7), which is found to cause only minor effects on the pressure and flow in the system. The hydromechanical effects attenuate upward along the specimen as the overpressure diminishes gradually (Kivi et al. 2021a). Figure 7 illustrates the temporal evolution of the CO 2 pressure at the upstream side, along with the cumulative outflow volume and brine pressure at the downstream side during the experiment on intact Opalinus Clay. In response to CO 2 penetration into the specimen, the upstream pressure continuously declines at a decreasing rate with time. The high differential CO 2 pressure of 14.5 MPa established across the intact shale induces higher downstream discharge volume than that of the remolded shale, with an excess of more than 1500 mm 3 after 2000 h (compare Figs. 4b, 7b). Overall, simulation results well match experimental observations in terms of the upstream pressure, but only qualitatively in terms of the measured outflow volume (Fig. 7a).

CO 2 Intrusion and Flow in Intact Opalinus Clay
Simulation results are basically similar to those described for the remolded shale. Distributions of capillary pressure, brine saturation, and CO 2 mass flux terms along the intact Opalinus Clay at several times are depicted in Figs. S2 and S3 in the Supporting Information. These distributions reveal that injecting CO 2 at the upstream side generates a brine pressure buildup of nearly 8.7 MPa when exposed to high-pressure CO 2 (t = 0) due to flow resistance imposed by the ultra-low permeability of the Opalinus Clay. After ceasing the injection, the brine and CO 2 flow and the respective pressure evolutions at the upstream remain correlated. Because the specimen is initially fully saturated with brine and provides higher conductivity to the aqueous phase, brine moves upward through it much easier than CO 2 , and the overpressure begins to dissipate. The decrease in brine pressure accelerates the CO 2 pressure decline at the upstream side during early times (Fig. 7a). However, the capillary pressure increases at the upstream side during the first 30 h of the experiment because the brine pressure undergoes a faster drop than that of the CO 2 (Fig. S2a). As a consequence, the desaturation of the upstream portion of the specimen continues for more than 30 h, during which the brine saturation decreases down to 0.71 (Fig. S2b). This progressive decrease at the beginning of the experiment reveals a much more pronounced effect of brine pressure buildup on CO 2 flow through intact shale compared to the remolded shale. As shown in Fig. 5b, the lowest brine saturation of approximately 0.7 was achieved for the remolded specimen only after 5 h, even though CO 2 was injected at a pressure of 8.0 MPa lower than that in the intact Opalinus Clay experiment. Furthermore, the uniform distribution of the pressure gradient across the specimen is reached earlier for intact shale. Precisely, the drainage path along the entire specimen is superseded by the imbibition of brine from the downstream back into the pores after 55 h.
The choice of a lower relative permeability exponent for intact Opalinus Clay leads to differences in terms of flow mechanisms with the remolded shale in which molecular diffusion thoroughly overwhelms the advective flux. Advection significantly contributes to the advancement of CO 2 deep into the intact specimen in the short term (Fig. S3b). Nevertheless, once the brine saturation begins to build up within the specimen following the imbibition path (t > 55 h), molecular diffusion of aqueous CO 2 gradually dominates the flow (Fig. S3a) and the slope of the outflow curve substantially decreases with time (Fig. 7b). Regarding even the short-term response, the advective flow is not able to break through the 1-cm-thick specimen and flow becomes diffusion-controlled in the downstream portion of the rock. Additionally, the advected brine carries CO 2 toward the downstream in dissolution and contributes to CO 2 advancement, especially at early times during which a substantial brine pressure gradient exists across the specimen (Fig. S3c).

Sensitivity Analyses
We present here results for sensitivity analyses on three parameters: specimen length, entry pressure of porous disks, and effective diffusion coefficient. All these parameters were found to impose notable impacts on the system response, but at the same time, their effect was very similar for the intact and remolded shale. Therefore, we only illustrate the results for the remolded shale and concentrate mainly on the pressure variations at the two sides of the specimen and the cumulative outflow volume.

Specimen Length
We show the simulation results for a specimen whose length is one order of magnitude larger than the base case, i.e., 100 mm. Injecting CO 2 at the upstream compartment at the same pressure of 14.5 MPa (recall Sect. 4.1) generates a moderately higher maximum brine pressure buildup of 2.5 MPa at the upstream in the longer specimen, in comparison with 1.8 MPa for the reference case. Given that the pressure and concentration gradients across the specimen are inversely proportional to its length, the flow is significantly affected in the longer specimen: The upstream portion desaturates with a minimum liquid saturation of 0.8, which is yet insufficient to form an effectively mobilized column of free-phase CO 2 through the pore network. Furthermore, CO 2 flow in the specimen remains primarily diffusion-dominated, occurring at remarkably lower rates. The diminished mass mobilities evidently delay the brine discharge at the downstream side for more than 3 h, while it takes CO 2 more than 100 h to reach the downstream side and about 3300 h to attain the maximum discharge rate. For the 100-mm-long specimen, the fluid discharge at the downstream and the corresponding CO 2 pressure decline at the upstream significantly slow down (Fig. 8).

Entry Pressure of Porous Disks
We compare the upstream pressure decline and the downstream outflow evolution for three testing setups with capillary entry pressures (p 0 ) of 0.01 (i.e., the reference case), 0.1, and 0.2 MPa for the porous disks (Fig. 9). Both variables present strong sensitivity to the entry pressure of porous disks. Maintaining all other parameters unaltered, injecting CO 2 into the upstream compartment with a higher entry pressure of the porous disks generates a larger brine pressure buildup both in the upstream portion of the rock specimen and the adjacent porous disk, leading to a more rapid brine advection, also carrying CO 2 in solution toward the downstream. The brine migration is accompanied by a volumetric expansion of CO 2 and the associated pressure drop in the upstream compartment, as already observed for the reference case, although with modest effects. The case with the largest entry pressure, i.e., 0.2 MPa, features a rapid upstream CO 2 pressure decline down to 9.1 MPa within the first 50 h brought on by a substantial upstream brine overpressure relief. The fluid discharge at the downstream consequently undergoes two stages: one in which the outflow volume increases sharply, provoking the rapid upstream pressure decline, and another one in which the outflow volume gently increases due to the reduced pressure gradients. During both periods, the induced capillary pressure at the upstream remains insufficient to drive CO 2 as a free phase along the specimen. As a result of that, the possibility of any twophase interferences in the pore system is thoroughly ruled out and CO 2 transport across the specimen is fundamentally governed by the molecular diffusion. For the intermediate entry pressure of 0.1 MPa, the system behavior is principally subjected to similar impacts, but at smaller temporal and spatial extents. Although an injection-induced brine pressure buildup of nearly 6 MPa at the upstream is detected (compared to 1.8 MPa pressure buildup for the reference case), it dissipates much sooner due to its smaller remaining saturation in the porous disk. Therefore, the CO 2 pressure decline at the upstream and the outflow volume only slightly increase at early times. Nevertheless, when looking at the long-term response of the system, the downstream discharge volume consistently declines with increasing the entry pressure of porous disks as they become less amenable to drain. Interestingly, the upstream pressure for all scenarios approaches almost the same value after 2000 h. The

Diffusion Coefficient
We analyze two different effective diffusion coefficients: D = 1.6 × 10 −9 (i.e., the aqueous diffusion coefficient considered as the reference case) and D = 1.6 × 10 −10 m 2 /s (Fig. 10). A lower diffusion coefficient slows down CO 2 penetration into the specimen and thereby maintains a higher upstream pressure for a longer period of time, notably decreasing the outflow volume. The reduced mobility of dissolved CO 2 also prompts a modest brine pressure increase of 0.5 MPa at the upstream, compared to 1.8 MPa for the base case. The combined effects of the lower brine pressure buildup and decelerated CO 2 pressure decline give rise to additional and long-lived induced capillary pressure at the upstream side of the specimen and associated further declines in brine saturation, relative to the reference case (Fig. S4).
Compared to the reference case (Figs. 5,S3), the desaturation of the specimen takes place to a larger extent at the upstream side and lasts for a longer period. However, the advancement of the free-phase CO 2 front is hindered by the low value of CO 2 relative permeability. Moreover, the CO 2 mass transport by brine advection is considerably diminished due to the minor brine pressure buildup. Thus, diffusion dominates the other concomitant flow mechanisms across the specimen. The CO 2 outflow rate at the downstream decreases remarkably and peaks after about 440 h. The reduced CO 2 flow decelerates drainage of the downstream compartment, whose degree of saturation decreases only by 0.06 during the first 1000 h in comparison with the saturation reduction of 0.32 achieved during the same period for the reference case.

Discussion
Characterization of clay-rich geomaterials is challenging due to their low permeability, difficulties in achieving full saturation, and strong stress and temperature dependence of their properties. The experimental assessment of CO 2 intrusion and flow through the shaly caprock is exposed to similar technical challenges and is exceedingly time-consuming, particularly when a conventional stepwise pressure injection method is employed (Thomas et al.  Li et al. 2005). Alternative experiments, among them the residual pressure technique (Hildenbrand et al. 2002;Espinoza and Santamarina 2017;Makhnenko et al. 2017), have thus been proposed to significantly decrease experimental times. Compared to the standard stepwise approach measuring the capillary breakthrough pressure on the drainage path, the residual technique provides estimates of this parameter during spontaneous imbibition of the wetting phase into the specimen. In spite of potential discrepancies associated with the hysteretic drainage-imbibition cycles (Hildenbrand et al. 2002;Egermann et al. 2006;Boulin et al. 2013), the technique was shown to provide the accuracy within 5% for lowpermeability shales (Makhnenko et al. 2017). The residual pressure testing procedure is utilized in this study, where the subsurface storage conditions are reproduced with supercritical CO 2 injection experiments conducted in a high-pressure oedometric cell. By taking advantage of a numerical interpretation scheme, we gain a detailed insight into the dynamics of the CO 2 flow process in the shale specimens and the pertinent alterations of capillary pressure, CO 2 saturation, and dissolved CO 2 into the brine. The numerical simulation of the experiments also enables exploring the long-term flow and its sensitivity to a variety of parameters, which is not possible in the laboratory. Nevertheless, to prevent introducing further uncertainties into the evaluations, an accurate parameterization of a two-phase flow model for the deformable porous rock is required. The model parameterization is tackled in this study first through laboratory measurements of rock properties and then verifying the model accuracy against CO 2 injection experiments in the short term.
Our simulations satisfactorily capture the general evolution trends of the conducted experiments. The numerical characterization approach allows identifying successive flow development stages, starting from CO 2 intrusion into the specimen, followed by its drainage and finally, a transition to spontaneous imbibition. This outcome is in line with the commonly outlined mechanisms (Hildenbrand et al. 2002;Egermann et al. 2006); however, neither the time at which the drainage-to-imbibition transition occurs nor the pressure difference across the specimens, i.e., the so-called breakthrough pressure, are material properties. Rather, as shown by our sensitivity analysis (Sect. 4.3), these variables depend on the specimen length, the initial pressure at the upstream compartment, capillary entry pressure of platens, and the dead volume in the downstream platen. Thus, it is essential to conduct two-phase flow simulations while considering relevant degrees of complexities to ensure a reliable assessment of the caprock sealing capacity. Analytical approaches founded purely on the basis of capillary threshold indices may be too simplistic and unrealistic.
A quantitative comparison between experimental records and numerical simulations does not yield a perfect match, in particular for the intact Opalinus Clay, where simulations underestimate the outflow volume. This inconsistency may give clues to the existence of preferential flow paths along interconnected pores with a relatively low entry pressure and relatively high intrinsic permeability. This conjecture has been corroborated by Damians et al. (2020), who numerically replicated gas flow experiments through compacted bentonite only after embedding representative flow discontinuities into the model. Similar flow discontinuity may be unfavorably encountered or unintentionally initiated in the caprock due to the overpressure-and/or thermally driven stress perturbations (Rutqvist and Tsang 2002;Vilarrasa et al. 2014;Wenning et al. 2021). Another factor added to uncertainties can be the relative permeability curves. Although the literature on intrinsic permeability measurements in shales has been eminently enriched and extended during the past two decades (Javadpour 2009;Civan 2010;Kim and Makhnenko 2020), experimental derivation of their relative permeability curves yet remains quite challenging and rare to find (some examples provided in Ghanbarian et al. 2020). In the only study to our knowledge that addresses these measurements for the CO 2 -brine system, Bennion and Bachu (2008) point to distinct relative permeability relationships for two low-permeability shales. One highlights immobilized CO 2 at brine saturations greater than 0.8, supporting the relative permeability model of our reference simulation for the remolded shale. The other recognized behavior is indicative of an almost linear enhancement of the CO 2 relative permeability to 0.2 on the drainage path toward the residual brine saturation of 0.6. Without questioning the credibility of the limited set of conducted experiments, the reported discrepancy refutes the existence of a universal power-law scaling for clay-rich caprock samples. Consequently, more promising reproduction of laboratory trends necessitates using different relative permeability exponents for the two examined shales, i.e., 6 for the remolded shale and 2 for intact Opalinus Clay. Decreasing the relative permeability exponents down to 2 provides a good simulation fit on the upstream pressure evolution, and at the same time, improves outflow predictions due to a lower resistance to advection by the increased relative permeability (more details in Kivi et al. 2021b). Partially improved by the relative permeability relationships, the system behavior displays a key and combined dependence on a set of contributing factors. The exploited approach appears to keep the obtained results within reasonable bounds of uncertainty and therefore meets the research goal as neither of the performed sensitivity analyses fundamentally change the dominant CO 2 leakage processes.
Our simulation results indicate that CO 2 mass transfer across the shale specimen is governed by three concomitant flow mechanisms: molecular diffusion in the aqueous phase, bulk volumetric advection of CO 2 , and brine advection conveying CO 2 in dissolution. While the advective CO 2 ‫‬ flow remains confined to the lowermost part of the remolded specimen and the contribution of brine ‫‬ advection to the total amount of CO 2 transported is negligible, the general behavior of the ‫‬ system is dominated by ubiquitous molecular diffusion. The immobilization of CO 2 in ‫‬ the specimen stems from the localized increase in the relative permeability to the non-wetting ‫‬ phase. Even the enhanced relative permeability (power-law model with an exponent of 2) with an injection-induced overpressure as high as 14.5 MPa is found to promote CO 2 advection exclusively in the short term and, for the most part, within the lower half of the intact Opalinus Clay, featuring a diffusion-controlled transport in the long term. This finding contrasts with the analytical interpretations by Hildenbrand et al. (2002), who denote a steep increase in the relative CO 2 permeability across the shale specimen ‫‬ following the breakthrough. The experimental evidence that CO 2 flow is established only after ‫‬ it is fully dissolved in brine further corroborates the dominance of CO 2 diffusion. Although ‫‬ actual subsurface storage sites imply complex geometries with varying degrees of ‫‬ heterogeneity (e.g., preferential high-permeability conduits), the sequence of caprock/aquifers may be idealized as one similar to the caprock ‫‬ specimen/ porous disks multilayered system considered here (Fig. 11). Therefore, the findings of this ‫‬ experimental study may have significant implications for the caprock sealing capacity at the ‫‬ field scale.
As shown in Fig. 11, if the capillary pressure at the aquifer-caprock interface ‫‬ exceeds the entry pressure of the caprock, CO 2 begins to penetrate into it and displaces brine. Both intact Opalinus Clay as a caprock representative for CO 2 storage and remolded shale, representing a possible shear zone within the caprock, exhibit the effective sealing capacities specified by their high capillary entry pressure and ultra-low permeability in the nano-Darcy range. Experimental results in the same range are documented for shaly and carbonate facies of Opalinus Clay by Minardi et al. (2021), who conducted long-term CO 2 injection experiments (for more than 100 days) and found negligible impacts of plausible chemical reactions on the intrinsic permeability of the samples. Their observations lend support to neglecting potential geochemical interaction in our simulations. In practice, notwithstanding the high caprock capillarity, the consensus is that the CO 2 bulk penetration is more likely to happen at the early stages of the injection due to the induced overpressure (Vilarrasa et al. 2014). Nevertheless, the free-phase CO 2 may appear as an inert bystander ‫‬ and may not be of a major concern, as pointed out at the lab scale by the near-zero vertical advection flux not too far ‫‬ from the injection interface. Additionally, the sensitivity analysis to the specimen length indicates that further increase of this parameter decreases the likelihood of capillary leakage. Note that the scale effect is more pronounced in the field because caprocks can have a thickness of tens to hundreds of meters, favoring safe storage of CO 2 with negligible leakage provided that no discontinuities or inhomogeneities with relatively high permeability and low entry pressure cross the caprock and excessively high injection overpressures are prevented.  Fig. 11 Conceptual representation of flow mechanisms controlling potential CO 2 leakage through the caprock: a at early times: CO 2 flow is advection-dominated if p c > p 0 , and b at later times: CO 2 flow is diffusion-dominated Not only the hydraulic properties of the caprock, but also those of the overlying and underlying aquifers do impose constraints on CO 2 leakage. Besides the intrinsic permeability of the aquifers affecting the vertical interlayer fluid migration (Birkholzer et al. 2009), their capillary characteristics can also be a contributing term, as evidenced by our sensitivity analysis to the entry pressure of the porous disks. Aquifers with relatively high entry pressure may accommodate CO 2 volumes at markedly reduced leakage risk. However, one should bear in mind that injection in such an aquifer builds up the brine pressure and does not favor high injectivities (Nordbotten et al. 2005;Vilarrasa et al. 2010). These pressure perturbations may instigate vertical or lateral brine flows and again inflict environmental concerns. On the other hand, the value of the effective diffusion coefficient D directly contributes to our ‫‬ estimates of the CO 2 leakage vertically out of the storage reservoir. Tewes and Boury (2005) reported diffusion coefficients of CO 2 in pure water ranging from 7.2 × 10 −10 to 1.8 × 10 −9 m 2 s −1 at 40 °C, increasing with pressure from 3 to 8 MPa. Laboratory experiments by Busch et al. (2008) returned an effective diffusion coefficient ‫‬ of 3.08 × 10 −11 m 2 s −1 for the Muderong Shale, whereas Wollenweber et al. (2010) reported a higher value of 7.8 × 10 −11 m 2 s −1 for the Emscher Marl. Equation (3) permits the estimation of tortuosity values of 0.014 and 0.036 for the Muderong Shale and Emscher Marl, respectively. As evidenced in the sensitivity analysis to the diffusion coefficient, the outlined tortuous character of the shaly caprock reasonably decelerates the diffusion front while retaining the flow behavior diffusion-dominated, suggesting an enhanced sealing capacity. Potential mineral reactions/adsorption processes may also play a relevant role by either increasing (self-enhancing behavior, Busch et al. 2008;Wollenweber et al. 2010) or decreasing (self-sealing behavior, Gherardi et al. 2007) the effective diffusion coefficient. Consequently, the long-term assessment of the interlayer flowpath connections necessitates the application of reactive transport modeling and careful estimates of the effective diffusion coefficient under realistic subsurface conditions.
Overall, deep sedimentary formations are targeted to effectively store CO 2 during geological time scales (i.e., orders of magnitude longer than the experiments). Yet, the question of whether the geologic CO 2 storage remains safe from the perspective of storage efficiency or may negatively affect the overlying strata due to diffusion-driven CO 2 leakage in the order of thousands of years remains open. To this end, field-scale simulations in multilayered caprock/aquifer systems over the time scale of interest are necessary but are computationally demanding, so they remain a topic of ongoing research.

Conclusions
We analyze the intrusion and flow of supercritical CO 2 in clay-rich low-permeability saturated rock, namely Opalinus Clay (representative caprock) and remolded shale (representative of a shear zone in the caprock). Numerical simulations allow capturing the laboratoryscale dynamics of CO 2 capillary intrusion and flow (i.e., a drainage-imbibition sequence) through caprock specimens of high capillary entry pressures, which are on the order of a few MPa, and ultra-low intrinsic permeability in the nD range. Simulation results provide a satisfactory replication of the experimentally observed evolution trends in terms of the CO 2 pressure and leaked volume. Thus, the parametrized model suggests a wide span of functionality, especially in addressing the caprock sealing capacity. Nevertheless, a perfect quantitative match is missing, highlighting uncertainties associated with possible preferential flow paths within the specimen and relative permeability curves. Therefore, special attention should be devoted to relevant characterization techniques.
This study reveals that the CO 2 viscous flow is principally localized to the lower part (the injection side) of the specimen, where the effective permeability to CO 2 partially increases following specimen desaturation. The CO 2 -rich brine advection also contributes to CO 2 transport, albeit at negligible extents. CO 2 flow is found to be specifically dominated by relatively slow molecular diffusion over much of the spatial and temporal domains. This finding contradicts a commonly held belief that treats CO 2 flow through the caprock as a capillary breakthrough phenomenon. Indeed, CO 2 injection has the potential to desaturate the lower zone of the caprock. However, CO 2 in free phase likely remains immobile and does not migrate deep into the caprock, which would not therefore endanger its sealing capacity. The thicker the caprock, the lower the extent of bulk CO 2 front advection and, consequently, the risk of leakage. The vertical migration of brine is far in advance of the CO 2 front. As a result, an excess storage volume is provided within the caprock, allowing for a notable CO 2 pressure drop that reduces its leakage risk. Disregarding the inherently tortuous nature of shaly caprocks, surface adsorption and chemical reactions, potentially retarding CO 2 migration, a small fraction of the stored CO 2 may scape out of the storage aquifer by diffusion at geological time scales.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s11242-022-01748-w. permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.