\jyear

2023

[3,4]\fnmMartin \surFernø

[5,4]\fnmJan M. \surNordbotten

[1,2]\fnmRuben \surJuanes

[1]\orgdivDepartment of Civil and Environmental Engineering, \orgnameMassachusetts Institute of Technology, \orgaddress\cityCambridge, \postcode02139, \stateMA, \countryUSA

2]\orgdivEarth Resources Laboratory, Department of Earth, Atmospheric and Planetary Sciences, \orgnameMassachusetts Institute of Technology, \orgaddress\cityCambridge, \postcode02139, \stateMA, \countryUSA

3]\orgdivDepartment of Physics and Technology, \orgnameUniversity of Bergen, \orgaddress\cityBergen, \postcode5020, \countryNorway

4]\orgnameNorwegian Research Center, \orgaddress\streetPostboks 22 Nygårdstangen \cityBergen, \postcode5020, \countryNorway

5]\orgdivCenter for Modeling of Coupled Subsurface Dynamics, Department of Mathematics, \orgnameUniversity of Bergen, \orgaddress\cityBergen, \postcode5020, \countryNorway

Direct Comparison of Numerical Simulations and Experiments of CO₂ Injection and Migration in Geologic Media: Value of Local Data and Predictability

\fnmLluís \surSaló-Salgado lsalo@mit.edu \fnmMalin \surHaugen Malin.Haugen@uib.no \fnmKristoffer \surEikehaug Kristoffer.Eikehaug@uib.no Martin.Ferno@uib.no Jan.Nordbotten@uib.no juanes@mit.edu * [ [ [ [

Abstract

Purpose: The accuracy and robustness of numerical models of geologic CO₂ sequestration are almost never quantified with respect to direct observations that provide a ground truth. Here, we conduct CO₂ injection experiments in meter-scale, quasi-2D tanks with porous media representing stratigraphic sections of the subsurface, and combine them with numerical simulations of those experiments.

Goals: We evaluate (1) the value of prior knowledge of the system, expressed in terms of ex-situ measurements of the tank sands’ multiphase flow properties (local data), to obtain an accurate simulation; and (2) the predictive capability of the matched numerical models, when applied to different settings.

Methods: We match three different simulation models —each with access to an increasing level of local data—to a CO₂ injection experiment in tank 1 (89.7 $\times$ 47 $\times$ 1.05 cm). Matching is based on a quantitative comparison of CO₂ migration at different times from timelapse image analysis. Next, we simulate a different injection scenario in tank 1, and, finally, a different injection scenario in tank 2 (2.86 $\times$ 1.3 $\times$ 0.019 m), which represents an altogether different stratigraphic section.

Results and conclusion: Our models can qualitatively match the CO₂ plume migration and convective mixing of the experimental truth. Quantitatively, simulations are accurate during the injection phase but their performance decreases with time. Using local data reduces the time required to history match. The predictive capability of matched models, however, is found to be similar. The sand-water-CO ${}_{2(\text{g})}$ system is very sensitive to effective permeability and capillary pressure changes; where heterogeneous structures are present, accurate deterministic estimates of CO₂ migration are unlikely.

keywords:

CO₂ storage, geologic carbon sequestration, two-phase flow, numerical simulations, history matching, FluidFlower

1 Introduction

CO₂ capture and subsequent geologic carbon sequestration (GCS) is a climate-change mitigation technology that can be deployed at scale to oﬀset anthropogenic CO₂ emissions during the energy transition (Marcucci et al., 2017; European Academies Science Advisory Council (EASAC), 2018; Celia, 2021; Intergovernmental Panel on Climate Change (IPCC), 2022). In GCS, reservoir simulation is the primary tool to assess subsurface CO₂ migration, which is necessary to understand and manage geologic hazards such as fault leakage (e.g., Caine et al., 1996; Ingram and Urai, 1999; Nordbotten and Celia, 2012; Zoback and Gorelick, 2012; Juanes et al., 2012; Jung et al., 2014; Vilarrasa and Carrera, 2015; Saló-Salgado et al., 2023) and induced seismicity (e.g., Cappa and Rutqvist, 2011; Zoback and Gorelick, 2012; Juanes et al., 2012; Ellsworth, 2013; Verdon et al., 2013; Alghannam and Juanes, 2020; Hager et al., 2021). In response to the inherent uncertainties associated with modeling and simulation of CO₂ storage (Nordbotten et al., 2012), building confidence in simulation models requires calibration (or, synonymously, history matching), a process that involves updating the reservoir model to match field observations as they become available (Oliver and Chen, 2011).

History matching is an ill-posed inverse problem (Oliver and Chen, 2011). This means that multiple solutions (i.e., parameter combinations) exist that approximate the data equally well. Automated techniques such as Markov chain Monte Carlo, randomized maximum likelihood or ensemble-based methods can be used to quantify uncertainty in history-matched models, especially in combination with surrogate models to reduce forward model computational time (see Aanonsen et al., 2009; Oliver and Chen, 2011; Jagalur-Mohan et al., 2018; Jin et al., 2019; Liu and Durlofsky, 2020; Santoso et al., 2021; Landa-Marbán et al., 2023, forthcoming, and references therein). In practice, however, it may be difficult to ensure that the chosen simulation model provides the best possible forecast. This is due to different subsurface conditions, the inability to include all sources of uncertainty in the models, incomplete field data and time limitations.

In the laboratory, intermediate-scale ( $\sim$ meter) experiments have been used to study the physics of petroleum displacement (e.g., Gaucher and Lindley, 1960; Brock and Orr, 1991; Cinar et al., 2006) and contaminant transport (e.g., Silliman and Simpson, 1987; Wood et al., 1994; Lenhard et al., 1995; Fernández-García et al., 2004). Similar 2D and 3D flow rigs have recently been applied to CO₂ storage, providing a link between core-scale measurements and field observations: Neufeld et al. (2010) studied the scaling of convective dissolution and found it to be an important mechanism in the long-term trapping of injected CO₂ in an idealized site. Wang et al. (2010) used a 3D setup to investigate the ability of electrical resistivity tomography to identify localized leaks. Trevisan et al. (2014); Trevisan et al. (2017) focused on the impact of structural and residual trapping. In homogeneous sands, they found that previous trapping models, such as the Land (1968) model, can approximate the residually trapped gas saturation (R ${}^{2}>0.6$ ). Studying an heterogeneous aquifer characterized by a log-normal distribution of six different sand facies, they report that trapping efficiency increased significantly due to structural trapping. Askar et al. (2021) used a $\sim$ 8 m-long tank to test a framework for GCS monitoring of CO₂ leakage. These studies employed homogeneous glass beads or sands, or focused on heterogeneities in the aquifer layer; structural complexity was minimal.

In this paper, we use quasi-2D, intermediate-scale experiments of CO₂ storage to evaluate, quantitatively, the predictive capability of history-matched simulation models against well-defined spatial data. An attempt was made to recreate realistic basin geometries, including stacking of storage reservoirs, faults, caprock and overburden. We simulate each of the three presented experiments with three different simulation models. This allows us to assess (1) the value of local information of the system, expressed in terms of sand petrophysical measurements, during history matching, and (2) transferability or predictive capability of our matched simulation models, when tested against a different experiment.

2 Physical Experiments

The physical experiments of CO₂ injection are conducted using the FluidFlower rigs. These rigs are meter-scale, quasi-2D tanks with transparent Plexiglass panels designed and built in-house at the University of Bergen (Fig. 1). Here, we used two tanks, with dimensions $89.9\times 47\times 1.05$ cm and $2.86\times 1.3\times 0.019$ m (referred herein to as tank 1 and tank 2, respectively) . Different geologic settings are constructed by pouring unconsolidated sands with desired grain sizes in the water-saturated rigs. The rigs have multiple ports which allow flushing out fluids after a given CO₂ injection, such that multiple injections can be conducted in the same setting. The location of the ports can be adjusted to accommodate different injection scenarios. A variety of techniques have been developed by UiB engineers in order to build complex structures such as folds and faults.

Refer to caption — Figure 1: Overview of the FluidFlower rigs and porous media used in the physical experiments. a Medium FluidFlower rig (tank 1). b Snapshot during sand pouring to build the porous medium used in experiments A1 and A2 in tank 1 (Haugen et al., 2023). c Front view of porous medium in tank 1, with lithologies in white and injector location shown with a red star. The length and height correspond to the porous medium. Note the fixed water table at the top. d Overview of the main FluidFlower rig (tank 2), showing the back panel with sensor network. e Porous medium in tank 2, used for experiment B1, with lithologies in white. Location of injectors and boxes A, B and C for analysis are shown with a red star and gray boxes, respectively. Length and variable height correspond to the porous medium.

Below, we summarize the petrophysical measurements, experimental setup, geologic model/porous media construction and experimental schedule. Details on the conceptualization of the FluidFlower rigs and technical information are given in Fernø et al. (2023); Eikehaug et al. (2023), this issue, while the full description of the physical experiment in tank 1 and ex-situ measurements are provided by Haugen et al. (2023), this issue. Further details on the experiment in tank 2 is provided by Flemisch et al. (2023), this issue.

2.1 Sand petrophysical properties

Measurements on the employed Danish quartz sands were conducted using specialized equipment to determine average grain size ( $d$ ), porosity ( $\phi$ ), permeability ( $k$ ), capillary entry pressure ( $p_{\text{e}}$ ) and drainage and imbibition saturation endpoints (denoted as connate water saturation, $S_{\text{wc}}$ , and trapped gas saturation, $S_{\text{gt}}$ ). The methodology is described by Haugen et al. (2023), this issue, and obtained values are provided in Tab. 1.

Table 1: Petrophysical properties for used quartz sands, as obtained from local, ex-situ measurements. Measured gas column heights for sands E-G were 0, so

p_{\text{e}}

could not be directly quantified. Methodology is provided by Haugen et al. (2023).

Sand type	$d$ [mm]	$\phi$ [-]	$k$ [D]	$p_{\text{e}}$ [mbar]	$S_{\text{wc}}$	$S_{\text{gt}}$
ESF	0.2	0.435	44	15	0.32	0.14
C	0.66	0.435	473	3	0.14	0.1
D	1.05	0.44	1110	1	0.12	0.08
E	1.45	0.45	2005	-	0.12	0.06
F	1.77	0.44	4259	-	0.12	0.13
G	2.51	0.45	9580	-	0.1	0.06

2.2 Experimental setup

The front and back panels of the FluidFlower are mounted on a portable aluminum frame, such that boundaries are closed on the sides and bottom (no flow). The top surface is open and in contact with fluctuating atmospheric pressure (Fig. 1). A fixed water table above the top of the porous medium was kept throughout the experiments conducted here. The experimental setup incorporates mass flow controllers to inject gaseous CO₂ at the desired rate, and a high-resolution photo camera with time-lapse function (Haugen et al., 2023, this issue).

Experiments were conducted in 2021 and 2022 in Bergen (Norway) at room temperature ( $\approx 23$ ^∘C) and fluctuating atmospheric pressure. Temperature changes were minimized as much as possible, but maintaining a constant temperature was not possible due to available laboratory space. The fluids and sands were set in the FluidFlowers using the following procedure:

1.

The silica sands are cleaned using an acid solution of water and HCl to remove carbonate impurities.
2.

The FluidFlower rig is filled with deionized water.
3.

Sands are manually poured into the rig using the open top boundary, in order to construct the desired porous medium.
4.

A pH-sensitive, deionized-water solution containing bromothymol blue, methyl red, hydroxide and sodium ions is injected through multiple ports until the rig is fully saturated. This enables direct visualization of CO₂ gas (white), dissolved CO₂ (orange to red), and pure water (blue).
5.

5.0 purity (99.999%) CO₂ is injected as gaseous phase at the desired rate.
6.

After the injection phase, injection ports are closed and CO₂ migration continues.
7.

Once the experiment is finished, the rig can be flushed with deionized water and the process can start again from step 4.

Full details on the fluids are given in Fernø et al. (2023); Eikehaug et al. (2023), this issue. Below, we refer to the pH-sensitive solution in the rigs as “dyed water”.

2.3 Porous media geometries

The geometries of the porous media used in this paper aim to recreate the trap systems observed in faulted, siliciclastic, petroleoum-bearing basins around the world, given the geometrical constraints of the FluidFlowers and manual sand pouring (Fernø et al., 2023; Eikehaug et al., 2023, this issue). Features such as folds, faults and unconformities were built in both tank 1 and 2. The construction of faults, shown in Fig. 1b and detailed in Haugen et al. (2023), requires a minimum thickness; hence, our fault structures are thicker than natural faults with the same displacement (Childs et al., 2009). Fine sands ( $d\approx 0.2$ mm) are used as sealing or caprock formation.

The geometry in tank 1 (Fig. 1c) contains three main high-permeability reservoirs (F sand). The bottom and middle F sand are separated by a seal (ESF sand), while the middle and top are separated by the C sand and connected through a higher permeability fault (refer to sect. 2.1 for pertrophysical properties). The fault separates the bottom section in two compartments. The bottom and top F sand provide anticlinal traps for the CO₂ to accumulate in.

The geometry in tank 2 (Fig. 1e) was specifically motivated by the structure of North Sea reservoirs and petroleum basins. From bottom to top, it contains two sections of decreasing-permeability reservoirs capped by two main sealing layers. A fault separates the bottom section in two compartments, while two faults separate the top section in three compartments. Each fault has different petrophysical properties: The bottom fault is a heterogeneous structure containing ESF, C, D, F and G sands, the top-left fault is an impermeable structure made of silicone and the top-right fault is a conduit structure containing G sand.

2.4 Experimental injection schedule

The injection schedule for experiments in tank 1 is provided in Tab. 2. Injection ports have an inner diameter of 1.8 mm.

Table 2: Schedules for the three CO₂ injection experiments simulated in this work.

I_{\text{R}}

is injection rate, while

I_{i}

denotes injector (port) number. A five-minute ramp-up and ramp-down was applied in experiments A1 and A2 in tank 1. Total duration of conducted experiments and simulations is 48h (A1), 5h (A2) and 120h (B1). Location of injection wells is provided in Fig. 1.

Experiment A1		A2		B1
$I_{\text{R}}$ [ml/min]	$t$ [hh:mm:ss]	$I_{\text{R}}$	$t$	$I_{\text{R}}$	$t$
0.1 ( $I_{1}$ )	00:00:00	0.1 ( $I_{1}$ )	00:00:00	10.0 ( $I_{1}$ )	00:00:00
2.0	00:05:00	2.0	00:05:00	10.0	05:00:00
2.0	00:50:00	2.0	04:43:44	0.0	05:00:01
0.0	00:55:00	0.0	04:48:33	10.0 ( $I_{2}$ )	02:15:00
0.1 ( $I_{2}$ )	01:09:11	0.0	05:00:00	10.0	05:00:00
2.0	01:14:11			0.0	05:00:01
2.0	02:29:11			0.0	120:00:00
0.0	02:34:00
0.0	48:00:00

3 Numerical simulations

3.1 Model setup

The isothermal simulations presented in this work were performed with the MATLAB Reservoir Simulation Toolbox, MRST (Krogstad et al., 2015; Lie, 2019; Lie and Møyner, 2021). Specifically, we used the black-oil module, which is based on fully implicit solvers with automatic differentiation, and assigned properties of water to the oleic phase, such that the gaseous phase (CO₂ only) can dissolve in it. Vaporization of water into the gas phase and chemical reactions are not considered, since they are not primary controls on fluid migration for our operational setup and analysis time.

In addition to structural and dissolution trapping, we also considered residual trapping (Juanes et al., 2006) to be consistent with local measurements showing nonzero trapped gas saturation (sect. 2.1). This is achieved through hysteretic relative permeability curves for the nonwetting (gas) phase (see sect. 3.2). Our implementation in MRST follows ECLIPSE’s technical description (Schlumberger, 2014), and Killough (1976)’s model is used to compute the scanning curves (Saló-Salgado et al., 2023, forthcoming). Physical diffusion was also included through the addition of a diffuse flux term with a scalar, constant coefficient in the computation of the total CO₂ flux (Bear, 1972).

The simulator requires minute timesteps due to the buoyancy of CO₂ at atmospheric conditions and high sand permeabilities (Tab. 1). Linear solver time was reduced by means of AMGCL (Demidov and Rossi, 2018; Lie, 2019), an external, pre-compiled linear solver. The greatest challenge was the convergence of the nonlinear solver, which required many iterations and timestep cuts. This is consistent with the experience reported from the groups working in the FluidFlower international benchmark study (Flemisch et al., 2023, this issue).

Next, we describe the computational grids for experiments in tanks 1 and 2, PVT properties and boundary conditions. Petrophysical properties are model-specific and are detailed in sect. 3.2.

3.1.1 Computational grids

A front panel image of the porous medium was used to obtain layer contact coordinates through a vector graphics software (Fig. 2a). These contacts were then imported into MATLAB to generate the computational grids using the UPR module (Berge et al., 2019, 2021)(Fig. 2b,d). The grids were generated in 2D and then extruded to 3D (using a single cell layer) to account for thickness and volume. Note that, in tank 1, where the porous medium has dimensions of $89.7\times 47\times 1.05$ cm, the thickness (space between the front and back panels) is constant (10.5 mm). Tank 2, which is significantly larger (porous medium dimensions $2.86\times 1.3\times 0.019$ m), has a thickness of 19 mm at the sides; however, it varies towards the middle due to forces exerted by the sand and water, to a maximum of 28 mm. A thickness map obtained after initial sand filling was used to generate our variable-thickness mesh via 2D interpolation (Fig. 2c). Also, the top surface of the porous medium is not flat (height = $130\pm 3$ cm).

Our composite Pebi grids (Heinemann et al., 1991) have a Cartesian background and are refined around face constraints (contacts and faults) as well as cell constraints (injection wells) (Berge et al., 2019, 2021). We generated multiple grids to test the finest grid we could afford to simulate experiment B1 in tank 2 with. Our grid has a cell size $h\approx 5$ mm and 151,402 cells (Fig. 2d). The grid used for tank 1 has a similar cell size ( $h\approx 4$ mm and 27,200 cells), which was chosen to reduce grid-size dependencies when applying our matched models to experiment B1.

3.1.2 PVT properties

Consistent with experimental conditions, our simulations are conducted at atmospheric conditions ( $T=25$ C), where the CO₂ is in gaseous state. We employed a thermodynamic model based on the formulations by Duan and Sun (2003) and Spycher et al. (2003); Spycher and Pruess (2005) to calculate the composition of each phase as a function of $p$ , $T$ . The implementation for a black-oil setup is described in Hassanzadeh et al. (2008) and references therein. Given the boundary conditions (sect. 3.1.3) and dimensions of our experimental porous media, pore pressure changes ( $\Delta p$ ) are very small in our simulations (max $\Delta p\ll 1$ bar). Hence, the fluid properties remain similar to surface conditions, where the water and CO₂ have, respectively, a density of 997 and 1.78 kg/m³, and a viscosity of 0.9 and 0.015 cP. The maximum concentration of CO₂ in water is $\approx 1.4$ kg/m³.

3.1.3 Initial, boundary and operational conditions

Our porous media are fully saturated in water at the beginning of CO₂ injection. No-flow boundary conditions were applied everywhere except at the top boundary, which is at constant pressure and includes a fixed water table a few cm above the top of the porous medium. Injection is carried out via wells completed in a single cell at the corresponding coordinates. The diameter of injection wells is 1.8 mm in both tank 1 and tank 2, which operate at a constant flow rate (see sect. 2). The simulation injection schedule follows the experimental protocol, provided in Tab. 2. Note that injection rates in our simulations of experiment A1 and A2 were slightly adjusted during the calibration procedure, as explained in sect. 3.3 and 4.

3.2 Simulation models

Three different models, I to III, are used throughout this study to evaluate the value of local data in predicting subsurface CO₂ migration. Each successive model was constructed based on access to an increasing level of local data, with model I having access to the least data and model III having access to the most data. The model-specific parameters are limited to the following:

•

Petrophysical properties (porosity, permeability, capillary pressure and relative permeability), which depend on available local data and are described in this section.
•

The molecular diffusion coefficient ( $D$ ). Models I-III were calibrated using the same value, $D=10^{-9}$ m²/s. Additionally, model III was also calibrated with $D=3\times 10^{-9}$ m²/s.
•

Injection rate. Experiments in tank 1 were conducted at a very low injection rate ( $I_{\text{R}}=2$ ml/min, see Tab. 2). Given that the mass flow controllers used in tank 1 may be inaccurate for this rate, the injection rate was also modeled as an uncertain parameter. Model calibration was achieved with $I_{\text{R}}\in[1.6,\,1.8]$ ml/min for all three models.

All other model characteristics, including the grid and numerical discretization, remain unchanged. Below, we describe the starting petrophysical values for each of our three simulation models. Note that the experimental geometry in tank 1, used for matching, only contained sands ESF, C, E and F. Properties for sands D and G are also provided because they were required to simulate the experiment in tank 2 (Fig. 1).

3.2.1 Model I

For this model, local petrophysical data was limited to a measure of the average grain size ( $d$ ; see sect. 2.1 and Tab. 1). Hence, petrophysical properties were estimated from published data in similar silica sands. Porosity was selected from data in Beard and Weyl (1973) and Smits et al. (2010) for moderately to well-sorted sands. Permeability was obtained from fitting a Kozeny-Carman model to data in Beard and Weyl (1973) and Trevisan et al. (2014). The resulting equation has the form $k=\alpha d^{2}\phi^{3}$ , where $\alpha$ equals 12,250.0 in our fit with $d$ in mm and $k$ in D. Obtained porosity and permeability values are provided in Table 3.

Table 3: Initial porosity and permeability for model I. See main text for estimation details.

Sand type	$d$ [mm]	$\phi$ [-]	$k$ [D]
ESF	0.2	0.37	25
C	0.66	0.38	290
D	1.05	0.40	930
E	1.45	0.39	1530
F	1.77	0.39	2280
G	2.51	0.42	5720

Capillary pressure curves were computed as described below:

1.

Capillary pressure measurements in a similar system were obtained from the literature. In this case, Plug and Bruining (2007) measured capillary pressure curves on the unconsolidated quartz sand-CO₂-distilled water system at atmospheric conditions. We used their measurements on sand packs with an average particle size between 0.36 and 0.41 mm, which are closest to the C sand in our experiments (Fig. 3a).
2.

A Brooks and Corey (1964) model of the form $p_{\text{c}}=p_{\text{e}}(S_{\text{w}}^{*})^{-\frac{1}{\lambda}}$ was fitted to these data, where $p_{\text{e}}$ is the nonwetting phase entry pressure at $S_{\text{w}}=1$ , $\lambda=2.6$ and $S_{\text{w}}^{*}=\frac{S_{\text{w}}-S_{\text{wc}}}{1-S_{\text{wc}}}$ is the normalized water saturation with irreducible or connate water saturation $S_{\text{wc}}$ . This fit led to our reference curve, $p_{\text{cr}}$ (Fig. 3a).
3.

The capillary pressure depends on the pore structure of each material, such that sands with different grain sizes require different $p_{\text{c}}$ curves. The capillary pressure variation can be modeled by means of the dimensionless $J$ -function proposed by Leverett (Leverett, 1941; Saadatpoor et al., 2010): $J(S_{\text{w}})=\frac{p_{\text{c}}}{\sigma\cos\theta}\sqrt{\frac{k}{\phi}}$ , where $\sigma$ is the surface tension and $\theta$ the contact angle. Assuming the same wettability and surface tension for different sand regions, and the same shape of the $p_{\text{c}}$ curve, the capillary pressure for any given sand ( $p_{\text{cs}}$ ) can be obtained from the reference curve as $p_{\text{cs}}(S_{\text{w}})=p_{\text{cr}}(S_{\text{w}})\sqrt{\frac{k_{\text{r}}\phi_{\text{s}}}{k_{\text{s}}\phi_{\text{r}}}}$ (Fig. 3b).

Drainage relative permeabilities were obtained from CO₂-water measurements by DiCarlo et al. (2000), who used water-wet sandpacks with 0.25 mm grain size. Specifically, we used the data reported in their Fig. 4 and 5, and fitted Corey-type functions (Corey, 1954; Brooks and Corey, 1964) of the form $k_{\text{rw}}=(S_{\text{w}}^{*})^{a}$ and $k_{\text{rg}}=c(1-S_{\text{w}}^{*})^{b}$ (Fig. 3c). The fitted exponents $a$ and $b$ are 4.2 and 1.4, respectively, while $c$ is 0.97. We assumed that the difference in relative permeability of different sands is the result of different irreducible water saturation only (see Fig. 3d,e). For each of our sands, $S_{\text{wc}}$ was obtained from Timur (1968) as $S_{\text{wc}}=0.01\times 3.5\frac{\phi^{1.26}}{k^{0.35}}-1$ , where $\phi$ is in percent and $k$ in mD. This model was used to compute $S_{\text{wc}}$ for both the $p_{\text{c}}$ and $k_{\text{r}}$ curves.

In CO₂ storage, secondary imbibition occurs where the water displaces buoyant gas at the trailing edge of the CO₂ plume, disconnecting part of the CO₂ body into blobs and ganglia and rendering them immobile (Juanes et al., 2006, and references therein). This means that the maximum water saturation that can be achieved during imbibition equals 1 - $S_{\text{gt}}$ (the trapped gas saturation). Here, we used measurements in sandpacks from Pentland et al. (2010) to determine $S_{\text{gt}}$ . In particular, we fitted Land (1968)’s model with the form $S_{\text{gt}}^{*}=\frac{S_{\text{gi}}^{*}}{1+CS_{\text{gi}}^{*}}$ , where $S_{\text{g}}^{*}=\frac{S_{\text{g}}}{1-S_{\text{wc}}}=1-S_{\text{w}}^{*}$ , $S_{\text{gi}}$ is the gas saturation at flow reversal, and $C$ is Land’s trapping coefficient with a value of 5.2 in our fit. Although Pentland et al. (2010) report that the best fit is achieved with the Aissaoui (1983) and Spiteri et al. (2008) models (cf. their Fig. 5), Land’s model was chosen here given that most relative permeability hysteresis models build on this one (see next paragraph).

Nonwetting phase trapping contributes to irreversibility of the relative permeability and capillary pressure curves (hysteresis). Here, we accounted for this mechanism in the gas relative permeability due to its importance in subsurface CO₂ migration (Juanes et al., 2006, and references therein). In particular, we used Land (1968)’s model to compute the bounding imbibition curve (see Fig. 3d), where $S_{\text{gt}}$ is obtained as described above, and Killough (1976)’s model to characterize the scanning curves. In Killough’s model, the scanning curves are reversible, such that the relative permeability at $S_{\text{g}}<S_{\text{gi}}$ no longer depends on the displacement type (Saló-Salgado et al., 2023, forthcoming).

3.2.2 Model II

This model had access to local, ex-situ measurements of single-phase petrophysical properties, i.e., porosity and intrinsic permeability (see sect. 2.1 and Tab. 1). Comparing with Tab. 3, it can be seen that our estimation for model I above was correct to the order of magnitude, but resulted in smaller values: porosity $\in[85,\,93]\%$ and permeability $\in[53,\,84]\%$ of the local measurements.

Capillary pressures and relative permeabilities were obtained using the same procedure described above for model I. The slight differences with respect to the curves shown in Fig. 3b,d,e come from the porosity and permeability values used in the Leverett scaling and to determine $S_{\text{wc}}$ , which were taken from Tab. 1 instead. The obtained curves for model II are provided in Fig 4.

3.2.3 Model III

This model was allowed access to all local, ex-situ measurements (see Tab. 1). Initial porosity and permeability remain unchanged with respect to model II. Capillary pressure curves were obtained by scaling the reference curve described in sect. 3.2.1 and shown in Fig. 3a using the measured entry pressure (sect. 2.1). The scaling followed the model $p_{\text{cs}}(S_{\text{w}})=p_{\text{cr}}(S_{\text{w}})\frac{p_{\text{e}}}{p_{\text{er}}}$ , where $p_{\text{e}}$ is the measured entry pressure for each sand, and $p_{\text{er}}$ is the reference curve entry pressure. The obtained curves are shown in Fig. 5a.

Relative permeabilities were computed following the same procedure described for model I above. In this case, however, each sand type was assigned the measured $S_{\text{wc}}$ and $S_{\text{gt}}$ values (see Tab. 1). This led to differences in both the drainage and imbibition curves, as shown in Fig. 5.

3.3 Model calibration

The mismatch between results obtained with each simulation model (I to III) and the validation experiment in tank 1 (A1, see sect. 2.4) is quantitatively assessed by comparing the following quantities (see Fig. 6):

1.

At $t=55$ min (end of injection in port $I_{1}$ ): Areas occupied by free-phase CO₂, and dyed water with dissolved CO₂ in the bottom F reservoir.
2.

At $t=154$ min (end of injection in port $I_{2}$ ): Areas occupied by free-phase CO₂, and dyed water with dissolved CO₂, in the middle and top F reservoirs.
3.

Time at which the first finger touches the tank bottom.
4.

Time at which the first finger (sinking from the top F reservoir) touches the middle C sand.

The experimental values for points 1-2 were obtained by computing areas from time-lapse images using a vector-graphics software. Careful visual inspection of color-enhanced images was used to distinguish between free-phase CO₂ (white) and dyed water with dissolved CO₂ (yellowish orange to red), and to identify the times for points 3-4 above. In the simulation models, the threshold gas saturation and CO₂ concentration in water used to compute areas were $S_{\text{g}}>10^{-3}$ and $C_{\text{CO}_{2}}>15\%(C_{\text{CO}_{2}}^{\text{max}})\approx 0.21$ [kg/m³], respectively. The $C$ value was chosen after a shape comparison of the region with dissolved CO₂. A smaller value of $C_{\text{CO}_{2}}>0.05$ [kg/m³] was selected to determine finger times for points 3 and 4 above. Fig. 6 shows an overview of the experimental values for points 2 and 3, while Fig. 11 in Sect. 4.2 shows the full comparison with the history-matched/calibrated simulation models.

The experiment was conducted first. Afterwards, the process consisted of running simulation models I to III, in parallel, starting with the petrophysical properties described in sect. 3.2. Given the number of uncertain variables (four petrophysical properties for each lithological unit, the diffusion coefficient and the injection rate) and the time required to complete a single simulation, a manual history matching method was employed. At the end of each run, quantities 1-4 above were compared and one or more properties were manually changed based on observed mismatch and domain knowledge. During the first few runs, only quantities 1 and 2 above were compared. After obtaining a satisfactory areal match, petrophysical properties were further adjusted to match quantities 3 and 4.

4 Results

In sect. 4.1, we present the results of the first simulation of experiment A1 with each model and property values detailed in sect. 3.2. Then, we detail the calibration of simulation models using experiment A1, and assess the value of local data to history-match CO₂ storage simulation models (sect. 4.2). Finally, we apply these matched models to experiment A2, analog for a longer injection in the same geology (sect. 4.3.1), and to experiment B1, analog for a larger-scale injection in a different geologic setting (sect. 4.3.2). We use simulations of experiments A2 and B1 to assess the predictive ability of simulation models in different conditions.

4.1 Initial model results

Fig. 7 shows the comparison between experiment A1 and the first run with each model, at times indicated in sect. 3.3. Numerous differences are evident between the experiment and models I and II, while model III is much closer to the experiment. In particular, models I and II overestimate the extent of CO₂-rich brine and underestimate the amount of gaseous CO₂ in all F reservoirs (refer to Fig. 6 for location). Model III approximates much better the areal extent of gaseous CO₂ in all regions, as well as the CO₂-rich brine in the middle and upper F reservoirs. Model II provides the closest finger migration times (points 3 and 4 in sect. 3.3), although this was not evaluated in the first run, as discussed below.

Petrophysical properties for models I and II were obtained from references in sect. 3.2, which also used silica sands with similar grain sizes. However, despite the relatively homogeneous nature of our quartz sands, model III performed significantly better. This result stems from natural sand variability and highlights the difficulty in establishing general, representative elementary volume-scale properties for porous media (see, for instance, Hommel et al., 2018; Schulz et al., 2019, for a discussion on intrinsic permeability). Additionally, results in Fig. 7 highlight the need for conducting sand/rock-specific measurements, even in the case of well-sorted, homogeneous sediments.

4.2 Manual history matching and value of local data

Fig. 8 shows convergence of areas occupied by free gas ( $A_{\text{g}}$ ) and water with dissolved CO₂ ( $A_{\text{d}}$ ), according to sect. 3.3. Each iteration corresponds to a successive model with updated parameters, and the different F sand regions evaluated in each panel (a) to (f) are provided in Fig. 6. With the exception of $A_{\text{d}}$ in the upper compartment, model III performs well since the beginning, and all areas were satisfactorily matched after four iterations. Conversely, model I and II were significantly off the experimental reference during the first few iterations. Model II, however, performed well after five iterations, while model I required seven iterations to give satisfactory areal estimates. The mean absolute error (MAE) over the six areal quantities presented in Fig. 8 is evaluated in Fig. 9, where it can be seen that, while all models are accurate towards the end (MAE $\in[5-10]$ cm²), that required a six-fold improvement in models I and II, but only two-fold in model III. As mentioned in sect. 3.3, $C_{\text{CO}_{2}}>15\%(C_{\text{CO}_{2}}^{\text{max}})\approx 0.21$ [kg/m³] was used as threshold to determine areas. While the absolute values and error would change with a different $C_{\text{CO}_{2}}$ threshold, we checked that the relative performance of our calibrated models does not with both $C_{\text{CO}_{2}}>0.01$ and 0.1 [kg/m³].

Convergence of quantities 3 and 4 in sect. 3.3, the times at which the first finger touches the rig bottom and the middle C sand, respectively, are provided in Fig. 10. These times were only evaluated after a satisfactory areal match for quantities in Fig. 8 was achieved. Therefore, areas no longer change much in the last few iterations in Fig. 8. In Fig. 10, it can be seen that model II and III, which incorporated local intrinsic permeability measurements, were significantly closer to our experimental reference than model I. Initially, however, we observed that sinking of gravity fingers in the experiment was faster than our model values by a factor of $\approx 2$ . A satisfactory match of all quantities evaluated was achieved after 11, 8, and 7 iterations for models I-III, respectively.

Overall, we find that model III, with access to local single-phase and multiphase flow properties, is closer to the experimental reference from the start. Model I started farthest, and required significantly more effort for calibration. After the calibration process, all models achieve a similar, accurate estimate of the evaluated quantities (Fig. 9).

The calibration shown in Fig. 8, 9, 10 employs $D=10^{-9}$ m²/s in all models. Injection rates ( $I_{\text{R}}$ ) started at 2.0 ml/min for all three models, and were 1.6 ml/min, 1.8 ml/min and 1.75 ml/min, respectively, at the end of the calibration. $I_{\text{R}}$ is slightly different because the goal was to obtain the best match with each model, considering $I_{\text{R}}$ to be an uncertain variable. In sect. 4.3 below, the same $I_{\text{R}}$ is used to make predictions with all three models.

Tab. 4 compares the starting and final (matched) key petrophysical variables for each model. The models were successfully calibrated by adjusting intrinsic permeability and the capillary pressure curves (same shape, but scaled to higher or lower $p_{\text{e}}$ ) only. It was found that CO₂ migration was most sensitive to the properties of the F sand, were most of the CO₂ migration occurs, as well as the ESF seal, which structurally traps the CO₂ plume. In our matched models, $p_{\text{e}}$ of ESF is about twice the measured value; this was required because the minimum saturation at which we can define $p_{\text{e}}$ and ensure convergence is $S_{\text{g}}\approx 10^{-4}$ . Reality, however, is closer to a jump in $p_{\text{c}}$ from 0 to $p_{\text{e}}$ at an infinitesimally small $S_{\text{g}}$ . Additionally, we found that our simulations performed better when using different values for the C and F sands in different model regions. In the case of the C sand, the explanation lies in the fault construction process, which may reduce porosity with respect to “natural” sedimentation of stratigraphic layers (Haugen et al., 2023, this issue). The increase in F sand permeability was required to match finger migration times, and is possibly compensating the absence of mechanical dispersion in the simulations. This is discussed in sect. 5. Our calibrated values are within the same order of magnitude of the ex-situ measurements (Tab. 4) and history-matched values for the porous medium in tank 2 (Landa-Marbán et al., 2023, forthcoming).

Table 4: Petrophysical properties for used quartz sands in experiment A1. Methodology for local measurements is provided by Haugen et al. (2023), while starting property modeling is described in sect. 3.2. For each sand, measured (first row), initial (superscript i) and final (superscript f) values for each of our models is shown. For sand C, the second permeability value refers to the fault, if different from the rest. For sand F, the second permeability value refers to the middle F layer, if different from the rest. For model III, where property values are different, m

{}_{\text{III},1}

refers to the calibration with

D=10^{-9}

m²/s and m

{}_{\text{III},3}

refers to

D=3\times 10^{-9}

m²/s.

Sand type / model	$\phi$ [-]	$k$ [D]	$p_{\text{e}}$ [mbar]	$S_{\text{wc}}$ [-]	$S_{\text{gt}}$ [-]
ESF	0.435	44	15	0.32	0.14
m ${}_{\text{I}}^{\text{i}}$	0.37	25	31.4	0.09	0.1468
m ${}_{\text{I}}^{\text{f}}$	0.37	6	31.4	0.09	0.1468
m ${}_{\text{II}}^{\text{i}}$	0.435	44	25.6	0.09	0.1468
m ${}_{\text{II}}^{\text{f}}$	0.435	44	25.6	0.09	0.1468
m ${}_{\text{III}}^{\text{i}}$	0.435	44	15	0.32	0.14
m ${}_{\text{III}}^{\text{f}}$	0.435	15	30	0.32	0.14
C	0.435	473	3	0.14	0.1
m ${}_{\text{I}}^{\text{i}}$	0.38	293	9.3	0.03	0.1565
m ${}_{\text{I}}^{\text{f}}$	0.38	293, 27	4.6	0.03	0.1565
m ${}_{\text{II}}^{\text{i}}$	0.435	473	7.8	0.03	0.1565
m ${}_{\text{II}}^{\text{f}}$	0.435	473, 158	2.6	0.03	0.1565
m ${}_{\text{III}}^{\text{i}}$	0.435	473	3	0.14	0.1
m ${}_{\text{III}}^{\text{f}}$	0.435	473, 118	4.5	0.14	0.1
E	0.45	2005	-	0.12	0.06
m ${}_{\text{I}}^{\text{i}}$	0.39	1528	4.1	0.01	0.16
m ${}_{\text{I}}^{\text{f}}$	0.39	1528	0.5	0.01	0.16
m ${}_{\text{II}}^{\text{i}}$	0.45	2005	3.86	0.01	0.16
m ${}_{\text{II}}^{\text{f}}$	0.45	3008	0.58	0.01	0.16
m ${}_{\text{III}}^{\text{i}}$	0.45	2005	0.33	0.12	0.06
m ${}_{\text{III},1}^{\text{f}}$	0.45	2406	0.33	0.12	0.06
m ${}_{\text{III},3}^{\text{f}}$	0.45	3208	0.33	0.12	0.06
F	0.44	4259	-	0.12	0.13
m ${}_{\text{I}}^{\text{i}}$	0.39	2277	3.3	0.01	0.16
m ${}_{\text{I}}^{\text{f}}$	0.39	6540, 2907	0	0.01	0.16
m ${}_{\text{II}}^{\text{i}}$	0.44	4259	2.62	0.01	0.16
m ${}_{\text{II}}^{\text{f}}$	0.44	6814, 4259	0	0.01	0.16
m ${}_{\text{III}}^{\text{i}}$	0.44	4259	0	0.12	0.13
m ${}_{\text{III},1}^{\text{f}}$	0.44	7240, 4685	0	0.12	0.13
m ${}_{\text{III},3}^{\text{f}}$	0.44	9796, 4259	0	0.12	0.13

Fig. 11 shows gas saturation ( $S_{\text{g}}$ ) and CO₂ concentration ( $C_{\text{CO}_{2}}$ ) maps at times at which quantities 1-4 described in sect. 3.3 are evaluated. Snapshots are provided for model III only, since all three calibrated models were qualitatively very similar. It can be seen that CO₂ migration is successfully approximated by our numerical model. In detail, however, some differences are apparent: Firstly, sinking of CO2-rich water from the bottom injector and horizontal migration along the bottom of the rig is faster in the model. This is due to the higher permeability that our numerical model requires in order to match the gravity fingering advance (cf. Tab. 4). Secondly, the experiment shows that denser, CO₂-rich water sinks with a rather compact front and closely spaced, wide fingers. Our model with constant $D=10^{-9}$ m²/s approximates all gravity-driven migration of the CO₂-rich water through thinner fingers, with the CO₂-saturated region receding with $S_{\text{g}}$ . To better represent fingering widths, we also matched model III with $D=3\times 10^{-9}$ m²/s, used in sect. 4.3.2.

4.3 Transferability: model predictions

A key question after history matching a flow simulation model is whether the physical description has actually been improved, or whether parameters have been modified to match a set of specific observations only. By applying the history-matched models to a different injection protocol (experiment A2 in tank 1; refer to Tab. 2), and subsequently to a different geometry (experiment B1 in tank 2), this can be assessed to some extent.

4.3.1 Analog for a longer CO₂ injection in the same geologic setting

This case illustrates performance of our history matched models in a much longer injection in the same geology (experiment A2). Before simulating this case, we observed that the trapped gas column against the fault in the experiment was different than what could be achieved with our previous $p_{\text{e}}$ for models I-III (Tab. 4). Because the capillary properties of the C sand in the fault were not directly involved in experiment A1, we increased $p_{\text{e}}$ in our calibrated models for that specific region ( $p_{\text{e}}$ = 5 mbar against the lower F sand, and 3.5 mbar against the middle F sand). All other parameters were taken from the values calibrated to match experiment A1.

Evaluation was performed at the end of injection, at $t=4$ h 48min, with a single run with models I-III. $I_{\text{R}}$ and $D$ were set to the same value in all three models: 1.7 ml/min and $10^{-9}$ m/s², respectively. The experimental result is shown in Fig. 12a, while the simulation with model III is depicted in Fig. 12b,c. We observe that the general distribution of CO₂ is close to the experimental truth. However, the experiment shows a compact sinking front of the CO₂-rich water without fingers; in our model, gravity fingering is apparent at this stage and fingers are close to the bottom of the rig. Additionally, CO₂-saturated brine touches the right boundary in the upper F reservoir, which does not occur in the experiment. This is due to capillary breach of the C sand above the middle F reservoir, as shown in Fig. 12b, and can be avoided by reducing the gas saturation value at which $p_{\text{e}}$ is defined, or by increasing $p_{\text{e}}$ .

The comparison of areal quantities is provided in Fig. 13, and shows good agreement between the experiment and simulation models. Model II (MAE = 16 cm²) and III (MAE = 14.54 cm²) perform similarly and slightly better than model I (MAE = 20.18 cm²), but there are no marked differences.

4.3.2 Analog for a larger-scale CO₂ injection in a different geologic setting

Finally, we compare the performance of our calibrated models against experiment B1, conducted in a larger-scale, more complex geology (Fig. 1e) (Flemisch et al., 2023, this issue). Similar to sect. 4.3.1, our goal is to assess the predictive capability of our calibrated models—without changing their properties. However, given that sand D controls migration in the lower fault (see Fig. 2e) and it was not present in our calibrated models, we allowed one change for models I and II, which did not have access to local $p_{\text{c}}$ measurements. This means that we ran an initial simulation of this experiment with model I and II, and then adjusted the $p_{\text{c}}$ curve of the D sand. The selected curve lies at $\approx\frac{1}{3}$ of the $p_{\text{c}}(S_{\text{w}})$ shown in Fig. 3 and Fig. 4, respectively.

Next, we evaluated the performance of models I-III by comparing them to the experimental truth after a single run. Evaluation is performed over the total duration of the experiment (120 h), which is simulated with the same $I_{\text{R}}$ (10 ml/min) and $D$ ( $10^{-9}$ m²/s) in all three models. Additionally, a run with $D=3\times 10^{-9}$ m²/s was completed with model III to better approximate finger widths, as noted in sect. 4.2.

Gas saturation and CO₂ concentration maps at the end of injection with model I are shown in Fig. 14a and Fig. 14b, respectively. The full visual comparison is provided in Fig. 15. We make the following observations:

•

At the end of injection ( $t=5$ h), all three models predict some migration of CO₂ into box B. Model II (Fig. 15c) and III (Fig. 15d) underestimate the amount of CO₂, while model I (Fig. 15b) overestimates the amount of CO₂ in the top C sand.
•

Also at the end of injection, all models predict faster sinking of the CO₂-charged water tongue arising from the lower injector. This is due to the higher F sand permeability required to match finger advance (see sect. 4.2), particularly in model III with $D=3\times 10^{-9}$ m²/s.
•

The speed at which CO2-rich fingers sink is slightly faster in our models, compared to the experiment. As expected, model III, with a higher diffusion coefficient, displays thicker fingers, with closer widths to the experiment. Similar to our previous observations, the numerical models cannot approximate the compact, CO₂-rich water front closely trailing the fingers.
•

Dissolution of CO₂ is underestimated by models I and II, while it is closer, but overestimated, by model III.

Consistent with our approach described in sect. 3.3, quantitative analysis is provided by means of areal quantities over time in Fig. 16. Experimental values were obtained via segmentation of timelapse images, and the data was reported on a 1 $\times$ 1 cm grid where 0 is pure water, 1 is water with dissolved CO₂, and 2 is gaseous CO₂. The segmentation procedure is explained in Nordbotten et al. (2023), this issue. We then obtained the areas of each phase within box A and B to generate Fig. 16 (refer to Fig. 14a for box location).

In box A, which contains the main F reservoir and ESF seal, we observe very accurate areas during injection. Afterwards, model III with $D=3\times 10^{-9}$ m²/s continues to follow the experiment closely, whereas the others overestimate gaseous CO₂. Note that the PVT properties of our fluids are the same in all models; differences arise due to (1) higher sand F $S_{\text{wc}}$ in model III, and higher sand F $k$ in model II and especially model III ( $D=3\times 10^{-9}$ m²/s), compared to model I, which allow greater convective mixing (Ennis-King and Paterson, 2005)(Tab. 4); and (2) lower $p_{\text{e}}$ and higher $k$ of sand ESF in model II (Tab. 4), which allows some CO₂ migration into the seal (Fig. 15). In box B (Fig. 16d-f), model I and model III with $D=10^{-9}$ m²/s are able to approximately track the experimental truth during injection. However, our models without dispersion cannot capture the areal increase of CO₂-rich water that occurs afterwards (cf. Fig. 15).

Further evaluation of simulation model performance, including comparison of mass quantities and error measures, is provided in Appendix A. From this analysis (sect. 4.3 and Appendix A), we find that all three calibrated models provide a reasonable prediction of the experiment. Qualitatively, they approximate well CO₂ migration, seal capacity, and onset of convective mixing. Quantitatively, our models are able to accurately estimate specific quantities during the injection phase, yet they accumulate higher errors later on (Fig. 16 and Appendix A). Overall, model III performs marginally better than model I and II (Appendix A), but, similar to sect. 4.3.1, the differences are not marked or consistent. The system is very sensitive to variations in petrophysical properties such as capillary pressure; therefore, in regions dominated by structures with uncertain properties (like box B in experiment B1), our estimates were less accurate.

In summary, calibrated models are transferable to a different operational setting or geologic structure, as long as sediments and trap systems remain similar (experiment A2 and box A in experiment B1). Where reservoir connectivity is provided by heterogeneous structures with uncertain properties, accurate deterministic estimates of CO₂ migration are unlikely; models calibrated elsewhere (experiment A1) did not perform well in our test (box B in experiment B1). Given unlimited computational time, the predictive ability of numerical models calibrated with published data appears similar to those having access to local measurements. The main value of local data lies in reducing the time required for history matching.

5 Discussion

Our experiments and numerical simulations of CO₂ injection and migration in unconsolidated sands at atmospheric $p,\,T$ conditions capture the CO₂-water system dynamics at short to intermediate timescales: buoyancy-driven flow and structural trapping (Bachu et al., 1994; Bryant et al., 2008; Hesse and Woods, 2010; Szulczewski et al., 2013), residual trapping (Juanes et al., 2006; Burnside and Naylor, 2014) and convective mixing and dissolution trapping (Weir et al., 1996; Ennis-King and Paterson, 2005; Riaz et al., 2006; Neufeld et al., 2010; Hidalgo et al., 2012; MacMinn and Juanes, 2013; Szulczewski et al., 2013). Due to the very large sand permeability ( $10^{2}-10^{4}$ D), the latter dominates CO₂ trapping.

The strong buoyancy and high permeability lead to persistent appearance and disappearance of fluid phases, as the gas migrates upward and dissolves in the water; this renders the nonlinear solver unstable, and, coupled with the nonlinearities introduced by relative permeability hystreresis, makes this problem difficult to solve numerically (e.g., Lie, 2019). As hinted in sect. 3.1, this was addressed by optimizing linear solver time, reducing the timestep length, increasing the number of timestep cuts and relaxing MRST’s maximum normalized residual where required.

In a 2D isotropic medium and assuming uniform flow, the hydrodynamic dispersion coefficient ( $\underline{\underline{D}}_{\text{h}}$ ) can be modeled as $\underline{\underline{D}}_{\text{h}}=\bigl{[}\begin{smallmatrix}\alpha_{\text{L}}\overline{u}&0\\ 0&\alpha_{\text{T}}\overline{u}\end{smallmatrix}\bigr{]}$ , where $\alpha_{\text{L}}$ and $\alpha_{\text{T}}$ are the longitudinal and transverse dispersivity, respectively, and $\overline{u}$ is the average Darcy velocity (Bear, 1972). Assuming dispersivities $\geq 10^{-2}-10^{-3}$ m (Garabedian et al., 1991; Gelhar et al., 1992; Schulze-Makuch, 2005) and $\overline{u}\approx 3\times 10^{-6}$ m/s (from our simulations), we get $\underline{\underline{D}}_{\text{h}}\in[3\times 10^{-9},\,3\times 10^{-8}]$ m²/s or larger; this means that $\underline{\underline{D}}_{\text{h}}\geq D$ for the timescales considered (Riaz et al., 2004; Rezk et al., 2022). Previous work suggested that hydrodynamic dispersion in homogeneous sediments can be accounted for by increasing $D$ (Riaz et al., 2004, 2006), as done here. However, our analysis shows that the spreading of CO₂-rich water during convective mixing can be loosely, but not accurately, represented by molecular diffusion. Given (1) the dominance of convective mixing on solubility trapping (Ennis-King and Paterson, 2005; Neufeld et al., 2010; MacMinn and Juanes, 2013); (2) heterogeneity of many natural reservoirs, which increases the importance of dispersion (Riaz et al., 2006; Bear, 2018); and (3) the acceleration of CO₂ dissolution due to dispersion, as observed here and by others (e.g., Hidalgo and Carrera, 2009), it is important to quantify the balance between diffusion and dispersion to estimate CO₂ trapping.

The models developed in this paper retained some error at the end of the calibration phase, which is a known problem of manual history matching (Oliver and Chen, 2011). This approach was chosen given the low geometric uncertainty; the availability of both local petrophysical measurements and the experimental truth; and the long computational times and number of uncertain variables. Consistent with previous results (e.g., Fisher and Jolley, 2007), results show that model II and III, which had access to local petrophysical measurements, achieved faster match to the experimental truth than model I (sect. 4.2). However, all models were similarly accurate when making predictions (sect. 4.3). Subsurface formations are significantly more heterogeneous than the silica sands employed here; together with time constraints, this may explain why existing literature emphasizes the importance of accurate geologic representation and reservoir model parametrization, including local data, to achieve successful history matching, and, especially, forecasting (e.g., Gosselin et al., 2003; Fisher and Jolley, 2007; Myers et al., 2007; Kam et al., 2015; Avansi et al., 2016).

Because of the availability of a ground truth to assess model predictions, we did not quantify uncertainty in the history matched models here. In general, however, this is necessary to manage reservoir operations (e.g., Aanonsen et al., 2009; Oliver and Chen, 2011; Jagalur-Mohan et al., 2018; Jin et al., 2019; Liu and Durlofsky, 2020; Santoso et al., 2021; Landa-Marbán et al., 2023, forthcoming, and references therein). It is also important to note that history-matched models may have grid-size dependencies (see Appendix B), which may require that the grid used to make predictions, if different or encompassing additional regions, maintain a similar resolution. Finally, multiphase flow in poorly-lithified sediments is non-unique (Haugen et al., 2023, this issue), which also contributes to uncertainty. Due to the sensitivity of the CO₂-brine system to relatively small reservoir parameter changes, as observed in sect. 4.3.2, it seems prudent to adopt a probabilistic perspective when estimating subsurface CO₂ migration (see Flemisch et al., 2023, this issue).

6 Conclusions

We performed experiments (sect. 2) and numerical simulations (sect. 3) of CO₂ migration in poorly-lithified, siliciclastic sediments at the meter scale. Three simulation models, with access to different levels of local data, were manually history-matched to the experiments (sect. 4.1, 4.2), and then used to make predictions (sect. 4.3). The main findings are:

1.

The time required to history match model III (access to both single-phase and multiphase measurements) is lower than model II (access to local single-phase measurements), which is lower than model I (no access to local petrophysical measurements).
2.

All simulation models achieve a satisfactory qualitative match throughout the experiments. Quantitatively, overall prediction accuracy of models I-III is similar: our models are close to the experimental truth during CO₂ injection, and accumulate larger errors afterwards, especially in regions where heterogeneous structures control CO₂ migration.
3.

The addition of a constant molecular diffusion coefficient allows matching convective finger widths to experimental observations. However, simulations without dispersion cannot approximate the compact, CO₂-rich sinking front closely trailing convective fingers in our experiments.

Simulation models were not always accurate. Given the degree of control in our study, it seems prudent to quantify uncertainty when assessing subsurface CO₂ migration in the field using numerical models. Obtained results suggest that confidence can be increased by obtaining local data (time is finite), quantifying petrophysical parameter uncertainty, testing sensitivity to petrophysical parameters in different model regions, including post-injection data when history matching, and incorporating multiple scenarios of CO₂ migration, particularly where heterogeneous structures are at play.

\bmhead

Acknowledgments LS gratefully acknowledges laboratory support provided by UiB Engineer Emil Bang Larsen and image processing support provided by UiB PhD student Benyamine Benali. LS and RJ are grateful to Olav Møyner and the MRST team at SINTEF for their guidance to implement new functionality and continuous support with the MATLAB Reservoir Simulation Toolbox. A special thanks goes to Robert Gawthorpe, Atle Rotevatn and Casey Nixon for helpful comments on the geology of North Sea reservoirs, which were key to build the geometry in tank 2. The authors also acknowledge support by the following organizations (refer to ‘Funding’ below for details): ExxonMobil, “la Caixa” Foundation, Research Council of Norway (RCN), Akademia.

Statements and Declarations

•

Funding This work was supported by ExxonMobil through the project “Modeling and Mitigation of Induced Seismicity and Fault Leakage during CO₂ storage”. LS gratefully acknowledges the support of a fellowship from “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/EU21/11890139. MH is funded by the Research Council of Norway (RCN) project no. 280341. KE, MH, JMN and MF are partly funded by the Centre for Sustainable Subsurface Resources, RCN project no. 331841. The work of JMN and MF is funded in part through the Akademia project “FracFlow”.
•

Competing interests The authors declare that they have no conflict of interest.
•

Availability of data and materials The dataset for experiment B1 can be obtained from https://github.com/fluidflower. The remaining experimental and simulation data are available from the corresponding authors on reasonable request.
•

Authors’ contributions JMN, MF and RJ designed the study and acquired the funding. JMN, MF and KE conceptualized, designed and built the FluidFlower rigs. MH, KE and MF conducted the ex-situ sand property measurements. LS designed the experiments in tank 1. JMN, MF, MH and KE designed the experiments in tank 2. MH, KE, LS and MF conducted the experiments. LS developed the simulation models and conducted the simulations. LS, JMN and RJ performed the simulation analysis. LS wrote the paper, with inputs from all authors.

Appendix A Additional analysis of simulation model performance against experiment B1

First, we provide the total mass of CO₂ in the computational domain in Fig. 17, and the mass in boxes A and B in Fig. 18.

Next, in Tab. 5, the following measures are compared with quantities estimated from the experiment via segmentation of timelapse images (Nordbotten et al., 2023). These measures correspond to the sparse data requested to participants of the FluidFlower international benchmark study (Flemisch et al., 2023):

1.

time of maximum mobile free phase in box A
2.

mass of mobile CO ${}_{2(\text{g})}$ , immobile CO ${}_{2(\text{g})}$ , dissolved CO₂, and CO₂ in the seal (in any phase), in box A, 72 h after injection start (2a-d)
3.

the same quantities as 2. for box B (3a-d)
4.

time at which $M$ (defined below) exceeds 110% of the width of box C
5.

total mass of CO₂ in the ESF seal, in box A, at $t=120$ h

Convective mixing in box C (see Fig. 1e) is reported as the integral of the magnitude of the gradient in relative concentration of dissolved CO₂ (Flemisch et al., 2023):

M(t)=\int_{\text{C}}\bigg{\lvert}\nabla\bigg{(}\frac{\chi_{\text{CO}_{2}}^{\text{w}}}{\chi_{\text{CO}_{2}}^{\text{w,max}}}\bigg{)}\bigg{\rvert}d\mathbf{x}

(1)

where $\chi_{\text{CO}_{2}}^{\text{w}}$ is the mass fraction of CO₂ in water, and the dissolution limit is $\chi_{\text{CO}_{2}}^{\text{w,max}}$ . Note that quantity 4, based on $M$ , cannot be provided with full accuracy based on experimental data, so an uncertain lower and upper bound is provided instead. Therefore, we assigned an error ( $\varepsilon$ ) = 0 in Tab. 5.

Relative error is evaluated with respect to the experimental mean ( $\overline{E}$ ) as $\varepsilon_{i}(\%)=100\times\frac{\lvert\overline{E_{i}}-\text{m}_{J,i}\rvert}{\overline{E_{i}}}$ , where $i$ is a given measure and $J$ refers to any of the models I-III. In Tab. 5, it can be seen that all models accumulate some error in most of the quantities reported. The maximum errors are $\approx 140\%$ for models I-II and $<100\%$ for model III. Model I performs better in the uncertain region (box B; see sect. 4.3.2 as well), while models II and III are more accurate in box A, the region where the calibration performed with experiment A1 is more meaningful. Overall, model III does marginally better.

Table 5: Sparse data comparison between experiment B1 in tank 2 and simulation results with models I-III. Experimental mean and standard deviation were obtained from six experimental runs with identical protocol, while the results for models I-III are for a single run with each matched model. For m

{}_{\text{III}}

, two cases are shown:

D=10^{-9}

m²/s (m

{}_{\text{III},1}

) and

D=3\times 10^{-9}

m²/s (m

{}_{\text{III},3}

). Experimental quantity 4 is reported using a lower and upper bound due to high uncertainty, so errors are not computed. See main text for measure description.

Measure	$\overline{E}$	$\sigma(E)$	m ${}_{\text{I}}$	$\varepsilon_{\text{I}}$ [%]	m ${}_{\text{II}}$	$\varepsilon_{\text{II}}$ [%]	m ${}_{\text{III},1}$	$\varepsilon_{\text{III},1}$ [%]	m ${}_{\text{III},3}$	$\varepsilon_{\text{III},3}$ [%]
1 [s]	14880	720	18000	21	17340	16.5	17880	20.2	17160	15.3
2a [g]	0.36	0.13	0.86	140.23	0.56	54.7	0.5	39.5	0.004	98.9
2b [g]	0	0	0	0	0	0	0	0	0	0
2c [g]	3.5	0.08	2.12	39.5	3.12	10.7	2.9	16.9	3.61	3.0
2d [g]	-	-	0.43	n/a	0.96	n/a	0.73	n/a	0.69	n/a
3a [g]	0	0	0	0	0	0	0	0	0	0
3b [g]	0	0	0	0	0	0	0	0	0	0
3c [g]	0.55	0.32	0.38	31.7	0.072	86.9	0.17	68.5	0.09	83.9
3d [g]	n/a	n/a	0.0015	n/a	0.0053	n/a	0.0021	n/a	0.0039	n/a
4 [s]	[12180, 17990]	[438, 2261]	15000	0	15600	0	19200	0	15600	0
5 [g]	0.38	0.047	0.52	36.2	0.92	142.2	0.72	89.1	0.62	62.7
$\overline{\varepsilon}$ [%]	n/a	n/a	n/a	29.8	n/a	34.6	n/a	26.0	n/a	29.3

Appendix B Comparison of simulation results with multiple grid resolutions

This section provides two comparisons of concentration maps obtained with model III after the calibration presented in sect. 4.2:

1.

For experiment A1, we compare two grid sizes: $h=4$ mm, as shown in the paper, and a coarser grid with $h=8$ mm (Fig. 19).
2.

For experiment B1, we compare three grid sizes: $h=5$ mm, as shown in the paper, and two coarser grids with $h=10$ mm and $h=20$ mm, respectively (Fig. 20).

It can be seen that, for the calibrated parameter set (Tab. 4), the coarser models maintain a general agreement with the finer ones (and the experimental solution). However, some differences are clear even in this qualitative comparison, including (1) smaller extent of the CO₂ plume, (2) lower dissolution, (3) lower number of fingers and finger widths, and (4) different CO₂-rich finger sinking speed. Therefore, the calibration process is somewhat cell-size dependent, which has implications for applying history matched models from e.g., pilot tests to field-scale CO₂ storage projects.

References

\bibcommenthead
Aanonsen et al. (2009) Aanonsen, S.I., G. Nœvdal, D.S. Oliver, A.C. Reynolds, and B. Vallès. 2009. The ensemble Kalman filter in reservoir engineering—a review. SPE Journal 14(03): 393–412. 10.2118/117274-PA .
Aissaoui (1983) Aissaoui, A. 1983. Etude théorique et expérimentale de l’hystérésis des pressions capillaries et des perméabilitiés relatives en vue du stockage souterrain de gaz. Ph. D. thesis, Ecole des Mines de Paris, Paris.
Alghannam and Juanes (2020) Alghannam, M. and R. Juanes. 2020. Understanding rate effects in injection-induced earthquakes. Nature Communications 11(1): 3053. 10.1038/s41467-020-16860-y .
Askar et al. (2021) Askar, A.H., T.H. Illangasekare, and A.M.C. Ilie. 2021. Monitoring brine leakage from deep geologic formations storing carbon dioxide: Design framework validation using intermediate-scale experiment. Water Resources Research 57(12): e2021WR031005. https://doi.org/10.1029/2021WR031005 .
Avansi et al. (2016) Avansi, G.D., C. Maschio, and D.J. Schiozer. 2016. Simultaneous history-matching approach by use of reservoir-characterization and reservoir-simulation studies. SPE Reservoir Evaluation & Engineering 19(04): 694–712. 10.2118/179740-PA .
Bachu et al. (1994) Bachu, S., W. Gunter, and E. Perkins. 1994. Aquifer disposal of CO₂: Hydrodynamic and mineral trapping. Energy Conversion and Management 35(4): 269–279. 10.1016/0196-8904(94)90060-4 .
Bear (1972) Bear, J. 1972. Dynamics of Fluids in Porous Media. Dover. 780 pp.
Bear (2018) Bear, J. 2018. Modeling Phenomena of Flow and Transport in Porous Media, Volume 31. Springer. 761 pp.
Beard and Weyl (1973) Beard, D. and P. Weyl. 1973. Influence of texture on porosity and permeability of unconsolidated sand. AAPG Bulletin 57(2): 349–369. 10.1306/819A4272-16C5-11D7-8645000102C1865D .
Berge et al. (2019) Berge, R.L., Ø.S. Klemetsdal, and K.A. Lie. 2019. Unstructured Voronoi grids conforming to lower dimensional objects. Computational Geosciences 23: 169–188. 10.1007/s10596-018-9790-0 .
Berge et al. (2021) Berge, R.L., Ø.S. Klemetsdal, B. Skaflestad, and K.A. Lie. 2021. The UPR module. https://github.com/rbe051/UPR.
Brock and Orr (1991) Brock, D. and F. Orr. 1991. Flow visualization of viscous fingering in heterogeneous porous media. SPE Annual Technical Conference and Exhibition (SPE-22614-MS). 10.2118/22614-MS .
Brooks and Corey (1964) Brooks, R.H. and A.T. Corey 1964. Hydraulic properties of porous media. Technical report, Colorado State University, Fort Collins, CO.
Bryant et al. (2008) Bryant, S.L., S. Lakshminarasimhan, and G.A. Pope. 2008. Buoyancy-dominated multiphase flow and its effect on geological sequestration of CO₂. SPE Journal 13(04): 447–454. 10.2118/99938-PA .
Burnside and Naylor (2014) Burnside, N. and M. Naylor. 2014. Review and implications of relative permeability of CO₂/brine systems and residual trapping of CO₂. International Journal of Greenhouse Gas Control 23: 1–11. 10.1016/j.ijggc.2014.01.013 .
Caine et al. (1996) Caine, J.S., J.P. Evans, and C.B. Forster. 1996. Fault zone architecture and permeability structure. Geology 24(11): 1025–1028. 10.1130/0091-7613(1996)024¡1025:FZAAPS¿2.3.CO;2 .
Cappa and Rutqvist (2011) Cappa, F. and J. Rutqvist. 2011. Impact of CO₂ geological sequestration on the nucleation of earthquakes. Geophysical Research Letters 38(17): L17313. 10.1029/2011GL048487 .
Celia (2021) Celia, M.A. 2021. The role of subsurface engineering in the net-zero energy transition. In IOP Conference Series: Earth and Environmental Science, Volume 861, pp. 072017. IOP Publishing.
Childs et al. (2009) Childs, C., T. Manzocchi, J.J. Walsh, C.G. Bonson, A. Nicol, and M.P. Schöpfer. 2009. A geometric model of fault zone and fault rock thickness variations. Journal of Structural Geology 31(2): 117–127. 10.1016/j.jsg.2008.08.009 .
Cinar et al. (2006) Cinar, Y., K. Jessen, R. Berenblyum, R. Juanes, and F.M. Orr. 2006. An experimental and numerical investigation of crossflow effects in two-phase displacements. SPE Journal 11(02): 216–226. 10.2118/90568-PA .
Corey (1954) Corey, A.T. 1954. The interrelation between gas and oil relative permeabilities. Producers Monthly 19: 38–41 .
Demidov and Rossi (2018) Demidov, D. and R. Rossi. 2018. Subdomain deflation combined with local AMG: A case study using AMGCL library. arXiv preprint. 10.48550/arXiv.1710.03940 .
DiCarlo et al. (2000) DiCarlo, D.A., A. Sahni, and M.J. Blunt. 2000. Three-phase relative permeability of water-wet, oil-wet, and mixed-wet sandpacks. SPE Journal 5(01): 82–91. 10.2118/60767-PA .
Duan and Sun (2003) Duan, Z. and R. Sun. 2003. An improved model calculating CO₂ solubility in pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar. Chemical Geology 193(3-4): 257–271. 10.1016/S0009-2541(02)00263-2 .
Eikehaug et al. (2023) Eikehaug, K., M. Haugen, O. Folkvord, B. Benali, E.L. Larsen, A. Tinkova, A. Rotevatn, J.M. Nordbotten, and M. Fernø. 2023. Engineering meter-scale porous media flow experiments for quantitative studies of geological carbon sequestration. Transport In Porous Media (this issue) .
Ellsworth (2013) Ellsworth, W.L. 2013. Injection-induced earthquakes. Science 341(6142): 1225942. 10.1126/science.1225942 .
Ennis-King and Paterson (2005) Ennis-King, J. and L. Paterson. 2005. Role of convective mixing in the long-term storage of carbon dioxide in deep saline formations. SPE Journal 10(03): 349–356. 10.2118/84344-PA .
European Academies Science Advisory Council (EASAC) (2018) European Academies Science Advisory Council (EASAC). Policy Report 35, 2018. Negative emission technologies: What role in meeting Paris Agreement targets? https://easac.eu/publications/details/easac_net/. Accessed on 01.20.2023.
Fernández-García et al. (2004) Fernández-García, D., T.H. Illangasekare, and H. Rajaram. 2004. Conservative and sorptive forced-gradient and uniform flow tracer tests in a three-dimensional laboratory test aquifer. Water Resources Research 40(10): W10103. 10.1029/2004WR003112 .
Fernø et al. (2023) Fernø, M., M. Haugen, K. Eikehaug, O. Folkvord, B. Benali, J.W. Both, E. Storvik, C.W. Nixon, R.L. Gawthorpe, and J.M. Nordbotten. 2023. Room-scale CO₂ injections in a physical reservoir model with faults. Transport In Porous Media (this issue) .
Fisher and Jolley (2007) Fisher, Q. and S. Jolley. 2007. Treatment of faults in production simulation models. Geological Society, London, Special Publications 292: 219–233. 10.1144/SP292.13 .
Flemisch et al. (2023) Flemisch, B., J.M. Nordbotten, M. Fernø, R. Juanes, H. Class, M. Delshad, F. Doster, J. Ennis-King, J. Franc, S. Geiger, D. Gläser, C. Green, J. Gunning, H. Hajibeygi, S.J. Jackson, M. Jammoul, S. Karra, J. Li, S.K. Matthäi, T. Miller, Q. Shao, C. Spurin, P. Stauffer, H. Tchelepi, X. Tian, H. Viswanathan, D. Voskov, Y. Wang, M. Wapperom, M.F. Wheeler, A. Wilkins, A.A. Youssef, and Z. Zhang. 2023. The FluidFlower international benchmark study: Process, modeling results, and comparison to experimental data. Transport In Porous Media (this issue) .
Garabedian et al. (1991) Garabedian, S.P., D.R. LeBlanc, L.W. Gelhar, and M.A. Celia. 1991. Large-scale natural gradient tracer test in sand and gravel, Cape Cod, Massachusetts: 2. Analysis of spatial moments for a nonreactive tracer. Water Resources Research 27(5): 911–924. 10.1029/91WR00242 .
Gaucher and Lindley (1960) Gaucher, D. and D. Lindley. 1960. Waterflood performance in a stratified, five-spot reservoir-a scaled-model study. Transactions of the AIME 219(01): 208–215. 10.2118/1311-G .
Gelhar et al. (1992) Gelhar, L.W., C. Welty, and K.R. Rehfeldt. 1992. A critical review of data on field-scale dispersion in aquifers. Water Resources Research 28(7): 1955–1974. 10.1029/92WR00607 .
Gosselin et al. (2003) Gosselin, O., S. Aanonsen, I. Aavatsmark, A. Cominelli, R. Gonard, M. Kolasinski, F. Ferdinandi, L. Kovacic, and K. Neylon. 2003. History matching using time-lapse seismic (HUTS). SPE Annual Technical Conference and Exhibition (SPE 84464). 10.2118/84464-MS .
Hager et al. (2021) Hager, B.H., J. Dieterich, C. Frohlich, R. Juanes, S. Mantica, J.H. Shaw, F. Bottazzi, F. Caresani, D. Castineira, A. Cominelli, M. Meda, L. Osculati, S. Petroselli, and A. Plesch. 2021. A process-based approach to understanding and managing triggered seismicity. Nature 595(7869): 684–689. 10.1038/s41586-021-03668-z .
Hassanzadeh et al. (2008) Hassanzadeh, H., M. Pooladi-Darvish, A.M. Elsharkawy, D.W. Keith, and Y. Leonenko. 2008. Predicting PVT data for CO₂–brine mixtures for black-oil simulation of CO₂ geological storage. International Journal of Greenhouse Gas Control 2(1): 65–77. 10.1016/S1750-5836(07)00010-2 .
Haugen et al. (2023) Haugen, M., L. Saló-Salgado, K. Eikehaug, B. Benali, J.W. Both, E. Storvik, O. Folkvord, R. Juanes, J.M. Nordbotten, and M. Fernø. 2023. Physical variability in meter-scale laboratory CO₂ injections in faulted geometries. Transport In Porous Media (this issue) .
Heinemann et al. (1991) Heinemann, Z., C. Brand, M. Munka, and Y. Chen. 1991. Modeling reservoir geometry with irregular grids. SPE Reservoir Engineering 6(2): 225–232. 10.2118/18412-PA .
Hesse and Woods (2010) Hesse, M. and A. Woods. 2010. Buoyant dispersal of CO₂ during geological storage. Geophysical Research Letters 37(1): L01403. 10.1029/2009GL041128 .
Hidalgo and Carrera (2009) Hidalgo, J.J. and J. Carrera. 2009. Effect of dispersion on the onset of convection during CO₂ sequestration. Journal of Fluid Mechanics 640: 441–452. 10.1017/S0022112009991480 .
Hidalgo et al. (2012) Hidalgo, J.J., J. Fe, L. Cueto-Felgueroso, and R. Juanes. 2012. Scaling of convective mixing in porous media. Physical Review Letters 109(26): 264503. 10.1103/PhysRevLett.109.264503 .
Hommel et al. (2018) Hommel, J., E. Coltman, and H. Class. 2018. Porosity–permeability relations for evolving pore space: a review with a focus on (bio-) geochemically altered porous media. Transport in Porous Media 124(2): 589–629. 10.1007/s11242-018-1086-2 .
Ingram and Urai (1999) Ingram, G.M. and J.L. Urai. 1999. Top-seal leakage through faults and fractures: the role of mudrock properties. Geological Society, London, Special Publications 158: 125–135. 10.1144/GSL.SP.1999.158.01.10 .
Intergovernmental Panel on Climate Change (IPCC) (2022) Intergovernmental Panel on Climate Change (IPCC). 2022. Climate change 2022: Mitigation of climate change, In Contribution of Working Group III to the Sixth Assessment Report of the IPCC, eds. Shukla, P.R., J. Skea, R. Slade, A. Al Khourdajie, R. van Diemen, D. McCollum, M. Pathak, S. Some, P. Vyas, R. Fradera, M. Belkacemi, A. Hasija, G. Lisboa, S. Luz, and J. Malley. Cambridge University Press, Cambridge, UK and New York, NY, USA. 10.1017/9781009157926.
Jagalur-Mohan et al. (2018) Jagalur-Mohan, J., B. Jha, Z. Wang, R. Juanes, and Y. Marzouk. 2018. Inferring fault frictional and reservoir hydraulic properties from injection-induced seismicity. Geophysical Research Letters 45(3): 1313–1320. 10.1002/2017GL075925 .
Jin et al. (2019) Jin, L., H. Lu, and G. Wen. 2019. Fast uncertainty quantification of reservoir simulation with variational U-Net. arXiv:1907.00718. 10.48550/arXiv.1907.00718 .
Juanes et al. (2012) Juanes, R., B.H. Hager, and H.J. Herzog. 2012. No geologic evidence that seismicity causes fault leakage that would render large-scale carbon capture and storage unsuccessful. Proceedings of the National Academy of Sciences 109(52): E3623–E3623. 10.1073/pnas.1215026109 .
Juanes et al. (2006) Juanes, R., E. Spiteri, F. Orr Jr, and M. Blunt. 2006. Impact of relative permeability hysteresis on geological CO₂ storage. Water Resources Research 42(12): W12418. 10.1029/2005WR004806 .
Jung et al. (2014) Jung, N.H., W.S. Han, Z. Watson, J.P. Graham, and K.Y. Kim. 2014. Fault-controlled CO₂ leakage from natural reservoirs in the Colorado Plateau, East-Central Utah. Earth and Planetary Science Letters 403: 358–367. 10.1016/j.epsl.2014.07.012 .
Kam et al. (2015) Kam, P., M. Nadeem, A. Novlesky, A. Kumar, and E.N. Omatsone. 2015. Reservoir characterization and history matching of the Horn River shale: an integrated geoscience and reservoir-simulation approach. Journal of Canadian Petroleum Technology 54(06): 475–488. 10.2118/171611-PA .
Killough (1976) Killough, J.E. 1976. Reservoir simulation with history-dependent saturation functions. SPE Journal 16(01): 37–48. 10.2118/5106-PA .
Krogstad et al. (2015) Krogstad, S., K. Lie, O. Møyner, H.M. Nilsen, X. Raynaud, and B. Skaflestad. 2015. MRST-AD – an open-source framework for rapid prototyping and evaluation of reservoir simulation problems. SPE Reservoir Simulation Conference: 1–20. SPE-173317-MS .
Land (1968) Land, C.S. 1968. Calculation of imbibition relative permeability for two-and three-phase flow from rock properties. Society of Petroleum Engineers Journal 8(02): 149–156. 10.2118/1942-PA .
Landa-Marbán et al. (2023) Landa-Marbán, D., T.H. Sandve, J.W. Both, J.M. Nordbotten, and S.E. Gasda. 2023. Performance of an open-source image-based history matching framework for CO₂ storage. (forhtcoming) .
Lenhard et al. (1995) Lenhard, R., M. Oostrom, C. Simmons, and M. White. 1995. Investigation of density-dependent gas advection of trichloroethylene: Experiment and a model validation exercise. Journal of Contaminant Hydrology 19(1): 47–67. 10.1016/0169-7722(94)00055-M .
Leverett (1941) Leverett, M.C. 1941. Capillary behavior in porous solids. Transactions of the AIME 142(01): 152–169. 10.2118/941152-G .
Lie (2019) Lie, K.A. 2019. An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User Guide for the MATLAB Reservoir Simulation Toolbox (MRST). Cambridge, UK: Cambridge University Press.
Lie and Møyner (2021) Lie, K.A. and O. Møyner eds. 2021. Advanced Modeling with the MATLAB Reservoir Simulation Toolbox. Cambridge, UK: Cambridge University Press.
Liu and Durlofsky (2020) Liu, Y. and L.J. Durlofsky. 2020. Multilevel strategies and geological parameterizations for history matching complex reservoir models. SPE Journal 25(01): 081–104. 10.2118/193895-PA .
MacMinn and Juanes (2013) MacMinn, C.W. and R. Juanes. 2013. Buoyant currents arrested by convective dissolution. Geophysical Research Letters 40(10): 2017–2022. 10.1002/grl.50473 .
Marcucci et al. (2017) Marcucci, A., S. Kypreos, and E. Panos. 2017. The road to achieving the long-term Paris targets: energy transition and the role of direct air capture. Climatic Change 144(2): 181–193. 10.1007/s10584-017-2051-8 .
Myers et al. (2007) Myers, R.D., A. Allgood, A. Hjellbakk, P. Vrolijk, and N. Briedis. 2007. Testing fault transmissibility predictions in a structurally dominated reservoir: Ringhorne Field, Norway. Geological Society, London, Special Publications 292(1): 271–294. 10.1144/SP292.16 .
Neufeld et al. (2010) Neufeld, J.A., M.A. Hesse, A. Riaz, M.A. Hallworth, H.A. Tchelepi, and H.E. Huppert. 2010. Convective dissolution of carbon dioxide in saline aquifers. Geophysical Research Letters 37(22): L22404. 10.1029/2010GL044728 .
Nordbotten et al. (2023) Nordbotten, J.M., B. Benali, J.W. Both, B. Brattekås, E. Storvik, and M. Fernø. 2023. DarSIA: An open-source python toolbox for two-scale image processing for porous media. Transport In Porous Media (this issue) .
Nordbotten and Celia (2012) Nordbotten, J.M. and M.A. Celia. 2012. Geological Storage of CO₂: Modeling Approaches for Large-scale Simulation. Hoboken, New Jersey: Wiley. 242 pp.
Nordbotten et al. (2012) Nordbotten, J.M., B. Flemisch, S.E. Gasda, H. Nilsen, Y. Fan, G.E. Pickup, B. Wiese, M.A. Celia, H.K. Dahle, G.T. Eigestad, and K. Pruess. 2012. Uncertainties in practical simulation of CO₂ storage. International Journal of Greenhouse Gas Control 9: 234–242. 10.1016/j.ijggc.2012.03.007 .
Oliver and Chen (2011) Oliver, D.S. and Y. Chen. 2011. Recent progress on reservoir history matching: a review. Computational Geosciences 15(1): 185–221. 10.1007/s10596-010-9194-2 .
Pentland et al. (2010) Pentland, C.H., E. Itsekiri, S.K. Al Mansoori, S. Iglauer, B. Bijeljic, and M.J. Blunt. 2010. Measurement of nonwetting-phase trapping in sandpacks. SPE Journal 15(02): 274–281. 10.2118/115697-PA .
Plug and Bruining (2007) Plug, W.J. and J. Bruining. 2007. Capillary pressure for the sand–CO₂–water system under various pressure conditions. application to CO₂ sequestration. Advances in Water Resources 30(11): 2339–2353. 10.1016/j.advwatres.2007.05.010 .
Rezk et al. (2022) Rezk, M.G., J. Foroozesh, A. Abdulrahman, and J. Gholinezhad. 2022. CO₂ diffusion and dispersion in porous media: Review of advances in experimental measurements and mathematical models. Energy & Fuels 36(1): 133–155. 10.1021/acs.energyfuels.1c03552 .
Riaz et al. (2006) Riaz, A., M. Hesse, H. Tchelepi, and F. Orr. 2006. Onset of convection in a gravitationally unstable diffusive boundary layer in porous media. Journal of Fluid Mechanics 548: 87–111. 10.1017/S0022112005007494 .
Riaz et al. (2004) Riaz, A., C. Pankiewitz, and E. Meiburg. 2004. Linear stability of radial displacements in porous media: influence of velocity-induced dispersion and concentration-dependent diffusion. Physics of Fluids 16(10): 3592–3598. 10.1063/1.1775431 .
Saadatpoor et al. (2010) Saadatpoor, E., S.L. Bryant, and K. Sepehrnoori. 2010. New trapping mechanism in carbon sequestration. Transport in Porous Media 82(1): 3–17. 10.1007/s11242-009-9446-6 .
Saló-Salgado et al. (2023) Saló-Salgado, L., J.S. Davis, and R. Juanes. 2023. Fault permeability from stochastic modeling of clay smears. Geology 51(1): 91–95. 10.1130/G50739.1 .
Saló-Salgado et al. (2023) Saló-Salgado, L., O. Møyner, K.A. Lie, and R. Juanes. 2023. Three-dimensional simulation of geologic carbon sequestration with miscibility, relative permeability hysteresis and diffusion, In Advanced Modelling with the MATLAB Reservoir Simulation Toolbox II, eds. Lie, K.A. and O. Møyner. Cambridge University Press (forthcoming).
Santoso et al. (2021) Santoso, R.K., X. He, M. Alsinan, H. Kwak, and H. Hoteit. 2021. Bayesian long-short term memory for history matching in reservoir simulations. SPE Reservoir Simulation Conference (SPE-203976-MS). 10.2118/203976-MS .
Schlumberger (2014) Schlumberger 2014. ECLIPSE Technical description (2014.1 ed.). Schlumberger Ltd.
Schulz et al. (2019) Schulz, R., N. Ray, S. Zech, A. Rupp, and P. Knabner. 2019. Beyond Kozeny–Carman: predicting the permeability in porous media. Transport in Porous Media 130(2): 487–512. 10.1007/s11242-019-01321-y .
Schulze-Makuch (2005) Schulze-Makuch, D. 2005. Longitudinal dispersivity data and implications for scaling behavior. Groundwater 43(3): 443–456. 10.1111/j.1745-6584.2005.0051.x .
Silliman and Simpson (1987) Silliman, S. and E. Simpson. 1987. Laboratory evidence of the scale effect in dispersion of solutes in porous media. Water Resources Research 23(8): 1667–1673. 10.1029/WR023i008p01667 .
Smits et al. (2010) Smits, K.M., T. Sakaki, A. Limsuwat, and T.H. Illangasekare. 2010. Thermal conductivity of sands under varying moisture and porosity in drainage–wetting cycles. Vadose Zone Journal 9(1): 172–180. 10.2136/vzj2009.0095 .
Spiteri et al. (2008) Spiteri, E.J., R. Juanes, M.J. Blunt, and F.M. Orr. 2008. A new model of trapping and relative permeability hysteresis for all wettability characteristics. SPE Journal 13(03): 277–288. 10.2118/96448-PA .
Spycher and Pruess (2005) Spycher, N. and K. Pruess. 2005. CO₂-H₂O mixtures in the geological sequestration of CO₂. II. Partitioning in chloride brines at 12–100 C and up to 600 bar. Geochimica et Cosmochimica Acta 69(13): 3309–3320. 10.1016/j.gca.2005.01.015 .
Spycher et al. (2003) Spycher, N., K. Pruess, and J. Ennis-King. 2003. CO₂-H₂O mixtures in the geological sequestration of CO₂. I. Assessment and calculation of mutual solubilities from 12 to 100 C and up to 600 bar. Geochimica et Cosmochimica Acta 67(16): 3015–3031. 10.1016/S0016-7037(03)00273-4 .
Szulczewski et al. (2013) Szulczewski, M., M. Hesse, and R. Juanes. 2013. Carbon dioxide dissolution in structural and stratigraphic traps. Journal of Fluid Mechanics 736: 287–315. 10.1017/jfm.2013.511 .
Thyng et al. (2016) Thyng, K.M., C.A. Greene, R.D. Hetland, H.M. Zimmerle, and S.F. DiMarco. 2016. True colors of oceanography: Guidelines for effective and accurate colormap selection. Oceanography 29(3): 9–13. 10.5670/oceanog.2016.66 .
Timur (1968) Timur, A. 1968. An investigation of permeability, porosity, & residual water saturation relationships for sandstone reservoirs. The Log Analyst 9(04): 8–17 .
Trevisan et al. (2014) Trevisan, L., A. Cihan, F. Fagerlund, E. Agartan, H. Mori, J.T. Birkholzer, Q. Zhou, and T.H. Illangasekare. 2014. Investigation of mechanisms of supercritical CO₂ trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions. International Journal of Greenhouse Gas Control 29: 35–49. 10.1016/j.ijggc.2014.07.012 .
Trevisan et al. (2017) Trevisan, L., R. Pini, A. Cihan, J.T. Birkholzer, Q. Zhou, A. González-Nicolás, and T.H. Illangasekare. 2017. Imaging and quantification of spreading and trapping of carbon dioxide in saline aquifers using meter-scale laboratory experiments. Water Resources Research 53(1): 485–502. https://doi.org/10.1002/2016WR019749 .
Verdon et al. (2013) Verdon, J.P., J.M. Kendall, A.L. Stork, R.A. Chadwick, D.J. White, and R.C. Bissell. 2013. Comparison of geomechanical deformation induced by megatonne-scale CO₂ storage at Sleipner, Weyburn, and In Salah. Proceedings of the National Academy of Sciences 110(30): E2762–E2771. 10.1073/pnas.1302156110 .
Vilarrasa and Carrera (2015) Vilarrasa, V. and J. Carrera. 2015. Geologic carbon storage is unlikely to trigger large earthquakes and reactivate faults through which CO₂ could leak. Proceedings of the National Academy of Sciences 112(19): 5938–5943. 10.1073/pnas.1413284112 .
Wang et al. (2010) Wang, S., M. Lee, M.K. Park, and J.M. Kim. 2010. Box experiments on monitoring the CO₂ migration in a homogeneous medium using electrical resistivity survey. Geosciences Journal 14: 77–85. 10.1007/s12303-010-0009-1 .
Weir et al. (1996) Weir, G.J., S.P. White, and W.M. Kissling. 1996. Reservoir storage and containment of greenhouse gases. Transport in Porous Media 23: 37–60. 10.1007/BF00145265 .
Wood et al. (1994) Wood, B.D., C.N. Dawson, J.E. Szecsody, and G.P. Streile. 1994. Modeling contaminant transport and biodegradation in a layered porous media system. Water resources research 30(6): 1833–1845. 10.1029/94WR00434 .
Zoback and Gorelick (2012) Zoback, M.D. and S.M. Gorelick. 2012. Earthquake triggering and large-scale geologic storage of carbon dioxide. Proceedings of the National Academy of Sciences 109(26): 10164–10168. 10.1073/pnas.1202473109 .

Direct Comparison of Numerical Simulations and Experiments of CO2 Injection and Migration in Geologic Media: Value of Local Data and Predictability

Abstract

keywords:

1 Introduction

2 Physical Experiments

2.1 Sand petrophysical properties

2.2 Experimental setup

2.3 Porous media geometries

2.4 Experimental injection schedule

3 Numerical simulations

3.1 Model setup

3.1.1 Computational grids

3.1.2 PVT properties

3.1.3 Initial, boundary and operational conditions

3.2 Simulation models

3.2.1 Model I

3.2.2 Model II

3.2.3 Model III

3.3 Model calibration

4 Results

4.1 Initial model results

4.2 Manual history matching and value of local data

4.3 Transferability: model predictions

4.3.1 Analog for a longer CO2 injection in the same geologic setting

4.3.2 Analog for a larger-scale CO2 injection in a different geologic setting

5 Discussion

6 Conclusions

Statements and Declarations

Appendix A Additional analysis of simulation model performance against experiment B1

Appendix B Comparison of simulation results with multiple grid resolutions

References

Direct Comparison of Numerical Simulations and Experiments of CO₂ Injection and Migration in Geologic Media: Value of Local Data and Predictability

4.3.1 Analog for a longer CO₂ injection in the same geologic setting

4.3.2 Analog for a larger-scale CO₂ injection in a different geologic setting