2023

[3,4]Martin Fernø

[5,4]Jan M. Nordbotten

[1,2]Ruben Juanes

[1]Department of Civil and Environmental Engineering, Massachusetts Institute
of Technology, Cambridge, 02139, MA, USA

2]Earth Resources Laboratory, Department of Earth, Atmospheric and Planetary
Sciences, Massachusetts Institute of Technology, Cambridge, 02139, MA, USA

3]Department of Physics and Technology, University of Bergen, Bergen, 5020,
Norway

4]Norwegian Research Center, Postboks 22 Nygårdstangen Bergen, 5020, Norway

5]Center for Modeling of Coupled Subsurface Dynamics, Department of
Mathematics, University of Bergen, Bergen, 5020, Norway

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

Lluís Saló-Salgado<EMAIL_ADDRESS>Malin Haugen<EMAIL_ADDRESS>Kristoffer
Eikehaug<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>*  [  [  [  [

###### Abstract

Purpose: The accuracy and robustness of numerical models of geologic CO2
sequestration are almost never quantified with respect to direct observations
that provide a ground truth. Here, we conduct CO2 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 CO2 injection experiment in tank 1
(89.7$\times$47$\times$1.05 cm). Matching is based on a quantitative
comparison of CO2 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 CO2 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 CO2 migration are
unlikely.

###### keywords:

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

##  1 Introduction

CO2 capture and subsequent geologic carbon sequestration (GCS) is a climate-
change mitigation technology that can be deployed at scale to oﬀset
anthropogenic CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2
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.

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 CO2 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. 1.

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

  2. 2.

The FluidFlower rig is filled with deionized water.

  3. 3.

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

  4. 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 CO2 gas (white),
dissolved CO2 (orange to red), and pure water (blue).

  5. 5.

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

  6. 6.

After the injection phase, injection ports are closed and CO2 migration
continues.

  7. 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 CO2 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 CO2 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 (CO2 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 CO2
flux (Bear, 1972).

The simulator requires minute timesteps due to the buoyancy of CO2 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).

Figure 2: Simulation grids overview. a front panel view of tank 1, where the
layer contacts have been highlighted in white. b front view of simulation grid
for experiments in tank 1, with lithologies indicated and colored based on
petrophysical properties (see sect. 3.2). Location of injection wells is shown
in red. c thickness map of simulation grid for experiments in tank 2. d front
view of simulation grid for experiments in tank 2, with lithologies indicated
and colored based on petrophysical properties. Location of injection wells is
shown in red.

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 CO2 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 CO2 have,
respectively, a density of 997 and 1.78 kg/m3, and a viscosity of 0.9 and
0.015 cP. The maximum concentration of CO2 in water is $\approx 1.4$ kg/m3.

####  3.1.3 Initial, boundary and operational conditions

Our porous media are fully saturated in water at the beginning of CO2
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 CO2 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}$ m2/s. Additionally, model III was also calibrated
with $D=3\times 10^{-9}$ m2/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. 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-CO2-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. 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. 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).

Figure 3: Multiphase flow properties for model I. a Capillary pressure
measurements and reference curve using a Brooks and Corey (1964) function. b
Initial capillary pressure curves, computed from the reference curve using
Leverett scaling (see main text). c Relative permeability data (squares and
$S_{\text{w}}^{5}$ model) and our fitted Corey model. d,e Relative
permeability of gas and water, respectively. The drainage curve is shown as a
solid line, while the bounding imbibition curve is shown for sands ESF and G
as a discontinuous line. No relative permeability hysteresis was considered
for the water phase.

Drainage relative permeabilities were obtained from CO2-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 CO2 storage, secondary imbibition occurs where the water displaces buoyant
gas at the trailing edge of the CO2 plume, disconnecting part of the CO2 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 CO2 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.

Figure 4: Multiphase flow properties for model II. b Initial capillary
pressure curves, computed from the reference curve using Leverett scaling (see
main text). b,c Relative permeability of gas and water, respectively. The
drainage curve is solid, while the bounding imbibition curve is shown for
sands ESF and G as a discontinuous line. No relative permeability hysteresis
was considered for the water phase.

####  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.

Figure 5: Multiphase flow properties for model III. b Initial capillary
pressure curves, computed according to the entry pressure determined
experimentally (see sect. 2.1). b,c Relative permeability of gas and water,
respectively, according to the endpoints determined experimentally (sect.
2.1). The drainage curves are solid, while the bounding imbibition curves are
shown as a discontinuous line. The inset in b is a zoom view around the
trapped gas saturation. No relative permeability hysteresis was considered for
the water phase.

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. 1.

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

  2. 2.

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

  3. 3.

Time at which the first finger touches the tank bottom.

  4. 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 CO2
(white) and dyed water with dissolved CO2 (yellowish orange to red), and to
identify the times for points 3-4 above. In the simulation models, the
threshold gas saturation and CO2 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/m3],
respectively. The $C$ value was chosen after a shape comparison of the region
with dissolved CO2. A smaller value of $C_{\text{CO}_{2}}>0.05$ [kg/m3] 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.

Figure 6: Front panel view of tank 1, showing quantities and times for history
matching of numerical models to experiment A1. a shows areas with gaseous CO2
(free-phase) and dyed water with dissolved CO2 at the end of injection.
Location of injection ports is shown with a star. b shows the time and
location where the first finger touches the bottom of the tank (white arrow),
as well as the different lithological units. Note the three F reservoirs
labeled ‘inf’, ‘mid’ and ‘sup’, mentioned in the text and other figures.

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 CO2 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
CO2-rich brine and underestimate the amount of gaseous CO2 in all F reservoirs
(refer to Fig. 6 for location). Model III approximates much better the areal
extent of gaseous CO2 in all regions, as well as the CO2-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.

Figure 7: Comparison between experiment A1 in tank 1 (left column) and first
run simulation results with models I-III. Color map in simulation plots refers
to CO2 concentration in water, according to color bar. The dotted red contours
in simulation plots indicate $S_{\text{g}}=10^{-3}$. a-d: end of injection in
port 1. e-h: end of injection in port 2. i-l: time at which the first finger
touches the tank bottom. m-p: time at which the first finger touches the
middle C sand.

###  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 CO2 ($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]$ cm2), 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/m3]
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/m3].

Figure 8: Convergence of areas occupied by free gas ($A_{\text{g}}$, left
column) and water with dissolved CO2 ($A_{\text{d}}$, right column), during
the calibration of models I-III with experiment A1. $A_{\text{d}}$ includes
area with gaseous CO2 (see Fig. 6). Each iteration represents a new simulation
run, and the experimental reference ($E$) is shown as a black line. Refer to
Fig. 6 for region location, and to sect. 3.3 for calibration procedure. a,b:
upper F sand. c,d: middle F sand. e,f: lower F sand. Figure 9: Convergence of
mean absolute error over the six areal quantities measured during the
calibration process. The error is computed with respect to experimental
values. See Fig. 8 for areas measured, and refer to sect. 3.3 for calibration
procedure.

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).

Figure 10: Convergence of times at which the first finger touches the bottom
of the rig (a) and the middle C sand (b), during the calibration of models
I-III with experiment A1. Refer to sect. 3.3 for calibration procedure.

The calibration shown in Fig. 8, 9, 10 employs $D=10^{-9}$ m2/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 CO2 migration was most
sensitive to the properties of the F sand, were most of the CO2 migration
occurs, as well as the ESF seal, which structurally traps the CO2 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}$ m2/s and m${}_{\text{III},3}$ refers to
$D=3\times 10^{-9}$ m2/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 CO2 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 CO2
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, CO2-rich water sinks with
a rather compact front and closely spaced, wide fingers. Our model with
constant $D=10^{-9}$ m2/s approximates all gravity-driven migration of the
CO2-rich water through thinner fingers, with the CO2-saturated region receding
with $S_{\text{g}}$. To better represent fingering widths, we also matched
model III with $D=3\times 10^{-9}$ m2/s, used in sect. 4.3.2.

Figure 11: Comparison between experiment A1 in tank 1 (left column) and
simulation results with model III after calibration (gas saturation shown in
middle column, and CO2 concentration shown in right column). Location of
injection ports shown by black stars in d. $D=10^{-9}$ m2/s. a-c: End of
injection in lower port. A contour at $C_{\text{CO}_{2}}=0.15\times
C_{\text{CO}_{2}}^{\text{max}}\approx 0.21$ kg/m3, the value used to compute
areas with dissolved CO2, is shown as a white solid line in c. d-f: End of
injection in upper port. g-i: Time at which the first finger touches the rig
bottom. j-l: Time at which the first finger touches the middle C layer.
Colormaps used in the simulations by Thyng et al. (2016).

###  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 CO2 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/s2, 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 CO2 is
close to the experimental truth. However, the experiment shows a compact
sinking front of the CO2-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, CO2-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}}$.

Figure 12: Comparison between experiment A2 in tank 1 (a) and simulation
results with model III (b,c) at the end of the injection phase ($t=4$h 48
min).

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

Figure 13: Comparison of areas occupied by free gas ($A_{\text{g}}$) and water
with dissolved CO2 ($A_{\text{d}}$) for experiment A2 in tank 1. Experimental
reference shown with a star ($E$). $A_{\text{g}}$ (F mid, left) not shown
because values are very close to 0. Refer to Fig. 6 or Fig. 12a for region
location.

####  4.3.2 Analog for a larger-scale CO2 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}$ m2/s) in all three models.
Additionally, a run with $D=3\times 10^{-9}$ m2/s was completed with model III
to better approximate finger widths, as noted in sect. 4.2.

Gas saturation and CO2 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
CO2 into box B. Model II (Fig. 15c) and III (Fig. 15d) underestimate the
amount of CO2, while model I (Fig. 15b) overestimates the amount of CO2 in the
top C sand.

  * •

Also at the end of injection, all models predict faster sinking of the
CO2-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}$ m2/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, CO2-rich water front closely trailing the fingers.

  * •

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

Figure 14: Comparison between experiment B1 in tank 2 (a) and simulation model
I (b,c) at the end of injection ($t=5$h). Circles in a denote the location of
injection ports. Figure 15: Comparison between experiment B1 in tank 2
(leftmost column) and CO2 concentration maps for simulation models I-III
(middle-left, middle-right and rightmost, respectively). $D=10^{-9}$ m2/s
(model I and II), $D=3\times 10^{-9}$ m2/s (model III). The red contours in
simulation plots indicate $S_{\text{g}}=10^{-3}$. a-d end of injection. e-h
$t=15$h. i-l $t=24$h. m-p $t=48$h. q-t $t=120$h.

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 CO2, and 2 is gaseous CO2. 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}$ m2/s continues to follow the experiment closely, whereas the others
overestimate gaseous CO2. 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}$ m2/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 CO2 migration into the seal (Fig. 15). In box B (Fig. 16d-f), model I and
model III with $D=10^{-9}$ m2/s are able to approximately track the
experimental truth during injection. However, our models without dispersion
cannot capture the areal increase of CO2-rich water that occurs afterwards
(cf. Fig. 15).

Figure 16: Comparison of areas occupied by each phase during the first 72h of
case B1. Experimental mean ($\overline{E}$) and standard deviation ($\sigma$)
obtained from four 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}$ m2/s
(m${}_{\text{III},1}$) and $D=3\times 10^{-9}$ m2/s (m${}_{\text{III},3}$).
Top row shows areas in box A, and bottom row shows areas in box B. a,d gaseous
CO2. b,e dissolved CO2 (includes area with gaseous CO2). c,f pure water.

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 CO2 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 CO2 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 CO2 injection and migration in
unconsolidated sands at atmospheric $p,\,T$ conditions capture the CO2-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 CO2 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}]$
m2/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 CO2-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 CO2 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 CO2 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 CO2-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 CO2 migration (see Flemisch et al.,
2023, this issue).

##  6 Conclusions

We performed experiments (sect. 2) and numerical simulations (sect. 3) of CO2
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. 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. 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 CO2 injection,
and accumulate larger errors afterwards, especially in regions where
heterogeneous structures control CO2 migration.

  3. 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, CO2-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 CO2
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 CO2
migration, particularly where heterogeneous structures are at play.

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 CO2 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 CO2 in the computational domain in Fig.
17, and the mass in boxes A and B in Fig. 18.

Figure 17: Total mass of CO2 for our simulations of experiment B1 presented in
sect. 4.3.2. Results are provided for models I to III. For m${}_{\text{III}}$,
two cases are shown: $D=10^{-9}$ m2/s (m${}_{\text{III},1}$) and $D=3\times
10^{-9}$ m2/s (m${}_{\text{III},3}$). Figure 18: Mass of CO2 in boxes A and B
defined in Fig. 1e, for our simulations of experiment B1 presented in sect.
4.3.2. Results are provided for models I to III. For m${}_{\text{III}}$, two
cases are shown: $D=10^{-9}$ m2/s (m${}_{\text{III},1}$) and $D=3\times
10^{-9}$ m2/s (m${}_{\text{III},3}$).

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. 1.

time of maximum mobile free phase in box A

  2. 2.

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

  3. 3.

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

  4. 4.

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

  5. 5.

total mass of CO2 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 CO2 (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 CO2 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}$ m2/s (m${}_{\text{III},1}$) and $D=3\times 10^{-9}$ m2/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. 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. 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 CO2 plume, (2) lower
dissolution, (3) lower number of fingers and finger widths, and (4) different
CO2-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 CO2 storage projects.

Figure 19: Concentration maps from our simulations of experiment A1 with model
III. Results with two grids are shown: $h=4$ mm (a, c, e, g) and $h=8$ mm (b,
d, f, h). Figure 20: Concentration maps from simulations of experiment B1 with
model III. Results with three grids are shown: $h=5$ mm (a, d, g, j, m),
$h=10$ mm (b, e, h, k, n) and $h=20$ mm (c, f, i, l, o).

## References

  *   * 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 CO2: 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 CO2.  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 CO2/brine systems and residual trapping of CO2.  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 CO2 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 CO2 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 CO2 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 CO2–brine mixtures for black-oil simulation of CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2: 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 CO2 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–CO2–water system under various pressure conditions. application to CO2 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.  CO2 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.  CO2-H2O mixtures in the geological sequestration of CO2. 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.  CO2-H2O mixtures in the geological sequestration of CO2. 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 CO2 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 CO2 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 CO2 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 CO2 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 .