University of Southern California Dissertations and Theses

Dynamics of CO₂ leakage through faults

(USC Thesis Other)

Dynamics of CO₂ leakage through faults

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content DYNAMICS OF CO
2
LEAKAGE THROUGH FAULTS
by
Saro Meguerdijian
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(PETROLEUM ENGINEERING)
December 2022
Copyright 2022 Saro Meguerdijian
To my family
ii
Acknowledgments
I am thankful to God for the many blessings I have enjoyed at USC during my Ph.D. studies, especially
the relationships I have had with those who contributed to my Ph.D. experience. Many have encouraged
me in my studies, and for this, I am very grateful.
I am grateful for the support and guidance provided by my advisor, Dr. Birendra Jha. His effort and
encouragement toward my growth as a scientist are much appreciated.
I would like to thank my dissertation committee, Dr. Birendra Jha, Dr. Iraj Ershaghi, Dr. Derek
Elsworth, Dr. Behnam Jafarpour, and Dr. Roger Ghanem, for taking the time to evaluate my defense
and dissertation.
I would also like to thank Dr. Aiichiro Nakano and Dr. Muhammad Sahimi for serving on my qualifying
exam committee.
I would like to acknowledge Dr. Fred Aminzadeh for serving on my screening exam committee. His
participation is appreciated.
Special thanks go to Dr. George V. Chilingar. His encouragement and advice many years ago helped
me to grow as a petroleum engineer.
Special thanks also go to Dr. Iraj Ershaghi. He has been a strong source of support throughout my
many years at USC, and his care for his students is legendary.
Dr. Donald G. Hill, under whom I worked as a teaching assistant, has followed up with me over the
years to encourage me to do my best. His enthusiasm and kind words are remembered.
iii
I am also grateful to have worked as a teaching assistant under Dr. Jincai Chang.
I would like to thank Dr. Joshua A. White from the Lawrence Livermore National Laboratory for his
mentorship and collaboration, both when I was at the Laboratory and afterward.
I would also like to thank Dr. Hari Viswanathan and Dr. Abigail Hunter for their mentorship during
my summer fellowship at the Los Alamos National Laboratory.
I would like to especially thank my collaborators from the Los Alamos National Laboratory with whom
I had the privilege to complete work toward my dissertation: my mentors Dr. Rajesh J. Pawar, Dr. Dylan
R. Harp, and Dr. Bailian Chen; and coauthors Dr. Carl W. Gable and Terry A. Miller. They helped me grow
as a scientist, too. I would also like to acknowledge Dr. Jha’s role in connecting me to them, an example
of his continual effort to seek the best for his students.
Many colleagues, peers, and friends contributed to my growth in my Ph.D. I express my gratitude to
them and hope to continue our mutually edifying friendships in the future.
I would like to acknowledge the support of the USC Provost Ph.D. Fellowship for my studies, including
the work done in Chapters 2, 3, and 4.
I would like to acknowledge the funding from the American Chemical Society Petroleum Research
Fund (PRF 58769-DNI9) for the study in Chapter 2.
The work in Chapters 3 and 4 was supported by the U.S. Department of Energy through the Los
Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security,
LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (Contract No.
89233218CNA000001). This work was completed as part of the National Risk Assessment Partnership
(NRAP) project. Support for this project came from the U.S. Department of Energy’s (DOE) Office of Fossil
Energy and Carbon Management’s Carbon Storage R&D program. (For clarity, during the time the study
in Chapter 3 was being conducted, the support was from the U.S. Department of Energy’s Office of Fossil
iv
Energy’s Coal Research program. The Office of Fossil Energy has since been renamed the Office of Fossil
Energy and Carbon Management.)
Finally, I would like to thank my family for their love and support.
v
TableofContents
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
TableofContents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
ListofTables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
ListofFigures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv
Chapter1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Background on Risk Assessment of GCS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Gaps in Knowledge and Scope of This Work . . . . . . . . . . . . . . . . . . . . . . . . . . 3
Chapter2: QuantificationofFaultLeakageDynamicsBasedonLeakageMagnitudeand
DipAngle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Mathematical and Numerical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.1 Leakage Magnitude Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.2 Physical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.3 Initial and Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3 Induced Fault Reactivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4 Induced Leakage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.5 Effect of Leakage and Fault Dip on Ground Deformation . . . . . . . . . . . . . . . . . . . 34
2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Chapter 3: Thermal and Solubility Effects on Fault Leakage during Geologic Carbon
Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2.1 Description of the Physics and Constitutive Equations . . . . . . . . . . . . . . . . 41
3.2.2 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.3 Results and Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3.1 Base Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3.2 Autocorrelation of Inputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.3.3 Solubility Effects on Brine and CO
2
Leakage Rates . . . . . . . . . . . . . . . . . . 56
3.3.4 Effects on Fault Reactivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
vi
3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.4.1 Solubility in the Fault-Leakage Context . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.4.2 Generalizability and Limitations for Real GCS Sites . . . . . . . . . . . . . . . . . . 65
3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
Chapter 4: Physics-Informed Machine Learning for Fault-Leakage Reduced-Order
Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.2.1 Description of the Mesh and Numerical Model . . . . . . . . . . . . . . . . . . . . 73
4.2.2 ROM Development Workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.3.1 Base Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.3.2 ROM Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.3.3 Sensitivity to Input Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.3.3.1 Delta Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.3.3.2 Sobol’ Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.3.3.3 PAWN Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.3.4 Prediction of Failed Runs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.4.1 Distribution of Parameters and Leakage Rates Among Complete vs. Incomplete
Reservoir Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.4.2 Physical Insights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.4.3 ROM Limitations and Applicability . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Chapter5: ConclusionsandFutureWork . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
AppendixA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
Supplemental Information for Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
vii
ListofTables
2.1 Input parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.1 Fixed parameters of the reservoir model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2 Varying parameters of the reservoir model. Square brackets denote a continuous, uniform
distribution while brackets denote a discrete distribution. *The x-axis is at the center of
the fault. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.3 The varying parameters for the base case. *The base case results are for the insoluble case
unless otherwise noted. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.1 The varying parameters of the model. Note that [] correspond to continuous ranges and
{} correspond to sets of discrete values.
∗ These correspond to distances from the well,
with larger values indicating larger distances. Because the meshes are unstructured, there
is some variation in the distance between dip angles. However, each well location value
corresponds roughly to a distance of 200 m from the intersection of the deep aquifer
with the fault.
†
These injection rates are for the half-reservoir due to the invocation of
symmetry along the injection well. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.2 Fixed parameters of the reservoir model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.3 The varying parameters of the base case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
viii
ListofFigures
2.1 (Left) A sketch of the faults and reservoirs in St. Johns Dome. CO
2
migrates from
magma across the Coyote Wash Fault (CWF) to the reservoir and the surface, where
travertines are formed. GWC indicates the gas-water contact. This is Figure 1(b) in Miocic
et al., 2019, used here with permission under the Creative Commons Attribution 4.0
International License CC, 2019. Text modified for readability. (Right) An example of how
a simultaneous, symmetric injection could occur in a real case of CO
2
sequestration. The
geological setting here is based on a CO
2
sequestration study in the Snøvhit field (Hansen
et al., 2013). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 A physical model of gas injection and leakage in faulted reservoirs. (a) Gas is injected
in the hanging wall block of a normal fault, which slips and begins to leak fluids into
aquifers A1 and A3 depending on the slip direction and gas buoyancy. (b) A conjugate
model similar to the model in (a) except that the fault dips in the opposite direction such
that the injection occurs in the footwall block. . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3 Pressure-dependent density and viscosity of gas and water phases. . . . . . . . . . . . . . . 22
2.4 Effect of the dip angle θ on overpressure, gas saturation and displacement after 149 days
of gas injection. Induced slip and subsequent leakage are shown. . . . . . . . . . . . . . . . 23
2.5 CFF profiles at three time steps showing the effect of the dip angle θ on the natural
selection of the slip nucleation points. The horizontal lines denote the boundaries of A2,
the injection aquifer. Slip begins at the depths where the CFF reaches zero. The nucleation
time, indicated by the timestamp of the dash-dot curve, is different for different θ . After
slip (solid black curve), CFF drops at the nucleation point due to slip-induced stress
relaxation. The x-axis scale is different for the four figures. . . . . . . . . . . . . . . . . . . 24
2.6 Fault dipθ affects the onset times of injection-induced slip (a) and leakage (b). (a) A2 Bot
and A2 Top indicate the hypocenter location in the hanging wall injection and footwall
injection cases, respectively. (b) A1 and A3 indicate the leakage location for the fourθ cases. Leakage of water is denoted by× symbols and leakage of gas by◦ symbols. (c)
The◦ symbol indicates the final state (last timestep) of the simulation. The lines show the
evolution of the leakage magnitude and seismic moment over time. M
0
andM
leak,g
have
units of N-m per unit length in the third dimension. (d) The dashed lines connecting the
circles are visual guides only. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
ix
2.7 Depth profiles of the fault overpressure, fault gas saturation (showing gas leakage into
aquifer), and slip along the fault are shown at three successive time steps (darker curves
indicate later time steps while lighter curves indicate earlier time steps). At early times,
the fault slips downdip in the hanging wall injection cases (θ = 60
◦ ,80
◦ ) and updip in
the footwall injection cases (θ = 100
◦ ,120
◦ ). Depths of the three aquifers–A1, A2, and
A3–are marked with black horizontal lines in each row. Up-dip slip is positive. . . . . . . . 26
2.8 Dynamics of fault leakage in the hanging wall injection (HWI) and footwall injection
(FWI) models in terms of evolution in fault permeability, fault slip (up-dip is positive),
and overpressure and gas saturation in aquifers A1, A2 and A3. The green vertical lines
denote the onset of leakage: t
leak
=125 days (HWI) andt
leak
=99 days (FWI). . . . . . . . . 29
2.9 A comparison of the evolution of the leakage magnitude (M
leak
) and the leakage mass
(m
leak
) as a function of the dip angle. The leakage magnitude and leakage mass have units
of N-m and kg, respectively, per unit length in the third dimension. . . . . . . . . . . . . . 31
2.10 The displacements at the upper boundary vary non-monotonically with time. Effect of
the fault dipθ on this evolution is evident. . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.11 The vertical displacements at the top boundary of the domain near the well’s horizontal
location are presented for each dip angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.1 A schematic of the fault-leakage model problem. . . . . . . . . . . . . . . . . . . . . . . . 44
3.2 The location of the corner cutout and permeable portion of the fault in the base case
figures. The corner cutout and permeable fault location are shown in solid colors. The
entire simulation domain is shown using a wireframe drawing. . . . . . . . . . . . . . . . . 52
3.3 Corner cutout views of the water saturation field along the fault in the base case at the
end of injection, where (a) does not include solubility effects and (b) includes solubility
effects. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.4 Corner cutout of CO
2
saturations along the fault in the base case at the end of injection,
where (a) shows the scCO
2
saturation without solubility effects, (b) shows the scCO
2
saturation with solubility effects included, (c) shows the gaseous CO
2
saturation without
solubility effects, (d) shows the gaseous CO
2
saturation with solubility effects included,
and (e) shows the dissolved CO
2
mass fraction with solubility effects included. . . . . . . 53
3.5 Corner cutout of pressure evolution in the fault in the base case, where (a) is the initial
pressure field and (b) is the pressure field at the end of injection. The permeable length of
the fault is evident in (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.6 Corner cutout of temperature evolution in the fault in the base case, where (a) is the
initial geothermal gradient, (b) is the temperature field at the end of the simulation, (c) is
difference in temperature in the fault between the end of the simulation and the initial
state, and (d) is the estimated thermal stress resulting from the temperature difference.
Note that normal stress is positive in tension. . . . . . . . . . . . . . . . . . . . . . . . . . . 54
x
3.7 The CO
2
phase state as a function of depth in the fault under initial conditions. . . . . . . 55
3.8 (a)∆ CFF and (b) pressure change from the initial pressure in the fault in the base case at
the end of injection. Note that positive values of∆ CFF indicate fault destabilization. . . . 55
3.9 Cumulative distribution functions of the brine and CO
2
maximum end-year leakage rates
(the maximum end-year leakage rate out of the 30 years of each simulation). The lower
bound for denoting a zero-leakage case for CO
2
is 10
− 4
t/yr. . . . . . . . . . . . . . . . . . 57
3.10 All changes shown are deviations from initial conditions. (a) Ensemble average of the
maximum pressure adjacent to the fault. Pressures taken at the time of the max(∆ CFF) for
each run in the ensemble. (b) Standard deviation of the ensemble pressures corresponding
to (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.11 A histogram and kernel-density plot of the maximum∆ CFF values encountered in each
run, both for runs with CO
2
solubility and those without CO
2
solubility. Note that the
soluble CO
2
histogram overlays the insoluble CO
2
histogram and a third color appears
where the two overlap. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.12 The locations of the points in the fault with the maximum ∆ CFF change for each run.
Note that while the runs including solubility and not including solubility were plotted
separately, the points overlap. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.13 The fraction of the permeable area in the fault with∆ CFF> 0.1 bar. . . . . . . . . . . . . 61
3.14 All changes shown are deviations from initial conditions. (a) Ensemble average of
the average temperature adjacent to the fault. Temperatures taken at the time of the
max(∆ CFF). (b) Standard deviation corresponding to (a). (c) Ensemble average temperature
increase in the fault at the end of the simulation. (d) Standard deviation of the ensemble
temperatures in (c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.15 All changes shown are deviations from initial conditions. (a) Ensemble average thermal
stresses corresponding to Fig. 3.14(a). (b) Standard deviation of the ensemble thermal
stresses in (a). (c) Ensemble average thermal stresses corresponding to Fig. 3.14(c). (d)
Standard deviation corresponding to (c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.16 Comparison of thermal stress and overpressure at the top and bottom of the fault. (a) The
thermal stress (positive in tension) and the overpressure in the fault relative to the initial
conditions at the top of the fault. (b) The thermal stress and the overpressure in the fault
relative to the initial conditions at the bottom of the fault but slightly above the storage
(deep) aquifer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.1 Physical setting of fault leakage. This figure is a modified version of Figure 1(b) from (Mio-
cic et al., 2019). Used under the Creative Commons Attribution 4.0 International License.
To view a copy of this license, please visit http://creativecommons.org/licenses/by/4.0/. . . 74
4.2 Schematic of the model problem. “BC” denotes boundary condition. . . . . . . . . . . . . . 74
xi
4.3 An example of the unstructured meshes used in these simulations. The dip angle (angle
from the horizontal) is 60
◦ . Note the additional refinement at the tops of the two aquifers
near the fault to more accurately capture CO
2
plume movement and the high refinement
of the fault to model flow accurately. Note also that the zoom-in figure of the fault has a
different color scheme to differentiate the fault from the surrounding caprock. . . . . . . . 75
4.4 Schematic of the deep neural network used in this study. . . . . . . . . . . . . . . . . . . . 76
4.5 Workflow for the generation of fault-leakage reduced-order models. . . . . . . . . . . . . . 77
4.6 The permeability and porosity distributions of the base case. Note that the fault appears
as a thin line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.7 The overpressure distribution of the base case at year 32. . . . . . . . . . . . . . . . . . . . 80
4.8 The temperature distribution of the base case at year 32. . . . . . . . . . . . . . . . . . . . 80
4.9 The base case’s (a) supercritical CO
2
saturation, and (b) water saturation at year 32 in the
simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.10 The performance of the ROM for predicting (a) instantaneous CO
2
leakage rate on the
training set, (b) instantaneous CO
2
leakage rate on the test set, (c) instantaneous brine
leakage rate on the training set, and (d) instantaneous brine leakage rate on the test set. . . 84
4.11 The prediction error distributions of the instantaneous CO
2
leakage rate for the test
set and the instantaneous brine leakage rate for the test set. The prediction error is the
predicted value of the leakage rate minus the actual value from the simulation. . . . . . . . 85
4.12 The performance of the ROM for predicting the following at 100 years: (a) cumulative
CO
2
leakage on the test set, and (b) cumulative brine leakage on the test set. . . . . . . . . 86
4.13 The performance of the ROM for predicting the onset time of CO
2
leakage. Note that the
threshold for determining the onset of CO
2
leakage is 275 t/yr. . . . . . . . . . . . . . . . . 87
4.14 The performance of the ROM for predicting (a) the maximum CO
2
leakage rate and (b)
the maximum brine leakage rate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.15 The performance of the ROM for predicting the maximum brine outflow rate from the
shallow aquifer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.16 The train and test errors of the ROM during training for the (a) CO
2
model, and (b) brine
model. Note that this training was performed excluding samples with time att=0 years
because the leakage value att=0 is set to 0.0 t/yr. . . . . . . . . . . . . . . . . . . . . . . . . 89
xii
4.17 Delta moment-independent measure values for (a) CO
2
-leakage DNN ROM and FEHM
simulations and (b) brine-leakage DNN ROM and FEHM simulations. The significant
variables for the ROM and FEHM simulations are those above their respective dummy
lines. The small black lines at the tops of the columns denote the 95% confidence intervals
for theδ values obtained from bootstrapping. . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.18 A Sobol’ sensitivity analysis of the ROM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.19 PAWN values for (a) CO
2
-leakage ROM and FEHM simulations and (b) brine-leakage ROM
and FEHM simulations. The significant variables for the ROM and FEHM simulations are
those above their respective dummy lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.20 Predicted vs. actual leakage rates for the available data from incomplete simulations for
(a) CO
2
and (b) brine. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.21 Stochastic regression-imputation estimates of R
2
of the ROM’s performance over missing
data for (a) CO
2
and (b) brine. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4.22 The distributional differences between the complete and incomplete simulations’
parameters using the test set of CO
2
and the set of incomplete runs. Note that due to the
use of a random selection of the CO
2
simulations’ test set, these parameter distributions
hold for the entire completed simulations’ dataset. . . . . . . . . . . . . . . . . . . . . . . . 102
S.1 Cumulative distribution functions of the input parameters for simulations with and
without solubility effects. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
S.2 Pearson correlation matrix of input parameters among simulations with solubility. . . . . . 126
S.3 Pearson correlation matrix of input parameters among simulations without solubility. . . 127
S.4 Spearman correlation matrix of input parameters among simulations with solubility. . . . 128
S.5 Spearman correlation matrix of input parameters among simulations without solubility. . 129
xiii
Abstract
Climate change due to rising anthropogenic greenhouse gas emissions has become a major environ-
mental issue over the past several decades. Geologic carbon storage is a promising technology to mitigate
these emissions by sequestering CO
2
in the subsurface and away from the atmosphere for long periods of
time. Geologic carbon storage requires the integrity of the subsurface storage system. Potential leakage
pathways, such as faults, may compromise this integrity, allowing CO
2
to migrate out of its storage reser-
voir to other locations in the subsurface and possibly even the atmosphere. Other risks also exist, such as
seismicity and groundwater impacts. Assessing the risk of a geologic carbon storage operation requires
understanding the physical phenomena involved, the leakage pathways which may exist, and the addi-
tional risks such as induced seismicity. These risks should be assessed in an integrated manner, which is
commonly done using full-physics simulations. However, these simulations can be computationally costly
to run in an ensemble, which is needed to quantify risk. In this dissertation, we seek to extend the physical
understanding of the dynamics of CO
2
leakage through faults and produce reduced-order models (ROM)
to reduce the computational cost of risk assessment. Over the course of three studies, we find that: the
fault dip angle is an important factor in the leakage behavior, yielding behavior such as downward gas
leakage and surface deformation effects that are dip-angle sensitive; thermal stresses during CO
2
leakage
can be sufficient to cause fault destabilization if CO
2
phase change is present; CO
2
solubility is impor-
tant for estimating CO
2
leakage rates but does not appear to be significant for fault stability estimates or
brine leakage-rate estimation; a deep learning approach yields a high-accuracy fault-leakage ROM; and
xiv
ROM performance can be evaluated over parameter spaces that have neither data nor simulation results
available. Future work includes extending the ROM with parameter ranges informed by field data and
incorporating fracturing, failure mechanisms, and CO
2
solubility.
xv
Chapter1
Introduction
1.1 Motivation
Climate change has grown over the past several decades from a scientific concern to a global societal
concern due to the substantial growth of anthropogenic greenhouse gas (GHG) emissions, particularly
CO
2
. One of the promising technologies to mitigate the rise in CO
2
levels in the atmosphere is carbon
capture, utilization, and storage (Hepburn et al., 2019; IPCC, 2005), notably geologic carbon storage (GCS),
which can potentially hold CO
2
underground and therefore separate it from the atmosphere for thousands
of years. The working principle of GCS is storing CO
2
underground so that it is not in the atmosphere and
therefore cannot influence the climate. GCS has particular appeal in that it allows for decoupling between
fossil fuel consumption and CO
2
emissions (cf. Chapters 2 and 5, Sinn, 2012).
Given the various challenges encountered in the energy transition to lower GHG emission sources,
GCS has long been considered a viable possibility for climate change mitigation (Bach, 1979) and has clear
and compelling precedents in nature with the existence, over geologic time, of natural gas reservoirs and
CO
2
reservoirs (such as St. Johns Dome, Miocic et al., 2019). With the expansion of the §45Q tax credits
in the United States (“26 USC 45Q: Credit for carbon oxide sequestration”, 2019), investment in GCS is
expected to increase substantially over the coming years. With the promise of decoupling GHG emissions
from energy consumption, clear policy incentives, and inspiration from nature, GCS seems to be near a
1
phase of rapid growth and positive impact on the environment, but these hopes can only be realized if
GCS is done safely and securely. Consequently, GCS requires risk assessment.
1.2 BackgroundonRiskAssessmentofGCS
GCS can involve certain risks, such as seismicity (White & Foxall, 2016), groundwater impacts (Keating
et al., 2010), leakage (Newell & Ilgen, 2019), and additional effects, necessitating a comprehensive risk
assessment strategy. The required components of such a strategy include:
1. Understanding the key subsurface physics and their couplings
2. Gathering data from experiments and field operations
3. Integrating this data and physics in a comprehensive model, such as a full-physics reservoir model
4. Running an ensemble of scenarios to quantify the uncertainties of the risks
Further background is provided on each of these points in this dissertation’s chapters as it is relevant to
the particular chapter. In this introduction, special focus will be given to the subject of the risk assessment
of fault leakage.
Fault leakage and its effects have been studied previously in the context of natural cases, but studies
with respect to fault leakage of CO
2
as a consequence of human activities are limited. Gas leakage in the
context of natural gas storage and naturally-occurring CO
2
leakage through faults has been studied (Khi-
lyuk et al., 2000; Miocic et al., 2019). Natural gas storage in underground reservoirs - a close analog of
GCS - was known to carry risks of fracturing and surface leakage, and in the case of the Leroy Gas Storage
Field in Wyoming, natural gas leakage through a fault was observed (Chapter 20, Khilyuk et al., 2000).
Consequently, the leakage of gas through faults is a known potential risk in gas storage operations. Case
studies on naturally-occurring CO
2
leakage through faults, such as Miocic et al., 2019, yield important
2
information about the nature of fault leakage but do not include elements essential to human-induced gas
storage, particularly wellbores and (in comparison to geologic time) fast fluid flow dynamics. A recent
field study has included both the CO
2
and gas storage aspects in the context of studying fault leakage but
in a shallow context (“The CSIRO In-Situ Laboratory," n.d.), while expected field applications are expected
to be significantly deeper (Newell & Ilgen, 2019). It is notable that the topic of GCS lies between natural
gas storage and naturally-occurring CO
2
reservoirs but is fully captured by neither. Data on CO
2
leakage
through faults is scarce because it is generally not the objective of GCS operations. Consequently, gaps in
knowledge regarding the risk assessment of GCS exist, which this dissertation aims to address partially.
1.3 GapsinKnowledgeandScopeofThisWork
The dynamics of CO
2
leakage through faults during GCS are increasingly important to understand
and quantify to ensure safe GCS operations, but knowledge of particular aspects of this leakage is needed.
These gaps in knowledge include the effect of fault dip angle on fault stability and CO
2
leakage behavior,
the consequences of thermal stresses induced by CO
2
leakage in faults, the importance of CO
2
solubility in
leakage and fault destabilization, and reduced-order models (ROM) for reducing the computational burden
of simulating fault leakage to make running an ensemble of scenarios for risk assessment computationally
viable. (In the CO
2
storage and fault-leakage literature, the term “reduced-order model” is used for both
a reduced-physics model and what in other fields is called a “statistical model.” This dissertation follows
the CO
2
storage and fault-leakage literature’s conventions in this naming but notes that other fields may
use different nomenclature. For clarity, the ROM developed in Chapter 4 is a deep neural network, not
a reduced-physics model.) The impetus for investigating the effect of fault dip angle is that the stability
of faults in a particular stress configuration changes with its dip angle, potentially making the fault more
likely to slip (Scholz, 2018). Slip can then potentially increase permeability (Dempsey et al., 2013), causing
3
further leakage. As faults in nature have a variety of dip angles (Scholz, 2018), it is necessary to take this
into consideration for improved GCS site selection.
Thermal stresses in the fault and CO
2
solubility effects could also exert influences on fault stability and
leakage. Thermal cooling around injection wells during CO
2
injection can cause stress changes that may
affect the stability of distant faults (Vilarrasa & Rutqvist, 2017). However, it is not certain if leaking CO
2
would cause fracturing, as the CO
2
, warmed by deeper formations, may actually impart heat to the fault as
it leaks, promoting compressive thermal stresses and consequently fault stabilization. Given the increased
computational difficulty with adding additional physical couplings to a simulator, knowing the magnitude
of thermal effects vs. other effects (such as poroelasticity) would help in guiding simulation studies and,
ultimately, ROM development. CO
2
solubility can also be difficult to incorporate in a simulator due to
its additional computational costs, but it is clearly important in the case of groundwater impacts (Keating
et al., 2010). Still, it is not clear if CO
2
solubility exerts a substantial impact on fault stability or on brine or
CO
2
leakage rates, effects which are important to know when planning a GCS operation. These physics,
therefore, require additional study.
There is also an effort to make ROMs for fault leakage to reduce the computational cost of evaluating
scenarios, and making these ROMs requires improved methods over previous ROM attempts. ROMs have
been made for several components of the GCS system (D. R. Harp et al., 2016; Pawar et al., 2016; Stauffer
et al., 2020; Vasylkivska et al., 2022; Vasylkivska et al., 2021). However, fault-leakage ROMs are notably
scarce: only a few models have been proposed over the last decade (Lindner, 2020; Lu et al., 2012; White
et al., 2017). One is readily usable for testing purposes as part of NRAP-Open-IAM, a U.S. DOE GCS risk
assessment software (Lindner, 2020; Vasylkivska et al., 2022). Fault-leakage ROMs also need to be validated
over their range of potential use, but field data is largely absent and even simulations may not be able to
simulate certain parameter spaces due to fast fluid flow physics (cf. Section 3.1, D. R. Harp et al., 2016).
4
There is a need to be able to comprehensively evaluate ROM accuracy even in the absence of data or
simulator results to build confidence in the ROM’s results for GCS risk assessment.
This dissertation makes progress in addressing these knowledge gaps through the studies described in
Chapters 2, 3, and 4. Chapter 2 describes a study examining the effect of fault dip angle on the leakage of
a generic gas (e.g., natural gas, CO
2
) and finds that while a relationship could not be identified between
seismicity and leakage, fault dip angle does impact the onset of leakage and even the direction of leakage,
with downward leakage being a possibility. Chapter 3 finds that thermal effects during fault leakage are
important if CO
2
passes a phase-change boundary, resulting in substantial cooling and tensile stresses,
which could cause fracturing. In the absence of phase change, the study does not show significant ther-
mal stresses resulting from CO
2
leakage. Furthermore, the study finds that CO
2
solubility is important
for estimating CO
2
leakage rates but not brine leakage rates or geomechanical fault stability. Chapter 4
shows work done to build a fault-leakage ROM using deep learning, which is different from previous
ROMs involving reduced physics models (Lindner, 2020), polynomial models (Lu et al., 2012), or response
surfaces (White et al., 2017). The ROM is extensively characterized, and a first-of-its-kind application of
the multiple imputation workflow to ROM performance evaluation is performed. This involves estimating
ROM performance over parameter spaces for which simulator data is unavailable. The collective conclu-
sions and future work of these studies are discussed in Chapter 5.
These studies have been published or are under review for publication as follows (with minor correc-
tions and changes being made as part of this dissertation):
1. Chapter 2: Meguerdijian, S., & Jha, B. (2021). Quantification of fault leakage dynamics based on
leakage magnitude and dip angle. International Journal for Numerical and Analytical Methods in
Geomechanics, 45(16), 2303–2320. https://doi.org/10.1002/nag.3267
5
2. Chapter 3: Meguerdijian, S., Pawar, R. J., Harp, D. R., & Jha, B. (2022). Thermal and solubility effects
on fault leakage during geologic carbon storage. International Journal of Greenhouse Gas Control,
116, 103633. https://doi.org/10.1016/j.ijggc.2022.103633
3. Chapter 4: Meguerdijian, S., Pawar, R. J., Chen, B., Gable, C. W., Miller, T. A., & Jha B. (2022).
Physics-informed machine learning for fault-leakage reduced-order modeling. International Journal
of Greenhouse Gas Control. Under review.
6
Chapter2
QuantificationofFaultLeakageDynamicsBasedonLeakageMagnitude
andDipAngle
Abstract
Subsurface leakage of fluids, such as CO
2
, along faults is a key quantity to estimate when model-
ing underground fluid storage. Using our coupled multiphase flow-geomechanics-fault slip simulator, we
quantify leakage dynamics using the fault dip angle and leakage magnitude, a proposed metric of gas leak-
age. We present novel leakage dynamics of faults to show that leakage is non-trivially coupled to induced
seismicity and multiphase flow along faults due to the effect of fault dip. The onset time of induced fault
slip, controlled by the initial shear-to-effective normal stress ratio on the fault, is a non-monotonic func-
tion of the fault dip. The leakage directions of the gas and liquid phases are determined by the directions
of slip propagation and the buoyancy vector, both of which depend on the dip. A consequence is that leak-
age evolution is non-monotonic in time for hanging wall injection-induced seismicity on a normal fault
because of the competition between up-dip oriented buoyancy and down-dip oriented induced slip. With
respect to monitoring, we note that subsidence at an injection well could be indicative of leakage.
Keywords: fault leakage, induced seismicity, underground gas storage
7
2.1 Introduction
Injection and storage of fluids in underground aquifers are important for many industrial applications.
Natural gas is stored underground during the summer months and produced during the winter months to
meet the seasonal demand of heating gas (Jha et al., 2015; Teatini et al., 2011). Radioactive waste (Alexan-
der et al., 2015; Ghanbarzadeh et al., 2015) and wastewater (Goebel et al., 2016; Horton, 2012; Hosseini
et al., 2018; Keranen et al., 2013) are candidates for injection underground for long-term disposal. CO
2
is injected underground as part of Carbon Capture and Sequestration (CCS) and enhanced oil recovery
projects (Bachu, 2000; Birkholzer et al., 2015; Birkholzer & Zhou, 2009; IPCC, 2005; Klusman, 2003; Nicot
et al., 2013; White & Foxall, 2016; Zhao & Jha, 2019). Underground storage is also a naturally occurring
process e.g. during creation of petroleum reservoirs in structural traps (Downey, 1984). Injection-induced
stresses can cause ground deformation (Y. Jung et al., 2013; Morris et al., 2011; Rucci et al., 2013; Sun &
Nicot, 2012; Vasco et al., 2010; Verdon et al., 2013) and induce shear and tensile failures along faults and
fractures (Cesca et al., 2014; Fialko & Simons, 2000; Jha & Juanes, 2014; Mazzoldi et al., 2012; Rutqvist et al.,
2008; Streit & Hillis, 2004; Zoback & Gorelick, 2012), which are ubiquitous in the subsurface. Fault slip is
often associated with changes in permeability that can lead to flow of stored fluids up-dip or down-dip
along a fault and consequently leakage into shallower or deeper formations, respectively (Ghanbarzadeh
et al., 2015; Guglielmi et al., 2015; Lewicki et al., 2003; Pruess, 2005; Vilarrasa et al., 2016). Leakage can
cause contamination of underground sources of drinking water (Llewellyn et al., 2015; Osborn et al., 2011)
and pose health and environmental risks on the ground surface (R. Hepple & Benson, 2002). In case of CO
2
storage, unless below a leakage limit (R. P. Hepple & Benson, 2005), even a small leakage flux can result
into a significant loss of volume over the hundred-year time scale of such projects. Such concerns have fu-
eled research in monitoring of leakage via various means e.g. geochemical (Romanak et al., 2012; L. Zheng
et al., 2012), hydrological (Lewicki et al., 2005; Shakiba & Hosseini, 2016; Sun & Nicot, 2012; Zeidouni &
8
Pooladi-Darvish, 2012), geophysical (Harbert et al., 2016; Luth et al., 2015), electrical (Harbert et al., 2016;
Strazisar et al., 2009; X. Yang et al., 2015), and isotopic (Jones et al., 2014; Riding & Rochelle, 2009).
Selection of a storage site requires evaluation of existing and potential leakage risks while accounting
for the site geology and dynamic processes of multiphase flow, rock deformation and fault slip during
injection (Birkholzer et al., 2015; White & Foxall, 2016). These evaluations inform decisions on where to
place injection wells with respect to known faults, what the maximum injection pressure should be, how
much fluid can be stored, what type of leaks can occur and how, and what parameters to measure during
leakage monitoring. For reservoirs containing dip-slip faults, an important question is which fault block—
hanging wall or footwall block—and which faulting environment—normal or reverse—offers a lower risk of
leakage (Pereira et al., 2014; Vilarrasa et al., 2016). While knowingly injecting near a fault under conditions
that cause leakage and seismicity is unlikely for a gas storage operation, a comprehensive risk assessment
should consider the possibility and consequences of leakage and induced seismicity and account for them
in storage site selection and operational plans. This requires numerical methods and metrics capable of
addressing the complexity of such operations.
In this paper, we show the capability of a coupled multiphase flow-geomechanics simulator (Jha &
Juanes, 2014; Zhao & Jha, 2019; Zhao & Jha, 2021) to simulate a scenario involving large permeability and
pressure changes in a fault during slip and suggest a new metric, the leakage magnitude, which when cou-
pled with the dip angle and leakage mass can be used to analyze such scenarios. We consider a numerical
study of an underground gas injection operation in a faulted reservoir with a vertical offset between its
hanging wall and footwall block compartments. This case extends a fault leakage study regarding St. Johns
Dome reservoir across the Arizona-New Mexico border (Miocic et al., 2019) by considering: the coupling
between fault leakage and seismicity, and downward gas migration across a fault (such as the CWF fault
in Figure 2.1). This paper also extends the CO
2
sequestration study in the Snøhvit field area (Hansen et al.,
2013) by considering the effect of simultaneous injection by two wells into an aquifer located in a fault
9
ramp relay pattern (see Figure 2.1). Examining the effect of symmetric injection provides a step toward
better understanding the coupling between fault leakage and seismicity for general multi-well injection
into an aquifer.
Injection-induced fault slip can cause permeability enhancement and so potentially lead to leakage
of the injected gas along the fault, a phenomenon noted in other studies, albeit limited in the amount of
leakage (Rinaldi et al., 2014). We characterize the leakage dynamics, i.e. evolution in gas phase pressure
and saturation inside aquifers and a fault, and relate it to the fault slip dynamics, i.e. evolution of slip
and permeability of the fault. We study these characteristics for a range of fault dip angles to capture the
effects of fault block and faulting environment on injection-induced seismicity and leakage. For a given
injection location, fault angles (measured counterclockwise from the horizontal axis) within0
◦ to90
◦ and
within 90
◦ to 180
◦ range can indicate injection in two different fault blocks (hanging wall and footwall
blocks) of the same faulting environment, or it can indicate two different faulting environments (normal
and reverse) of the same fault block.
We find that leakage resulting from injection into the hanging wall block of a normal fault is delayed
and smaller in flux than the leakage from the footwall block injection due to the combined effect of gas
buoyancy and tectonic stress regime. A delayed and smaller leakage leads to larger magnitudes of over-
pressure and induced seismicity in the hanging wall injection case compared to the footwall injection case.
This is analogous to delayed and larger magnitude induced events reported around a horizontal well com-
pared to a vertical well (Rinaldi et al., 2015). We discuss the implication of this two-way coupling between
slip (seismic or aseismic) and leakage—slip inducing leakage, and leakage controlling the magnitude of
seismicity—on site selection and risk assessment for underground storage projects. We note the utility of
comparing the leakage magnitude with other metrics as a method of quantifying fault leakage dynamics.
10
St. Johns Dome
Springerville
Volcanic Field
CWF
Travertines
GWC today
Paleo
GWC
Buttes
Fault
Reservoir
W E
CO
2
Simulation Domain
Symmetric CO
2
Injection
CO
2
Figure 2.1: (Left) A sketch of the faults and reservoirs in St. Johns Dome. CO
2
migrates from magma across
the Coyote Wash Fault (CWF) to the reservoir and the surface, where travertines are formed. GWC indi-
cates the gas-water contact. This is Figure 1(b) in Miocic et al., 2019, used here with permission under the
Creative Commons Attribution 4.0 International License CC, 2019. Text modified for readability. (Right)
An example of how a simultaneous, symmetric injection could occur in a real case of CO
2
sequestration.
The geological setting here is based on a CO
2
sequestration study in the Snøvhit field (Hansen et al., 2013).
2.2 MathematicalandNumericalModel
We simulate injection-induced fault activation and leakage using our coupled multiphase flow and
geomechanics simulation framework (Jagalur-Mohan et al., 2018; Jha & Juanes, 2014; Panda et al., 2018;
Zhao & Jha, 2019; Zhao & Jha, 2021) that solves for gas pressure, saturation, displacement and stress in the
domain and the effective normal traction σ ′
n
(positive in tension), shear tractionτ (positive down-dip) and
slip (positive down-dip) on the fault. The numerical model is 2-D with the number of space dimensions as
n
dim
=2. We consider isothermal conditions.
The governing equations for solid momentum balance under quasi-static equilibrium are the following:
∇· σ +ρ b
g =0, (2.1)
whereσ is the Cauchy total stress tensor,g is the gravity vector, andρ b
=ϕ P
n
phase
β ρ β S
β +(1− ϕ )ρ s
is the bulk density,ρ β andS
β are the density and saturation of fluid phase β (=w,g for water and gas),ρ s
is the density of the solid phase,ϕ is the porosity,g is gravitational acceleration, andn
phase
is the number
11
of fluid phases. The porosity is defined as the ratio of the pore volume to the bulk volume in the original
configuration.
Assuming that the fluids are immiscible, the mass-conservation equation for each phase α (=w,g for
water and gas) is
dm
α dt
+∇· w
α =ρ α f
α , (2.2)
where the accumulation termdm
α /dt describes the time variation of fluid mass relative to the motion
of the solid skeleton,w
α is the phase mass-flux relative to the solid skeleton, f
α is the volumetric source
term, andρ α is the phase density. We assume a multiphase fluid system and use the incremental formula-
tion of poromechanics for multiphase systems (as opposed to linearizing stress from the initial to current
state) (Coussy, 1995). Increment in total stressδσ is given in terms of the increments in effective stress
δσ ′
and pore pressureδp as follows:
δσ =δσ ′
− X
β b
β δp
β 1, δσ ′
=C
dr
:δε (2.3)
whereC
dr
is the rank-4 drained elasticity tensor given in terms of Young’s modulusE and Poisson’s
ratioν for a homogeneous isotropic linear elastic material,1 is the rank-2 identity tensor,ε is the linearized
strain tensor, defined as the symmetric gradient of the displacement vector u, andb
β are the Biot coeffi-
cients for individual phases such that
P
β b
β = b, whereb is the Biot coefficient of the saturated porous
material. We assume thatb
β =S
β b following previous work (Coussy et al., 1998; Lewis & Schrefler, 1998;
Lewis & Sukirman, 1993). For a two-phase system (waterw and gasg), the expression for capillary pressure
becomes:
P
c
(S
w
)≡ P
wg
(S
w
)=p
g
− p
w
. (2.4)
12
We choose the van Genuchten capillary pressure model (van Genuchten, 1980) such that:
P
wg
=P
o
((S
n
)
− 1/m
− 1)
1− m
,
S
n
=(S
w
− S
wic
)/(1− S
wic
),
S
wic
=0.1S
wi
(2.5)
whereP
o
is the capillary entry pressure andS
wi
is the irreducible water saturation. The final mass balance
equation for fluid phase α is (Jha & Juanes, 2014):
∂
∂t


ρ α X
β
N
αβ +
b
α b
β K
dr

p
β 

+
1
K
dr
∂
∂t
(ρ α b
α σ v
)+∇· w
α =ρ α f
α , (2.6)
where K
dr
is the drained bulk modulus, σ v
= tr(σ )/n
dim
is the volumetric stress. Expressions of N
αβ are given elsewhere (Jha & Juanes, 2014; Lewis & Schrefler, 1998). The mass flux of phase α is defined as
w
α =ρ α v
α wherev
α is the Darcy velocity given by
v
α =− kk
r
α µ α (∇p
α − ρ α g), (2.7)
whereµ α andk
r
α are the dynamic viscosity and the relative permeability of phaseα . The relative perme-
abilities are modeled using Corey-type relations (Brooks & Corey, 1964): k
rw
= (
ˆ
S
w
)
p
, k
rg
= (1− ˆ
S
w
)
p
and
ˆ
S
w
=(S
w
− S
wi
)/(1− S
wi
− S
gr
).
Faults are represented as surfaces of displacement discontinuity within the domain (Zhao & Jha, 2021)
such that
(u
+
− u
− )=d on S
fault
, (2.8)
13
whereu
+
andu
− are the displacements on the two sides of the fault surface, denoted byS
fault
, andd is
the fault slip vector. Here, the fault slip vector has one component, d, which is slip along the fault. The
following fault constitutive model is used to compute the frictional stressτ f
along the fault:
τ f
=











τ c
− µ f
σ ′
n
, σ ′
n
<0,
τ c
, σ ′
n
≥ 0,
(2.9)
where τ c
is the cohesive strength of the fault and µ f
is the coefficient of friction. Note that stress is
considered positive in tension. The following slip-weakening model is used to calculateµ f
:
µ f
=











µ s
− (µ s
− µ d
)
|d|
dc
, |d|≤ d
c
,
µ d
, |d|>d
c
.
(2.10)
whereµ f
is a function of the slipd, and drops linearly from its static value,µ s
, to its dynamic value,µ d
,
over a critical slip distanced
c
. The Mohr-Coulomb theory is used to define the stability criterion for the
fault (Jaeger & Cook, 1979; Reasenberg & Simpson, 1992; Scholz, 2018):
CFF=|τ |+µ f
σ ′
n
(2.11)
where CFF is the Coulomb failure function and fault cohesion is assumed to be zero. Initialization of the
hydromechanical model ensures that CFF is negative initially everywhere on the fault. Later, injection-
induced overpressure and stress changes around some points on the fault causes the CFF to increase at
14
those points. Slip begins at a point when the CFF there exceeds zero. The permeability of the faultk
f
is
assumed to increase with shear slipd according to the following formula (Dempsey et al., 2013):
∆ k
f
=
∆ k
f,max
1+exp

ln(19)

1− 2
|d|− d
5
d
95
− d
5
(2.12)
where∆ k
f,max
is the maximum increase in permeability,d is the shear displacement,d
5
is shear displace-
ment at which 5% of the permeability increase occurs, andd
95
is the shear displacement at which 95% of the
permeability increase occurs. The form of (Eq. 2.12) reflects general features observed in experiments (Lee
& Cho, 2002) and can be adapted to different lithologies (Dempsey et al., 2013).
The weak form of the mechanics problem is as follows (Jha & Juanes, 2014):
Z
Ω ∇
s
η :(σ ′
− bp
E
1)dΩ+
Z
Γ f
+
η · (l− bp
f
n)dΓ − Z
Γ f
− η · (l− bp
f
n)dΓ − Z
Ω η · ρ b
gdΩ − Z
Γ σ η · tdΓ=0 .
(2.13)
Z
Γ f
+
η · u
+
dΓ − Z
Γ f
− η · u
− dΓ − Z
Γ f
η · ddΓ=0 , (2.14)
where Ω denotes the domain, η denotes the test functions belonging to the functional space satisfy-
ing η = 0 on Γ u
(the Dirichlet boundaries of the mechanics problem having fixed displacements),
p
f
=max(p
+
,p
− ) is the fault pressure (where p
+
and p
− are the pressures on the plus and minus sides
of the fault, respectively) as shown in grain-scale simulations of fault gouge poromechanics (Z. Yang &
Juanes, 2018),Γ denotes a boundary integration variable,Γ f
+
is the interface on the plus side of the fault
surface,Γ f
− is the interface on the minus side of the fault surface,Γ f
is the fault surface,Γ σ are Neumann
boundaries with prescribed tractions,t are the prescribed tractions,l are the Lagrange multipliers and also
the effective fault tractions, n is the outward unit normal to a given boundary, andp
E
=
P
β S
β p
β is the
equivalent pressure. The fault pressure definition above yields the most conservative estimate against the
15
fault reactivation hazard, because a maximum value of the pressure implies a larger CFF and hence a faster
destabilization.
The weak form of the flow problem is as follows for the water phase:
∂
∂t
Z
Ω i
ρ w

N
ww
+N
wg
+
bb
w
K
dr

p
g
−
N
ww
+
b
2
w
K
dr

P
wg

dΩ +
1
K
dr
∂
∂t
Z
Ω i
ρ w
b
w
σ v
dΩ − Z
∂Ω i
w
w
· n
i
dΓ=
Z
Ω i
ρ w
f
w
dΩ , (2.15)
where Ω i
is the region enclosed by element i within the domain and n
i
is the outward unit normal to
the boundary of element i. We have a similar equation for the gas phase. The spatial discretization of
the equations involves the use of the finite volume method for the flow problem and the bilinear finite
element method for the mechanics problem. The Backward Euler method is used to integrate the semi-
discrete equations in time. The details of the above formulations and their solution methods have been
discussed in detail previously (Jha & Juanes, 2014). We solve the two problems implicitly and sequentially
using the fixed stress method (Jha & Juanes, 2014; Zhao & Jha, 2019; Zhao & Jha, 2021).
Here, we discuss our extension of the previous formulation to include flow along and across discrete
fractures represented as lower dimensional objects (in this case as 1D segments) in the mesh. This ex-
tension is novel and allows us to model fault leakage simultaneously with fault slip and slip-dependent
permeability, a capability the previous formulation did not have.
We model the fault as a zero-thickness 1D segment with an associated fault aperture, prescribed as a
fault property, to scale the equations. We adapt previous work (Karimi-Fard et al., 2004; Karimi-Fard &
Firoozabadi, 2003) to our formulation. We calculate the volumetric flux between the fault and adjacent
reservoir mesh elements using two-point flux approximation as:
Q
AB,β
=T
AB
λ β (p
A,β
− p
B,β
) (2.16)
16
where p
A,β
is the pressure of phase β in Cell A, p
B,β
is the pressure of phase β in Fault Segment B,
Q
AB,β
is the flow rate of phase β from CellA to Fault SegmentB,λ β is the fluid mobility of phase β , and
T
AB
is the phase-independent geometric portion of the face transmissibility, whereT
AB
is computed from
harmonic mean of half-transmissibilities as:
T
AB
=
α A
α B
α A
+α B
, whereα i
=
A
i
k
i
D
i
n
i
· f
i
, (2.17)
whereA
i
is the area of the interface between the cell and fault segment;k
i
is the absolute permeability of
either the cell or fault segment;D
i
is the distance between the centroid of the cell or fault segment and the
centroid of the interface between the cell and fault segment;f
i
is the unit vector in the direction of the
line connecting the cell or fault segment centroid and the centroid of the interface between the cell and
the fault segment; andn
i
is the unit normal vector of the interface. Given the small thickness (aperture)
of the fault segment relative to the reservoir cell size,f
i
andn
i
are co-linear: n
i
· f
i
= 1. In the case
of a fault segment connecting to another fault segment, the situation is similar except thatA
i
is the fault
aperture (times the dimension in the third direction) for both fault segments.
2.2.1 LeakageMagnitudeDefinition
Analogous to the concept of seismicity quantification, we define location, time, and magnitude as three
quantities that can quantify a leakage event. We define them as follows: the leakage location is the offset
aquifer into which the injected fluid leaks, the leakage onset time, t
leak
, is the time when leakage begins,
and the leakage magnitude is a quantity related to the increase inenergy of the offset aquifer due to leakage.
We note that, in contrast to seismicity, leakage is a multiphase flow event for which a magnitude definition
does not exist in the literature. Here, we propose a definition of the leakage magnitude that is analogous
to the moment of an earthquake (Hanks & Kanamori, 1979). First, we recognize that all three quantities
associated with a leakage event–location, time, and magnitude–should be defined for a specific fluid phase.
17
Parameters Value
Young’s modulus,E (GPa) 10
Poisson’s ratioν (-) 0.25
Biot’s coefficient, b (-) 1
Aquifer permeability (m
2
) 1.97e-13
Non-aquifer permeability (m
2
) 9.87e-22
Initial fault permeability in aquifer zones (m
2
) 9.87e-15
Initial fault permeability outside aquifer zones (m
2
) 9.87e-20
Porosity,ϕ (-) 0.1
Relative permeability exponent,p (-) 2
Irreducible water saturation,S
wi
0.12
Residual gas saturation,S
gr
0
Van Genuchten shape parameter,m (-) 0.5
Van Genuchten entry pressure,P
0
(kPa) 13.8
Maximum permeability increase at fault,∆ k
f,max
(m
2
) 9.87e-8
5% of maximum permeability increase distance,d
5
(mm) 1
95% of maximum permeability increase distance,d
95
(mm) 10
Static coefficient of friction, µ s
0.5
Dynamic coefficient of friction, µ d
0.2
Critical slip distance,d
c
(m) 0.5
Fault cohesive strength, (kPa) 0
Fault aperture (m) 0.3
Fault porosity 0.9
Solid density,ρ s
(kg/m
3
) 2260
Fluid-to-bulk density ratio,r
den
0.47
Gas injection rate,f
g
(standard m
3
/day/m) 67,960
Standard density of injected gas,ρ g,0
(kg/m
3
) 1.794
Standard density of water,ρ w,0
(kg/m
3
) 1000
Horizontal-to-vertical total stress ratio,r
tec
(-) 0.7
Injection location,(x,z) (m) (378,1506)
Aquifer thickness (m) 100
Injection aquifer length (m) 500
No. of grid cells in the x-direction, N
x
52
No. of grid cells in the z-direction, N
z
58
Average grid cell size in the x-direction, dx (m) 38.5
Average grid cell size in the z-direction, dz (m) 34.5
Table 2.1: Input parameters.
18
For example, water can start to leak earlier, into a deeper offset aquifer, and create a smaller magnitude
leakage event than gas. In this paper, we focus on leakage of the injected gas phase into an offset aquifer.
Such an event is associated with increases in pressure and gas saturation of the offset aquifer. Therefore,
we define the leakage magnitude as
M
leak,g
(t)=
Z
V
off
ϕ (p
g
(x,z,t)− p
0
(x,z))s
g
(x,z,t)dV, (2.18)
whereV
off
is the bulk volume of the offset aquifer into which the gas is leaking, p
g
is the gas phase pressure,
p
0
is the initial pressure, ands
g
is the leaked gas saturation in the offset aquifer. M
leak,g
quantifies the size
of a leakage event in analogy with the seismic momentM
0
(Hanks & Kanamori, 1979), which quantifies
the size of a seismic event. Here, M
0
(t) =
R
S
fault
G(x,z)|d(x,z,t)| dS, where S
fault
is the fault surface
area, andG is the shear modulus. Note that the leakage magnitude is a complementary numerical measure
to the leakage amount: while the leakage amount measures the mass of leakage, the leakage magnitude
measures the propensity of leakage by taking into account both the amount of leakage and its pressure.
Essentially, the leakage magnitude measures the tendency of leaked gas to mobilize further, such as from
the leakage aquifer to the ground surface.
2.2.2 PhysicalModel
We consider a faulted aquifer system with three confined aquifer compartments, A1, A2, and A3, offset
across a through-going fault, as shown in Figure 2.2. We consider four conjugate models with the fault
dipping at θ =60
◦ ,80
◦ ,100
◦ , and 120
◦ angles measured counterclockwise from the ground surface. A2
is located in the hanging wall blocks of the 60
◦ and 80
◦ dip models and footwall blocks of the 100
◦ and
120
◦ dip models. The top and right boundaries of the domain are under compression with a horizontal-to-
vertical total stress ratio equal tor
tec
=0.7, which results in a normal faulting stress regime implying that
the hanging wall block of the fault prefers to slip down-dip. r
tec
>1 will result in a reverse faulting stress
19
0 2000
2500
500
x (m)
z (m)
Fault

Injection
A3
A2
A1
Footwall
Hanging
wall

Footwall
Hanging
wall
No-ow
No-ow
No-ow
A3
A2
A1
(a) (b)
θ = 80 ◦
100 ◦
Figure 2.2: A physical model of gas injection and leakage in faulted reservoirs. (a) Gas is injected in
the hanging wall block of a normal fault, which slips and begins to leak fluids into aquifers A1 and A3
depending on the slip direction and gas buoyancy. (b) A conjugate model similar to the model in (a) except
that the fault dips in the opposite direction such that the injection occurs in the footwall block.
regime. Gas is injected at a point in the A2 aquifer at a constant rate. We denote the 80
◦ dip model as
the HWI (Hanging Wall Injection) model and the100
◦ dip model as the FWI (Footwall Injection) model.
The porosity is assumed to be uniform throughout the domain. The permeabilities of the aquifer and
non-aquifer caprock are assumed to be uniform, isotropic, and constant in time. The fault permeability
is assumed to be higher across the aquifer depth intervals than across caprock, and it is assumed to vary
with time-dependent slip. Before slipping, the fault acts as a hydraulic barrier to flow along and across
its surface. After slipping, the permeability along the fault is assumed to increase with the amount of
shear slip (Lee & Cho, 2002; Sibson, 1994) following an empirical model (Eq. 2.12) (Dempsey et al., 2013).
See Table 2.1 and Figure 2.3 for the input parameters to the model problem. The numerical values of the
parameters, apart from those of the fault permeability model, are representative of the model problem
and have been described in previous publications (Cappa & Rutqvist, 2011; Jagalur-Mohan et al., 2018;
Jha & Juanes, 2014). We do not consider parameter uncertainty in this study. We choose to model with
the possibility of large fault permeability changes over seven orders of magnitude, which allow us to
investigate the extremes of leakage changes and better demonstrate our proposed metric, the leakage
20
magnitude. We choose geomechanical properties in accordance with field data. The poroelastic properties,
Young’s modulus and Poisson’s ratio, are representative of a shallow (∼ 1500 m depth) sandstone aquifer,
e.g. the Morrow-B Sandstone in the Farnsworth field in Texas (Zhao & Jha, 2019). The d
c
value falls within
the range of critical slip weakening distance reported in multiple studies Ide and Takeo, 1997; Ma and
Elbanna, 2015; Mikumo et al., 2003; Olsen et al., 1997; Scholz, 1988. All the values are chosen to realize
induced slip and leakage within the period of simulation (approximately six months), which is aligned
with the goals of this study.
2.2.3 InitialandBoundaryConditions
The domain is initially fully saturated with water, pressures are hydrostatic, and initial total stresses
are in equilibrium with boundary stresses and gravitational load. The boundary conditions (shown in
Figure 2.2) are no-flow on all boundaries, tectonic compression on the right boundary, overburden stress
on the top boundary, and roller boundary conditions (fixed in the direction perpendicular to the boundary,
zero traction parallel to the boundary) on the left and bottom boundaries. Symmetry boundary conditions
are used to model injection into a fault relay ramp pattern reservoir with less computational cost. Note
that the left-hand side imposes a symmetry condition that mimics simultaneous injection by two wells
into the A2 aquifer. The arrows on the top and right boundaries of the left figure indicate overburden and
tectonic compressions, respectively. The stresses increase linearly with depth under a lithostatic gradient.
The hypocenter locations of injection-induced seismicity, shown with the star symbols, are located on the
aquifer boundary where injection-induced shear is aligned with the tectonic shear.
21
Viscosity
Density
Viscosity
Density
Figure 2.3: Pressure-dependent density and viscosity of gas and water phases.
2.3 InducedFaultReactivation
The locations of the hypocenters during fault reactivation due to injection depend on the dip angle.
Injection-induced overpressure in the aquifer A2 results in volumetric expansion of the aquifer and re-
duction in the effective compression inside the aquifer (Jagalur-Mohan et al., 2018; Jha & Juanes, 2014).
Decrease in the effective compression, − ∆ σ ′
n
<0, and expansion-induced increase in the shear traction
magnitude,∆ |τ |>0, culminates into slip nucleation when the Fault Stress Ratio, FSR=|τ/σ
′
n
|, also known
as the slip tendency, at a point increases to be equal to µ s
. The slip nucleation points, or hypocenters
of induced seismicity (denoted as stars in Figure 2.4), are locations along the A2 aquifer–fault boundary
where the tectonic shear and the injection-induced shear act in the same direction. This is at the bottom
of A2 (depth 1550 m) forθ< 90
◦ and at the top of A2 (depth 1450 m) forθ> 90
◦ (Figure 2.4, top row). We
note that the overpressures can reach up to 15 MPa, which may lead to fracturing, but as we are interested
in investigating and quantifying the limiting cases of leakage, we proceed accordingly because the high
overpressures allow us to better investigate the limits of leakage. We show more detail regarding the slip
nucleation points vs. the fault dip angle in Figure 2.5. The CFF’s time evolution forθ =60
◦ ,80
◦ shows that
the CFF initially increases in the fault at the bottom of A2 despite the bottom of the aquifer being initially
more stable (Figure 2.5(a) and (b)). After slip nucleates at the bottom, adjacent sections of the fault begin
22
Figure 2.4: Effect of the dip angle θ on overpressure, gas saturation and displacement after 149 days of gas
injection. Induced slip and subsequent leakage are shown.
to slip and the slip nucleation point stabilizes. A non-intuitive phenomenon also occurs at the top of A2:
the CFF decreases at the top of A2, indicating that the fault region at the top of A2 is actually stabilized in
the early stages of injection. The opposite pattern is observed forθ =100
◦ ,120
◦ : the CFF increases and slip
nucleates in the fault region near the top of A2 while the fault region at the bottom of A2 is initially stabi-
lized (Figure 2.5(c) and (d)). The same post-slip phenomenon is observed at these angles: regions adjacent
to the slip nucleation points begin to slip while the slip nucleation point stabilizes after slip. These results
pertain to a normal faulting regime: in case of reverse faulting regime withr
tec
>1, the hypocenters will
be located at the top of A2 forθ< 90
◦ and at the bottom of A2 forθ> 90
◦ because the tectonic shear vector
on the injection side of the fault points up-dip.
23
A2
Depth (m)
CFF (MPa)
0 day
5 day
35 day
Depth (m)
CFF (MPa)
(a) (b)
(c) (d)
0 day
10 day
36 day
0 day
7 day
36 day
0 day
15 day
36 day
θ = 60 ◦
80 ◦
100 ◦
120 ◦
CFF (MPa)
CFF (MPa)
Figure 2.5: CFF profiles at three time steps showing the effect of the dip angle θ on the natural selection
of the slip nucleation points. The horizontal lines denote the boundaries of A2, the injection aquifer. Slip
begins at the depths where the CFF reaches zero. The nucleation time, indicated by the timestamp of the
dash-dot curve, is different for different θ . After slip (solid black curve), CFF drops at the nucleation point
due to slip-induced stress relaxation. The x-axis scale is different for the four figures.
24
(a) (b)
(c) (d)
Hanging Wall
FootWall
Figure 2.6: Fault dip θ affects the onset times of injection-induced slip (a) and leakage (b). (a) A2 Bot
and A2 Top indicate the hypocenter location in the hanging wall injection and footwall injection cases,
respectively. (b) A1 and A3 indicate the leakage location for the fourθ cases. Leakage of water is denoted
by× symbols and leakage of gas by◦ symbols. (c) The◦ symbol indicates the final state (last timestep)
of the simulation. The lines show the evolution of the leakage magnitude and seismic moment over time.
M
0
andM
leak,g
have units of N-m per unit length in the third dimension. (d) The dashed lines connecting
the circles are visual guides only.
25
0 10 20 30
0 0.5
0 10 20 30
1000
1500
2000
0 0.5
1000
1500
2000
0 10 20 30
0 0.5
0 10 20 30
0 0.5
-10 0 10 -10 0 10 -10 0 10 -10 0 10
1000
1500
2000
Slip (cm)
Vertical Depth (m)
A1
A2
A3
time
time
Slip (cm) Slip (cm) Slip (cm)
Gas saturation Gas saturation Gas saturation Gas saturation
Overpressure (MPa) Overpressure (MPa) Overpressure (MPa) Overpressure (MPa)
Vertical Depth (m) Vertical Depth (m)
A1
A2
A3
A1
A2
A3
0.25 0.75
0.25 0.75 0.25 0.75 0.25 0.75
Figure 2.7: Depth profiles of the fault overpressure, fault gas saturation (showing gas leakage into aquifer),
and slip along the fault are shown at three successive time steps (darker curves indicate later time steps
while lighter curves indicate earlier time steps). At early times, the fault slips downdip in the hanging wall
injection cases (θ = 60
◦ ,80
◦ ) and updip in the footwall injection cases (θ = 100
◦ ,120
◦ ). Depths of the
three aquifers–A1, A2, and A3–are marked with black horizontal lines in each row. Up-dip slip is positive.
26
The slip onset time,t
slip
, depends on the following quantities: the initial horizontal-to-vertical tectonic
stress ratior
tec
, fluid and bulk densities, the Biot coefficient b, the dip angleθ , and the rate at which FSR
increases with time to reach µ s
. The effect of dip on t
slip
can be understood from the application of the
Coulomb theory to an elastic medium, where the optimum angles for failure are given bytan2θ =∓ 1/µ
s
,
assuming that the maximum principal stress is vertical. As a result, θ =60
◦ is closer to Coulomb failure
compared toθ =80
◦ , andθ =120
◦ is closer to Coulomb failure compared toθ =100
◦ . Also, two models with
faults dipping atθ and180
◦ − θ are equivalent in terms of their proximity to failure. However, the rate of
increase in FSR is different at different depths, because it is a ratio of stresses. So, the two models with
faults atθ and180
◦ − θ will have different t
slip
because their hypocenter depths are different. This is what
we observe in Figure 2.6(a), where angles closer to90
◦ have later slip times andt
slip
changes depending
on whether or not injection occurred in the hanging wall or footwall. The rate of increase of FSR during
fluid injection or production is a function of the pressurization rate ( dp/dt) and the traction increase rates
(dσ
′
n
/dt anddτ/dt ). Under the assumption of undrained deformation (Jaeger & Cook, 1979), which may be
applicable for a low permeability locked fault prior to slip-induced permeability enhancement, the traction
increase rates can be expressed in terms of the pressurization rate and the fluid injection rate. Since the
fault stresses increase with depth due to the lithostatic gradient, the rate of increase of FSR is larger for
θ< 90
◦ than forθ> 90
◦ . As a result, FSR reachesµ s
earlier in HWI than in FWI, as shown by comparingt
slip
in Figure 2.6(a) forθ =60
◦ vs. 120
◦ andθ =80
◦ vs. 100
◦ . The slip propagates away from the hypocenters
in the HWI and FWI models in both up-dip and down-dip directions with preferential slip direction being
down-dip in HWI and up-dip in FWI model, respectively (Figure 2.7, bottom row).
After slip nucleation at the hypocenter, the direction of slip at other points along the fault is determined
by the competition between: (1) the shear stress transferred from the hypocenter slip (Jagalur-Mohan
et al., 2018; Jha & Juanes, 2014), (2) the shear from aquifer expansion induced by injection, and (3) the
leakage-induced overpressures in the leakage aquifer. In the HWI model, the top of A2 slips up-dip because
27
expansion-induced up-dip shear is larger than the down-dip shear transferred from the hypocenter (Figure
2.7, bottom row). In the FWI model, the bottom of A2 slips down-dip because the expansion-induced down-
dip shear is larger than the up-dip shear transferred from the hypocenter. As a result, the cumulative slip
region expands outward from the injection aquifer A2 in both HWI and FWI models. The cumulative slip
region is larger in case of HWI in part because slip begins earlier (t
slip
is smaller) than in FWI (Figure
2.6(a)). However, this does not completely account for the difference in the size of the cumulative slip
region because the effect of leakage on slip must be investigated.
2.4 InducedLeakage
A delayed and smaller leakage amount leads to larger magnitudes of overpressure and induced slip in
HWI compared to FWI. Att
leak
, the fault permeability increases rapidly with slip (Figure 2.8) to establish
hydraulic communication between A2 and A3 (in case of HWI) and A2 and A1 (in case of FWI). t
leak
depends on the following quantities: t
slip
, the slip weakening model, the permeability-slip model, and the
hydraulic potential difference for multiphase flow through the fault. We call A3 the leakage aquifer in HWI
and A1 the leakage aquifer in FWI. The general characteristics of leakage in both models are as follows: a
rapid decline in the A2 pressure due to water exiting A2 (Figure 2.7, top row,80
◦ ,100
◦ ), increase in the gas
saturation in the fault zone (Figure 2.7, middle row), and increase in the gas saturation in the leakage aquifer
(Figures 2.4 and 2.7, middle rows). Forθ = 60
◦ and120
◦ , the overpressures do not decline substantially
(Figure 2.7, top row). For HWI and FWI, the aquifer gas saturation increases rapidly aftert
leak
is reached,
which is 125 days in the HWI model and 99 days in the FWI model (Figure 2.8). As a result of largert
leak
in
HWI, the overpressure reaches a higher value (Figure 2.8, third from top) that contributes to a larger slip
area and seismic magnitude in HWI. Again, we note that Figure 2.8 shows aquifer overpressures reaching
30 MPa. This may lead to fracturing. However, the objective is to study the limits of fault slip and leakage
using our methods, not fracturing. Consequently, we proceed with our analysis using these overpressures.
28
Figure 2.8: Dynamics of fault leakage in the hanging wall injection (HWI) and footwall injection (FWI)
models in terms of evolution in fault permeability, fault slip (up-dip is positive), and overpressure and gas
saturation in aquifers A1, A2 and A3. The green vertical lines denote the onset of leakage:t
leak
=125 days
(HWI) andt
leak
=99 days (FWI).
29
Leakage occurs earlier in FWI (Figure 2.6(b) and (d)) because both the gas buoyancy and preferential
slip vectors point in the same up-dip direction. In other words, the multiphase flow potential difference
overwhelms the difference in t
slip
between HWI and FWI for the dip angles investigated. Buoyancy also
plays a role in determining the leakage magnitude, which is smaller in HWI even though the pressure drop
during leakage and induced slip are larger in HWI. This leads to our first key observation: the leakage
magnitude can be smaller in the case with larger induced slip. This is shown explicitly in Figure 2.6(c),
which shows that the relationship between M
0
and M
leak,g
varies with the dip angle. Note that these
quantities could have no correlation depending on the circumstances. As the maximum M
0
increases
from approximately 2.0e10 N-m (per unit length in the third dimension) in θ = 80
◦ to 3.4e10 N-m in
θ =60
◦ , the maximumM
leak,g
increases from 3.0e9 N-m (per unit length in the third dimension) to 1.4e10
N-m. A similar observation can be made for FWI. However, if we consider HWI and FWI cases together,
there is no universal relationship between the quantities; the maximumM
leak,g
forθ =100
◦ is larger than
that forθ = 80
◦ even thoughM
0
forθ = 100
◦ is smaller than that forθ = 80
◦ . We also make another
observation regarding the relationship between the two quantities: forθ = 60
◦ ,120
◦ , the two quantities
evolve proportionately, and forθ =80
◦ ,100
◦ ,M
leak,g
evolves almost independently ofM
0
.
The second key observation results from the time evolution of gas saturation during leakage (Figure 2.8
top). The evolution profile is monotonic in case of FWI because both the buoyancy and slip velocity vectors
point upward. A similar observation has been made earlier in case of leakage through a vertical fault with
constant permeability (Pruess, 2005). However, the prior work did not account for fault slip or coupling
from flow to deformation. Our coupled multiphase flow-deformation-slip model improves on the prior
work and provides a mechanistic explanation for the evolution of the leakage event magnitude M
leak,g
,
which is defined in terms of s
g
.
The third observation is that the gas saturation can evolve non-monotonically in time due to acompe-
tition between the upward direction of the fluid buoyancy force and the down-dip direction of the induced
30
slip velocity. Figure 2.8 top shows the non-monotonic gas saturation in the HWI model’s leakage aquifer
(A3). This competition evolves in time with the evolving slip magnitude in the HWI model. Also, the mag-
nitude of the buoyancy force changes mildly with time due to the gas density increasing with pressure
(the fluid compressibility effect) as leakage continues into the HWI model’s leakage aquifer.
(a) (b)
(c) (d)
Figure 2.9: A comparison of the evolution of the leakage magnitude (M
leak
) and the leakage mass (m
leak
)
as a function of the dip angle. The leakage magnitude and leakage mass have units of N-m and kg, respec-
tively, per unit length in the third dimension.
We compare the leakage magnitude with the more traditional measure, the leakage mass, to show how
it can provide additional information about leakage. First, we define the leakage mass as follows:
m
leak,g
(t)=
Z
V
off
ϕρ g
(x,z,t)s
g
(x,z,t)dV, (2.19)
31
where m
leak,g
(t) is the leakage mass at time t and ρ g
(x,z,t) is the gas density at a particular x and z
location in the leakage aquifer at timet. Note that the differences between the leakage magnitude and the
leakage mass are due to the differences between density and overpressure in the reservoir. This difference
can be used to obtain additional information about leakage. Figure 2.9 shows the variation in leakage
magnitude and leakage mass over time for different dip angles. We choose axis values such that relative
changes in the leakage magnitude and leakage mass show as equally large on the page. This enables
easier comparisons. Figure 2.9(a) shows that while the leakage magnitude and leakage mass are generally
correlated, the leakage mass increases faster than the leakage magnitude. Since this difference must be due
to density vs. overpressure, it indicates that that flow initially occurred in the reservoir without generating
a proportional increase in overpressure. The leakage magnitude then rises steeply while the leakage mass
exhibits no change in slope, indicating that the fluid in the reservoir is pressurizing without a substantial
change in gas flow rate. Figure 2.9(b) shows that both the leakage magnitude and leakage mass flatline
soon after leakage begins. This shows that gas flow ceases into the reservoir and that the reservoir is
no longer pressurizing, suggesting that the reservoir has stabilized. Figure 2.9(c) shows that both the
leakage magnitude and leakage mass are increasing. This indicates that gas is flowing into the reservoir
and the reservoir is pressurizing. Figure 2.9(d) shows that leakage mass increases linearly but the leakage
magnitude initially increases sublinearly and then increases to a superlinear slope. This shows that the
reservoir initially pressurizes slowly as gas flows in and then pressurizes more quickly for each unit of
entering gas as the leakage mass increases. This means that while an operator using leakage mass would
see a linearly rising risk, using leakage magnitude would show that the risk rises faster than the leakage
mass indicates in the later stages of leakage due to the pressurization of fluids in the reservoir. These
comparisons show the value of the leakage magnitude as a complementary measure to leakage mass by
demonstrating the additional information that can be obtained about the pressurization of the reservoir.
32
65 day
θ = 60 ◦
θ = 80 ◦
θ = 100 ◦
θ = 120 ◦
U z U z U z U z (a) (b)
(c) (d)
Figure 2.10: The displacements at the upper boundary vary non-monotonically with time. Effect of the
fault dipθ on this evolution is evident.
Time (days)
0 20 40 60 80 100 120 140 160
0
0.01
0.02
0.03
0.04
0.05
0.06
(m)
60 ◦
120 ◦
80 ◦
U z θ = 100 ◦
Figure 2.11: The vertical displacements at the top boundary of the domain near the well’s horizontal loca-
tion are presented for each dip angle.
33
2.5 EffectofLeakageandFaultDiponGroundDeformation
The bottom row of Figure 2.4 and their time evolution in the supplementary videos (available online
as part of Meguerdijian and Jha, 2021) show new results with regard to the idea of monitoring leakage
using ground surface displacement data from InSAR (Y. Jung et al., 2013; Morris et al., 2011) and GPS.
First, leakage of both water and gas along low-angle faults (θ =60
◦ ,120
◦ ) leads to a non-uniform vertical
displacement profile across the fault, which may imply an easier leakage detection than in θ = 80
◦ ,100
◦ cases with relatively uniform displacement profiles. Second, it is non-trivial to relate the displacement
profile anomaly to the leakage magnitude in θ = 60
◦ ,120
◦ , because the displacement above the offset
aquifers is a superposition of displacements due to leakage of both water and gas phases.
The effect of leakage on top boundary displacements can be seen in Figure 2.10. Initially, as injection is
performed, the top boundary rises due to the increase in pressure (see supplementary videos dip80_uz and
dip100_uz). After some time, leakage from the injection aquifer occurs. The leakage of gas or water from
an aquifer reduces pressure within the aquifer when the leakage rate exceeds the injection rate, causing
subsidence. The continuation of injection eventually restores pressure in the aquifer despite gas/water
leakage, resulting in a rise at the top boundary again. Consequently, non-monotonic surface displacements
above injection aquifers may indicate leakage, especially in case of an injection setting that is similar to
the one in Figure 2.2, i.e. symmetric well injection into a largely homogeneous reservoir near a fault.
The effects of leakage on vertical displacements can be quite pronounced near the well. Figure 2.11
shows that in theθ = 80
◦ and100
◦ cases, the rate of increase in vertical displacements was nearly linear
until leakage began, after which the vertical displacements began to decrease. In the θ = 60
◦ and 120
◦ cases, leakage occurred much earlier (see the local maximum around 35 days) and the decline was less
pronounced. Higher-angle faults appear to have larger vertical displacements prior to leakage, but all
faults exhibit a change from a rising trend to a subsiding trend during injection when leakage begins.
34
This suggests that if a previously rising surface at the site of a gas injection well begins to subside during
injection, leakage may occur.
The maximum displacements at the upper boundary occur at the side boundaries (Figure 2.10). There
may be concerns that boundary effects in these simulations may cause pressure to increase faster than is
physical. However, considering the scenario of symmetric injection into an aquifer in a fault ramp relay
pattern, the boundary effects that appear from these simulations are in fact realistic. A physical explana-
tion for the maximum displacements occurring at the side boundaries would be that the presence of sealing
faults restricts the flow of fluid, and fast pressure diffusion leads to the reservoir acting like a gas storage
vessel. This would result in the aquifer experiencing a “balloon”-like deformation, with maximum defor-
mations occurring at its midpoint, equidistant from the sealing faults on either side. Notably, this manner
of deformation would be observed in the case of a homogeneous, isotropic aquifer without fractures or
permeability changes.
There are other possible factors that could show a similar displacement response in a field or exper-
imental setting. Material property heterogeneities could smooth out the surface displacement response
above reservoirs. Deformations of the surface can be due to mechanisms other than leakage (Rinaldi &
Rutqvist, 2013; Rucci et al., 2013; Vasco et al., 2018). Large changes in the fault permeability amplify the
deformation effects we observe here, which is on the order of a few centimeters. Smaller values are likely
to be observed in practice. Field applications should use the results from this study in combination with
other data, such as well data or geophysical measurements, to ascertain the causes of a given deformation
pattern.
35
2.6 Conclusion
Injection-induced fault slip can lead to an increase in the along-fault permeability, which can cause
leakage into and pressurization of the offset aquifers across the fault. Using our coupled multiphase flow–
geomechanics simulator, we simulate leakage scenarios with large changes in the fault permeability. We
determine the leakage characteristics of a fault as a function of fluid buoyancy, faulting environment, and
induced slip on the fault. Our conclusions are as follows:
1. The propensity of gas to leak can be quantified using the proposed numerical measure of gas leakage,
the leakage magnitude (Eq. (2.18)). Combined with leakage mass, the leakage magnitude can indicate
additional information about a reservoir’s pressurization.
2. We observe an asymmetry in the leakage profile along the fault between down-dip and up-dip prop-
agating slip events due to the buoyancy of the leaking fluid and the asymmetry in the induced slip
profile between hanging wall and footwall injection scenarios for symmetric injection into a fault
ramp relay pattern aquifer (Figure 2.8).
3. The leakage magnitude associated with a smaller induced slip event can be larger than the leakage
magnitude associated with a larger induced slip event. This suggests that we cannot equate the risk
of induced seismicity with the risk of leakage (Figure 2.6).
4. With respect to injection monitoring and leakage detection via geodetic means (InSAR, GPS), we
conclude that the magnitude of the leakage-induced ground deformation anomaly across the fault is
sensitivetothedipangleandcombinestheeffectofbothwaterandgasleakage (Figures 2.10 and 2.11).
We also note that other causes of ground deformation could be influential in other scenarios.
Finally, coupled modeling of multiphase flow and deformation processes while accounting for the fault
geometry and regional stress regime is essential for site selection and operation of underground fluid
36
storage projects. Future work includes incorporating the effect of fracturing and failure mechanics in the
aquifers with fault leakage and fault dip effects. For faults within layers of different lithologies (carbonates
vs. sandstone), faults with different mineral compositions (e.g., clay content), or faults at different stages
of their earthquake cycle, slip-induced alteration of the fault permeability may be different from the model
used here (Eq. (2.12)) and may require analysis of lab or field data (Chaves et al., 2020; Duan et al., 2020;
Okazaki et al., 2013).
37
Chapter3
ThermalandSolubilityEffectsonFaultLeakageduringGeologic
CarbonStorage
Abstract
Geologic carbon storage (GCS) is a promising method for reducing anthropogenic CO
2
emissions to the
atmosphere. To safely deploy GCS in the field, it is necessary to assess risks and the effect of uncertainty
on safe storage. The effect of uncertainty can be quantified using batches of simulations, but the high com-
putational costs of high-resolution simulations necessitate use of reduced-order models (ROMs). Previous
work involves ROMs for quantifying the risk of different potential leakage paths from storage reservoirs to
shallow formations. However, previous studies on development of fault-leakage ROMs have limited num-
bers of uncertain parameters and do not explicitly examine impacts of CO
2
solubility and thermal stresses
on fault reactivation, which can generate high permeability pathways and compromise CO
2
storage. In
this study, we analyze an ensemble of simulations considering CO
2
leakage from a storage reservoir to a
shallow aquifer through a fault while varying a number of uncertain parameters related to thermo-hydro-
mechanical properties and CO
2
injection. We show the effects of solubility on: free-phase CO
2
-leakage
rates, brine-leakage rates, and poroelastic fault destabilization. We find that CO
2
solubility is more im-
portant for estimating free-phase CO
2
-leakage rates compared to brine-leakage rates or poroelastic fault
38
destabilization. We also find that thermal stresses and overpressures have different spatial distributions
within the fault, indicating that the spatial variability of overpressures due to variation in flow parameters
does not necessarily make the spatial variability of thermal stresses negligible. We suggest the use of the
CO
2
phase-change path as a variable in future fault-leakage ROMs.
Keywords: geologic carbon storage, reduced-order model, fault leakage, CO
2
, fault reactivation
3.1 Introduction
Geologic carbon storage (GCS) is considered a promising method for reducing CO
2
emissions to the
atmosphere (IPCC, 2005). With the enactment of the recent §45Q statute (26 USC §45Q, 2019), many billions
of dollars of tax credits are set to incentivize the growth of carbon capture and storage significantly over
the next decades. However, concerns exist over the safety of CO
2
storage sites (Newell & Ilgen, 2019;
White & Foxall, 2016). Risk assessment is essential to both the technical and economic viability of GCS
projects. Probabilistic ensembles of simulations are a prominent way to quantify risk under uncertainty
because of their ability to integrate the various complexities of the subsurface; however, they can often be
too computationally expensive or numerically challenging to execute. Reduced-order models (ROM) are
one way to overcome this problem and so are favored for applications such as site-screening (Middleton
et al., 2020). ROMs reduce the computational burden of simulations by capturing only the physics essential
to a given risk assessment problem. ROMs have been applied to field cases (Onishi et al., 2019).
Comprehensive risk assessment includes identification of uncertain variables and impact of uncertainty
on risks. It involves the risk of fault reactivation and subsequent change in fault permeability, although
these effects are not straightforward to capture (Rinaldi et al., 2014; F. Zheng et al., 2021). A fault may reac-
tivate and cause substantial changes in permeability (Guglielmi et al., 2017; Tran & Jha, 2021), resulting in
a high permeability pathway to the surface, potentially causing loss of containment of stored CO
2
and as-
sociated geomechanical deformation patterns across the fault (Meguerdijian & Jha, 2021). CO
2
is partially
39
miscible in brine. The dissolution process depends on CO
2
pressure and temperature, both of which can
change with time as CO
2
migrates within the reservoir and potentially through a fault. It is believed that
monitoring dissolved CO
2
in brine can assist with monitoring CO
2
plume migration, solubility trapping,
and groundwater contamination risk associated with CO
2
storage (Keating et al., 2010). However, a limita-
tion of the high computational cost of storage simulations is that, for the purposes of producing ROMs, the
sample size of simulations to adequately cover the parameter space of even a few uncertain variables can be
quite large. When combined with coupled flow and geomechanics, the computational burden can become
substantial, such that even training ROMs for use in uncertainty quantification becomes challenging.
In previous research, ROMs have been generated for simulations of wellbore leakage (D. R. Harp et
al., 2016), aquifer impacts (Bianchi et al., 2016; Keating et al., 2016), fault leakage without geomechanics
(Lu et al., 2012), and preliminary work on fault leakage as either a plane (Lindner, 2020) or a continuum
with aseismic slip (White et al., 2017). These fault-leakage ROM studies include those involving a limited
number of uncertain variables (around five, though possibly after a sensitivity analysis of potential inputs)
or simplified physics, such as the assumption of parallel plate flow. While these studies provide important
information regarding the effects of fault slip or initial steps toward development of a fault-leakage ROM
that can be incorporated into risk assessment frameworks, such as NRAP-IAM-CS (Pawar et al., 2016),
they do not explore the effects of solubility or thermal stresses in detail.
In this paper, we identify insights for future ROM development by investigating the physical process
of fault leakage. We present an analysis of brine and CO
2
leakage rates through faults and their association
with fault destabilization induced by CO
2
injection, with particular focus on CO
2
solubility and thermal
stresses. We use a 3D model to simulate leakage from a storage reservoir through a fault to a shallow
aquifer, and we consider 19 uncertain variables associated with fluid flow through faults. We compute
poroelastic destabilization, temperature changes, and thermal stress changes within the fault across an
ensemble of realizations. We consider the effects of solubility on fault stability. We show that in spite
40
of substantial variations in input parameters, thermal stresses can still exert significant effects on fault
destabilization and thereby on changes in fault permeability and resulting leakage rates. For the parameter
space used in this study, we note that the spatial distribution of thermal stresses in the fault is different
from that of overpressures. The two effects occur in different portions in the fault in such a manner that the
spatial variability of one does not render the spatial variability of the other negligible. Moreover, we show
that CO
2
dissolution affects free-phase CO
2
leakage rates, where “free-phase” refers to liquid, vapor, and
supercritical phases combined but not CO
2
dissolved in brine. The effect of solubility on brine leakage rates
and poroelastic fault destabilization is minimal. We draw conclusions as to how to generate fault-leakage
ROMs which incorporate these findings.
3.2 Methods
3.2.1 DescriptionofthePhysicsandConstitutiveEquations
We used the Finite Element Heat and Mass Transfer (FEHM) code to perform fully-coupled 3D multi-
phase simulations with heat and mass transfer (Zyvoloski et al., 2015). The relevant equations have been
described previously (Zyvoloski et al., 1999). We used the Pylith-FlowSim simulator (Jha & Juanes, 2014;
Tiwari et al., 2021; Zhao & Jha, 2019; Zhao & Jha, 2021) to conduct geomechanical analyses using the out-
puts of FEHM and one-way flow-to-mechanics coupling. Given that this study is concerned with the risks
to fault stability due to injection in a deep aquifer and leakage from a deep aquifer to a shallow aquifer, it
is safer to assume that the mechanics-to-flow coupling is negligible, as the presence of mechanics-to-flow
coupling would tend to produce lower pressure estimates and consequently lower leakage and risk of fault
reactivation.
Mechanics-to-thermal coupling is considered insignificant because heat flow occurs by two main mech-
anisms in the reservoir: (1) conduction and (2) convection. Regarding conduction, the thermal conductivity
41
of a cell depends on the relative fractions of brine, rock, and supercritical CO
2
(scCO
2
). Mechanics would
change this by changing the porosity. However, the porosity changes are expected to be small, as the
change in porosity due to pressure is proportional to the inverse of the drained bulk modulus (Biot and
Willis, 1957; Coussy, 2004, Eqns. 4.21c, 4.35, 4.44; as given by White et al., 2019, Eqn. 4). The drained bulk
modulus of a reservoir can be on the order of GPa, leading to small porosity changes in the response to
pressure changes on the order of MPa. Thermal contraction would also be small, considering that the vol-
umetric thermal expansion coefficient of many common minerals is on the order of 10
− 5
/
◦ C (Robertson,
1988, Table 9). This would limit the amount of convection which would occur as a result of thermome-
chanical deformation. Consequently, in this study, one-way coupling is safer for the estimation of fault
reactivation risk than two-way coupling.
The numerical formulation of (Jha and Juanes, 2014) is used to calculate poroelastic effects while ther-
mal stresses in the fault are calculated separately. The change in the Coulomb Failure Function (CFF) (also
known as the Coulomb failure stress) (Reasenberg & Simpson, 1992; Scholz, 2018) is used to estimate fault
destabilization due to poroelastic effects:
∆ CFF=∆ |τ |+µ f
∆ σ ′
n
(3.1)
whereτ is the shear stress on the fault, µ f
is the fault’s static coefficient of friction, and σ ′
n
is the effec-
tive normal stress (assumed positive in tension) on the fault. We consider only changes in the CFF in our
work to avoid specifying an initial stress state, which can have substantial uncertainties due to its methods
of measurement (Burghardt and Appriou, 2021) and is unnecessary for establishing this study’s conclu-
sions concerning fault destabilization. The thermal stresses are volumetric in nature and are calculated as
(Coussy, 2004, Eqn. 4.24):
σ T
=− 3α s
K
dr
(T − T
0
)1 (3.2)
42
where1 is the rank-2 identity tensor,K
dr
is the drained bulk modulus,T is the current temperature,T
0
is the initial temperature, andα s
is the linear thermal expansion coefficient of the solid phase. Note that
cooling of a fault can lead to destabilization depending on the reservoir conditions (Jacquey et al., 2015).
3.2.2 TheModel
Our conceptual model is shown in Figure 3.1. The model consists of a deep aquifer, a caprock, a shallow
aquifer, overburden, and a two-part fault. The fault consists of a fault core and damage zone and follows the
description given by (Johri, 2012). The fixed parameters of the reservoir model are listed in Table 3.1. The
varying parameters of the reservoir model are shown in Table 3.2 along with their respective ranges. The
values of the varying parameters are sampled uniformly (uniformly in the logarithmic domain for perme-
abilities) or with equal probability for the discrete values using Latin Hypercube Sampling (LHS) (Helton
& Davis, 2002). Note that a reduced-order model is not generated in this study, but rather an ensemble of
realizations across the range of these varying parameters is analyzed. The FEHM flow simulation model is
constructed using a mesh (termed the “FEHM mesh”) that has evenly spaced cells in all directions, except
for a geometrically spaced region in the x-direction near x=0 m to resolve the fault. The cell centroids of
this mesh are then used to make another mesh (termed the “Pylith mesh”). The Pylith mesh is used with
Pylith-FlowSim, and pressure values from the FEHM mesh are assigned to the centroids of the Pylith mesh.
This allows us to couple FEHM to Pylith to perform geomechanical simulations using pressures directly
from FEHM. Constructing the Pylith mesh in this manner avoids any loss in accuracy due to interpolation
of pressure, which will be required in case of two arbitrary meshes for flow and geomechanics.
The reservoir is initialized to a steady-state geothermal and hydrostatic gradient. The boundary con-
ditions for fluid flow and heat transfer are no-flow on all boundaries, except for the shallow aquifer, which
has constant pressure outflow-only boundaries for both brine and CO
2
. This increases the pressure gradi-
ent from the injection well to the shallow aquifer and so promotes leakage, enabling a conservative leakage
43
Overburden
Shallow Aquifer
Caprock
Deep Aquifer
Damage Zone Fault Core
Overburden
Shallow Aquifer
Caprock
Deep Aquifer
Slice Across Y-Axis
Side View
x
z
y
z
Map View Slice Across Shallow/Deep Aquifer
y
x
CO
2
Injection Well
Free surface
4000 m
No flow
No flow
No flow
4000 m
200 m
225 m
225 m
Deep/Shallow Aquifer
4000 m
4000 m
Damage Zone Fault Core
No flow
(Shallow Aquifer)
No flow
(Deep Aquifer)
CO
2
Injection Well
(Deep Aquifer)
200 m
225 m
225 m
3750 m
3750 m
p = constant
p = constant
(a)
(b) (c)
Figure 3.1: A schematic of the fault-leakage model problem.
44
Parameters Value
Deep aquifer thickness (m) 225
Caprock thickness (m) 3750
Shallow aquifer thickness (m) 200
Overburden thickness (m) 225
Density of the rock (solid phase) (kg/m
3
) 2650
Specific heat capacity of the rock (solid phase) (J/(kg · K)) 920
Surface temperature (at the top of the domain)(
◦ C) 10
Geothermal gradient (
o
C/km) 8.4
Temperature at the bottom of the domain (
◦ C) 47
Pressure at the top of the domain (MPa) 0.1
Porosity of the caprock (dimensionless) 0.1
Porosity of the overburden (dimensionless) 0.1
Permeability of the caprock (m
2
) 10
− 18
Permeability of the overburden (m
2
) 10
− 18
CO
2
injection temperature (
◦ C) 32
Biot coefficient of the domain (dimensionless) 0.8
Biot coefficient of the fault ( b, dimensionless) 0.8
Coefficient of static friction for the fault ( µ f
, dimensionless) 0.6
Height of the domain (m) 4400
X-width of the domain (m) 4000
Y-width of the domain (m) 4000
Well top z-coordinate (m) 225
Well bottom z-coordinate (m) 0
Total injection time (years) 10
Total simulation time (years) 30
Relative permeability model linear
Fault x-coordinate (m) 0
Fault y-cooordinates (m) [-1700, 1700]
Fault z-coordinates (m) [75, 4175]
Table 3.1: Fixed parameters of the reservoir model.
45
Parameters Range/Set
Fault core permeability (log
10
(m
2
)) [-18,-12]
Fault core porosity (dimensionless) [0.001,0.1]
Damage zone permeability (log
10
(m
2
)) [-15,-12]
Damage zone porosity (dimensionless) [0.001,0.1]
Shallow aquifer permeability (log
10
(m
2
)) [-14,-12]
Shallow aquifer porosity (dimensionless) [0.05,0.5]
Deep aquifer permeability (log
10
(m
2
)) [-14,-12]
Deep aquifer porosity (dimensionless) [0.05, 0.35]
Fault core thickness (m) [0.0, 0.5]
Damage/fault zone thickness (includes fault core thickness) (m) {33.33,70.85}
Fault permeable length (m) {800,1600,2400}
Rock conductivity (W/(m· K)) [1.0, 6.0]
Well x-coordinate (m)* {154.26,248.10,388.86}
Well y-coordinate (m) {-200,0,200}
Well rate (kg/s) [7.5, 25]
Solubility of CO
2
{soluble, insoluble}
Volumetric thermal expansion coefficient (
◦ C
− 1
) [10
− 5
,5× 10
− 5
]
Poisson’s ratio (dimensionless) [0.2,0.35]
Young’s modulus (GPa) [5,20]
Table 3.2: Varying parameters of the reservoir model. Square brackets denote a continuous, uniform dis-
tribution while brackets denote a discrete distribution. *The x-axis is at the center of the fault.
estimate. Note that pressure does increase within the deep aquifer during injection under these boundary
conditions. Injection occurs through the well at depths of 4175-4400 m. The CO
2
injection rate is constant
during the injection period and zero after the injection period ends. The depths are chosen such that the
CO
2
passes from a supercritical state to a gas state as it migrates from the deep aquifer to shallow aquifer
through the fault. The flow and transport solvers are coupled. Density changes with pressure, tempera-
ture, and dissolved CO
2
concentration. CO
2
solubility changes as a function of pressure, temperature, and
salinity. As capillary pressure would tend to inhibit CO
2
leakage, we neglect capillary pressures to obtain
conservative estimates of leakage and more easily draw conclusions from our simplified model.
The fault consists of a vertical fault core and damage zone whose properties are considered equivalent
properties (e.g. an upscaled permeability) and whose parameter variations are meant to reflect the variety
of conditions which can occur in the subsurface. Faults can exhibit complex, time-dependent behavior
46
with slip, such as slip-dependent permeability followed by leakage and fracture self-healing (Chester et al.,
1993), and numerous mechanisms can be involved in a fault’s history (Sibson, 1977). Fault core composition
varies with distance from the slip surface (Scholz, 2018), with asymmetry potentially across the fault (e.g.
pulverization on one side of the fault when rock stiffness values differ on each side of the fault (Dor et al.,
2006; Mitchell et al., 2011)). Fault zone conceptual models include differences in the architecture of the
fault (e.g. a single fracture vs. a deformation zone, fractured rock around the core vs. solid rock around
the core) (Caine & Forster, 1999). In this work, we opt for a fault model which focuses on the leakage
amounts and selected stresses, incorporating the uncertainties in the fault architecture through a two-part
fault model (fault core and damage zone), the use of equivalent permeabilities and porosities for the fault,
and a range of fault property values to capture the fault’s properties’ uncertainties. This is a simplified
model and is not applicable in all cases. We choose to increase the parameter uncertainties instead of the
model complexity to obtain clear causative relationships.
The permeabilities and porosities of the fault core and damage zone are obtained by LHS. These values
are then upscaled, with the resulting upscaled permeabilities and porosities assigned to the mesh vertices
in the x-direction nearest to the fault. The fault is modeled as a zone consisting of a plane of cells one-cell
thick on each side of the fault’s center plane when the damage zone is 33.33 m thick and a region two-
cells thick on each side of the fault when the damage zone is 70.85 m thick. When the damage zone is
70.85 m thick, the properties of the second sheet of cells furthest from the fault’s center plane are those
of the damage zone. In both cases, the first sheet of cells (those adjacent to the fault’s center plane, here
called the “inner fault zone cells”) have upscaled properties incorporating both damage zone and fault core
properties. These upscaled properties are calculated using the following equations:
k
∗ fc,x
=
x
ib
− x
c
L
fc
k
fc
+
x
ib
− L
fc
k
dz
(3.3)
47
k
∗ fc,y/z
=
k
fc
L
fc
+k
dz
(x
ib
− L
fc
)
x
ib
− x
c
(3.4)
ϕ ∗ fc
=
ϕ fc
L
fc
+ϕ dz
(x
ib
− L
fc
)
x
ib
− x
c
(3.5)
where:
• k
∗ fc,x
is the upscaled permeability in the x-direction of the inner fault zone cells,
• k
∗ fc,y/z
is the upscaled permeability in the y- and z-direction of the inner fault zone cells,
• ϕ ∗ fc
is the upscaled porosity of the inner fault zone cells,
• x
ib
is the thickness of the inner fault zone cells in the x-direction (equal to 33.33 m),
• x
c
is the x-coordinate of the center of the fault (x
c
=0 m in this case),
• L
fc
is the fault core thickness,
• k
fc
is the fault core permeability,
• k
dz
is the damage zone permeability,
• ϕ fc
is the fault core porosity, and
• ϕ dz
is the damage zone porosity.
A low geothermal gradient is chosen (the value in Table 3.1 is roughly that of an old continental craton
(Lowell et al., 2014)); this allows for a more conservative estimate of CO
2
solubility effects on ∆ CFF and
allows phase change to occur within the fault. Phase change produces both temperature changes (which
induce thermal stresses) and substantial drops in density in the fault.
3.3 ResultsandAnalysis
3.3.1 BaseCase
We performed a base case simulation to demonstrate the physics of the model problem. Table 3.3
shows the parameters used in the base case. The case does not take into account CO
2
solubility in brine,
48
but the effect of solubility is shown through comparison with results from the model with CO
2
solubility.
The damage zone permeability is about 7 mD while the deep aquifer permeability is about 389 mD. The
shallow aquifer permeability is about 405 mD. A corner of the simulation domain is plotted in the following
figures to show the properties from the simulation in both the fault and aquifers. Figure 3.2 shows this
cutout of a corner of the domain withx = 0 m corresponding to the fault. The injection rate is roughly
706,000 metric tons of CO
2
per year. Figure 3.3(a) shows that water movement occurs around the well and
within the fault both vertically up the fault and laterally across the fault, with greater water movement
occurring in the fault areas closest to the well and vertically above. Such 3D patterns of leakage have been
missing in some of the earlier 2D modeling studies of leakage. Water movement also occurs in the fault
section within the shallow aquifer. Figure 3.3(b) shows that water movement is substantially reduced both
laterally in the fault and within the shallow aquifer when solubility effects are included. This indicates
changes in the pattern of CO
2
spread, which we analyze further below.
While the CO
2
plume travels a smaller distance laterally in the fault when solubility effects are present,
the fault is nevertheless pressurized across its permeable length. A substantial amount of CO
2
dissolution
occurs within the fault. Figure 3.4(a) and (b) show that the supercritical CO
2
(scCO
2
) is present near the
injection well and within the fault, but the lateral spread of scCO
2
in the fault is reduced when solubility
is included. The same is observed with the length of the plume across the top of the fault, in the shallow
aquifer zone, in Figure 3.4(c) and (d). Figure 3.4(e) shows that when solubility is present, dissolution of
CO
2
limits the lateral spread of the free-phase CO
2
plume in the fault. This is important because it shows
how CO
2
may leak in both the free-phase (undissolved) and the aqueous phase (dissolved), and why CO
2
monitoring efforts should consider CO
2
solubility in cases where dissolved CO
2
may act as a signature of
CO
2
presence or risks such as elevated metal concentrations. The pressure has a different distribution than
the CO
2
saturation: Figure 3.5 shows that, while initially pressurized according to the hydrostatic gradient
in (a), the fault is laterally pressurized nearly equally along the fault’s permeable length (y-direction), as
49
shown in (b). This can be attributed to the relatively high deep-aquifer permeability in comparison to the
damage zone permeability.
The thermal stresses show a substantially different spatial distribution than the poroelastic ∆ CFF but
with a similar magnitude, suggesting that thermal effects are important for the evaluation of fault stability
in this case. Figure 3.6(a) and (b) show that the thermal gradient stays largely the same in the fault except
at the top, below the overburden. Figure 3.6(c) shows that the fault itself warms slightly, but the top of
the fault cools much more. This is due to CO
2
vaporizing from a liquid phase to a gas phase, absorbing
its latent heat of vaporization from it surroundings and causing cooling (see Figure 3.7 for the phase tran-
sition depths), and Joule-Thomson cooling from gaseous CO
2
expanding. In terms of thermal stresses,
Figure 3.6(d) shows that the thermal stresses within the fault are slightly stabilizing (compressive), but
those at the top of the fault are significantly tensile and could potentially reactivate the fault depending
on the proximity of the initial stress state to the Mohr-Coulomb failure envelope. In contrast, Figure 3.8(a)
shows that the poroelastic∆ CFF is greatest slightly above the bottom of the fault. This corresponds to the
overpressure (shown in Figure 3.8(b)) in the fault and above the deep aquifer. A difference appears between
the∆ CFF and pressure fields at the bottom of the fault, in the deep aquifer region and immediately above
it. As the pore pressures on both sides of the fault are more equal in the deep aquifer than in the fault
(where permeable damage zone is surrounded by impermeable caprock), the reduction in effective normal
stress is greater in the region within the caprock than in the deep aquifer. Therefore, while poroelastic
effects promote destabilization deep in the domain but above the deep aquifer in this case, thermal stresses
promote destabilization at shallow depths.
50
Parameters Value
Fault core permeability (log
10
(m
2
)) -13.8114
Fault core porosity (dimensionless) 0.0134959
Damage zone permeability (log
10
(m
2
)) -14.1377
Damage zone porosity (dimensionless) 0.00549412
Shallow aquifer permeability (log
10
(m
2
)) -12.3988
Shallow aquifer porosity (dimensionless) 0.444103
Deep aquifer permeability (log
10
(m
2
)) -12.4161
Deep aquifer porosity (dimensionless) 0.315762
Fault core thickness (m) 0.483355
Damage/fault zone thickness (includes fault core thickness) (m) 33.333
Fault permeable length (m) 2400
Rock conductivity (W/m*K) 5.99166
Well x-coordinate (m) 154.265
Well y-coordinate (m) 0
Well rate (kg/s) 22.3783
Solubility of CO
2
insoluble*
Volumetric thermal expansion coefficient (
o
C
− 1
) 2.15876*10
− 5
Poisson’s ratio (dimensionless) 0.20312
Young’s modulus (GPa) 7.91083
Table 3.3: The varying parameters for the base case. *The base case results are for the insoluble case unless
otherwise noted.
3.3.2 AutocorrelationofInputs
We sampled the uncertain parameters using LHS. A portion of numerical simulations did not complete
due to numerical instabilities. We used 217 simulations and verified that the sampling is unbiased to a first-
order approximation. Cumulative distribution functions of the inputs for the completed runs are shown in
Figure S.1 in the supplementary information. Figure S.1 shows that the input distributions are sufficiently
uniform.
It is apparent that the inputs do not have any substantial linear correlations between themselves. Fig-
ures S.2 and S.3 show that all the correlations have correlation coefficient values less than 0.3. Because of
this, first-order linear effects and interactions estimated from the completed runs can be taken as represen-
tative of the simulation physics. While linear relationships can be detected using Pearson coefficients, they
do not necessarily show monotonic nonlinear effects. We use Spearman coefficients to estimate first-order
51
Corner cutout
Shallow Aquifer
Permeable
Fault Location
Deep Aquifer
Caprock
Overburden
Figure 3.2: The location of the corner cutout and permeable portion of the fault in the base case figures.
The corner cutout and permeable fault location are shown in solid colors. The entire simulation domain
is shown using a wireframe drawing.
(a) Without solubility (b) With solubility
Shallow
Aquifer
Fault
Deep
Aquifer
Well
Figure 3.3: Corner cutout views of the water saturation field along the fault in the base case at the end of
injection, where (a) does not include solubility effects and (b) includes solubility effects.
52
Without solubility
(a)
With solubility
(b)
(c) (d)
(e)
Shallow
Aquifer
Fault
Deep
Aquifer
Well
Figure 3.4: Corner cutout of CO
2
saturations along the fault in the base case at the end of injection, where
(a) shows the scCO
2
saturation without solubility effects, (b) shows the scCO
2
saturation with solubility
effects included, (c) shows the gaseous CO
2
saturation without solubility effects, (d) shows the gaseous CO
2
saturation with solubility effects included, and (e) shows the dissolved CO
2
mass fraction with solubility
effects included.
53
(b) Pressure field at the end of injection. (a) Initial pressure field.
Shallow
Aquifer
Fault
Deep
Aquifer
Figure 3.5: Corner cutout of pressure evolution in the fault in the base case, where (a) is the initial pressure
field and (b) is the pressure field at the end of injection. The permeable length of the fault is evident in (b).
(a) Initial temperature field. (b) T emperature field at the end of simulation.
(c) T emperature difference
(end of simulation - initial).
(d) Change in thermal stress in the fault
(end of simulation - initial).
Shallow
Aquifer
Fault
Deep
Aquifer
Deep
Aquifer
Shallow
Aquifer
Figure 3.6: Corner cutout of temperature evolution in the fault in the base case, where (a) is the initial
geothermal gradient, (b) is the temperature field at the end of the simulation, (c) is difference in temperature
in the fault between the end of the simulation and the initial state, and (d) is the estimated thermal stress
resulting from the temperature difference. Note that normal stress is positive in tension.
54
Depth (m)
Gaseous CO
2
Liquid CO
2
scCO
2
Shallow
Aquifer
Deep
Aquifer
Figure 3.7: The CO
2
phase state as a function of depth in the fault under initial conditions.
Shallow
Aquifer
Deep
Aquifer
(a) (b)
Figure 3.8: (a)∆ CFF and (b) pressure change from the initial pressure in the fault in the base case at the
end of injection. Note that positive values of∆ CFF indicate fault destabilization.
55
nonlinear correlations between the input variables to ensure that significant nonlinear relationships do not
exist between the variables. Figures S.4 and S.5 show that the nonlinear monotonic relationships between
the variables are not significant, regardless of CO
2
solubility. We, therefore, proceed with our analysis of
the outputs and their first-order relationships to the inputs, noting that the inputs are uncorrelated to a
first-order and are generally uniformly distributed regardless of CO
2
solubility and thus can be used to
draw conclusions about first-order effects on CO
2
and brine leakage rates and CFF.
3.3.3 SolubilityEffectsonBrineandCO
2
LeakageRates
Brine and CO
2
leakage rates and their variations show that CO
2
dissolution does not significantly
affect brine leakage rates but does strongly affect CO
2
leakage rates, suggesting that solubility trapping,
while potentially taking hundreds of years to reach its full potential (IPCC, 2005, Fig. 5.9), can change
the rate of horizontal movement of free-phase CO
2
substantially in the first few years and decades of a
CO
2
injection operation. As Figure 3.9 shows, the maximum brine leakage rates are largely clustered at
low leakage rates for the range of parameters used in our study, with a moderate number reaching high
leakage rates. The shape of the distribution is the same regardless of whether or not solubility effects
are included. The distribution of maximum CO
2
leakage rates is starkly different with and without CO
2
solubility effects: when solubility effects are included, a substantial proportion of simulations yield no CO
2
leakage, as indicated by a jump in the cumulative distribution function (CDF) near zero. The CDF of CO
2
leakage rates reaches 1 at lower leakage rates when solubility is included. Using the maximum leakage
rates as proxies for overall leakage, Figure 3.9 suggests that the inclusion of solubility effects increases
the share of brine leakage relative to free-phase CO
2
leakage. This can be due to the free-phase CO
2
’s
horizontal mobility being reduced in the reservoir due to CO
2
dissolution in the brine. After injection,
the CO
2
rises toward the caprock and spreads thinly as a plume due to buoyancy. Since the edges of the
plume consist of a thin layer of CO
2
in a large contact area with brine, dissolution may occur faster than
56
Reduction in free-phase CO
2
leakage-rate extremes due to solubility
0 1.0 2.0 3.0 4.0 5.0 6.0
Maximum leakage rate (t/yr)
0.0
0.2
0.4
0.6
0.8
1.0
CDF
Brine, No Solubility
Brine, Solubility
CO
2
, No Solubility
CO
2
, Solubility
x1e5
Jump in zero-leakage cases due to solubility
Figure 3.9: Cumulative distribution functions of the brine and CO
2
maximum end-year leakage rates (the
maximum end-year leakage rate out of the 30 years of each simulation). The lower bound for denoting a
zero-leakage case for CO
2
is 10
− 4
t/yr.
during vertical migration of CO
2
toward the caprock. This dissolution reduces the relative permeability of
the edge of the CO
2
plume, resulting in substantially inhibited plume movement toward the fault. These
results, along with results showing the impact of vertical grid resolution on CO
2
dissolution (Pickup et al.,
2010), suggest that solubility must be carefully considered when making fault-leakage ROMs with respect
to free-phase CO
2
leakage rates.
3.3.4 EffectsonFaultReactivation
The poroelastic changes in the fault are sufficient to cause fault reactivation, and the driving factor is
the overpressure in the fault. Figures 3.10, 3.14, and 3.15 show results from the combined ensemble of 217
simulations, including those with and without solubility effects. Note that these figures show the entire
fault plane (the permeable fault location in Figure 3.2 and adjacent aquifer and caprock). Figure 3.10(a)
shows that the overpressures in the lower portion of the fault (below z = 1000 m) can reach above 10
MPa relative to initial conditions. Figure 3.11 shows that the mode of the max(∆ CFF) across all runs is
approximately 4 MPa. As ∆ CFF = ∆ |τ | + µ f
∆ σ ′
n
(Eqn. 3.1) and σ ′
n
= σ n
+ bp, ∆ CFF varies as µ f
bp,
57
where p is the pore pressure and σ n
is the total normal stress. In this case, µ f
b = 0.6× 0.8 = 0.48,
such that a 10 MPa overpressure causes a ∆ CFF change of 4.8 MPa, which is close to the average of the
values in Figure 3.11, with the difference being due to changes induced in the total traction and the shear
traction. The overpressure values in Figure 3.10(a) in the lower portion of the fault also correspond to the
maximum∆ CFF locations in Figure 3.12. Note that Figure 3.11 shows that the presence of CO
2
solubility
does not exert a clear impact on the distribution of max(∆ CFF), and Figure 3.12 shows that solubility does
not affect the areal distribution of max( ∆ CFF) in the fault. Figure 3.13 shows that virtually the entire fault
permeable region experiences some degree of destabilization, which Figure 3.10(a) strongly suggests is
due to pressurization. Pore pressure, regardless of CO
2
solubility, is the driving factor behind the fault
destabilization.
The greatest uncertainties in overpressure are deep within the fault and within the uncertainty in the
spatial extent of the permeable fault region. The standard deviation of overpressure (Figure 3.10(b)) shows
that the greatest uncertainties in overpressure are at around z = 425-675 m and y =± 700-900 m, which
reflect variability in the fault permeable length. These standard deviations reach about 6 MPa across the
ensemble of runs (both soluble and insoluble). The standard deviations are also larger immediately above
the aquifer (which terminates at a depth of 225 m). Both the average overpressures and standard deviations
of overpressure are substantially less significant at the top of the fault. The overpressures in the fault tend
to be higher the closer the y-coordinate is to those of the wells (0 m, ± 200 m), although the standard
deviations of overpressure do not change significantly in the y-direction except at fault permeable region
values. This shows that uncertainty in the fault permeable length introduces significant uncertainty in the
spatial distributions of overpressures, especially deep within the fault.
The thermal stresses caused by temperature changes in the fault exhibit significant destabilizing tensile
stresses at the top of the fault, distinctly different from the pressure distribution. Figure 3.10(a) shows the
largest overpressures below 1000 m depth, and Figure 3.10(b) shows the largest uncertainties also below
58
this depth. In contrast, the average temperatures in Figure 3.14(a) and (c) show that there is mild heating
in the fault but substantial cooling at the top of the fault, with the amount of cooling increasing between
the time of the maximum∆ CFF (Figure 3.14(a)) and the end of the simulation (Figure 3.14(c)). The stan-
dard deviations of these temperature changes are largest at the top of the fault and increase with time
(Figure 3.14(b) and (d)). This corresponds to the phase transition from liquid CO
2
to gaseous CO
2
, during
which the latent heat of vaporization is absorbed from the environment and Joule-Thomson expansion oc-
curs, causing cooling (cf. Figure 3.7). The thermal stresses corresponding to these temperature changes are
mildly compressive in the fault and tensile (destabilizing) at the top (Figure 3.15(a) and (c)). The standard
deviations of these thermal stresses are largest at the top and increase with time (Figure 3.15(b) and (d)).
The thermal stresses can potentially reactivate the fault at a shallow depth, under the overburden, which
is different from pressure-induced stresses, which may result in fault reactivation slightly above the deep
aquifer. Thermal stresses seem benign or small outside of the region where CO
2
phase change occurs.
These observations are confirmed by a comparison of the magnitudes of the thermal stresses versus the
pore pressure increase in the fault, shown in Figure 3.16. At the top of the fault (shown in Figure 3.16(a)),
the thermal stresses and fault’s pore pressure increase are comparable in magnitude, albeit different in
timing (year of max(∆ CFF) versus the end of the simulation). At the bottom of the fault (shown in Fig-
ure 3.16(b)), the distributions of thermal stresses and overpressures do not overlap, and the thermal stresses
are negligible compared to the overpressures. These results suggest that fault-leakage ROMs should cap-
ture the phase-change path in the fault, possibly as a variable, to incorporate the effects of thermal stresses
on fault reactivation and possible leakage pathway creation. The CO
2
phase-change path is the sequential
order of CO
2
phase states as CO
2
leaks up the fault based on the initial conditions. For example, Figure 3.7
shows that the CO
2
phase-change path in the base case is: scCO
2
, liquid CO
2
, and then gaseous CO
2
. This
path can, for example, be represented in a ROM by an indicator variable e.g., the path of scCO
2
, liquid
59
(b)
1500 1000 500 0 500 1000 1500
Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Z-Coordinate (m)
Overpressure in Fault
2
4
6
8
10
12
Overpressure (MPa)
1500 1000 500 0 500 1000 1500
Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Z-Coordinate (m)
Standard Deviation of Overpressure in Fault
1
2
3
4
5
6
Std. dev. (MPa)
(a)
Figure 3.10: All changes shown are deviations from initial conditions. (a) Ensemble average of the max-
imum pressure adjacent to the fault. Pressures taken at the time of the max(∆ CFF) for each run in the
ensemble. (b) Standard deviation of the ensemble pressures corresponding to (a).
0 2 4 6 8 10 12 14
Maximum CFF (MPa)
0.00
0.05
0.10
0.15
0.20
0.25
Probability Density Function
Soluble CO 2
Insoluble CO 2
Insoluble CO 2
Soluble CO 2
Figure 3.11: A histogram and kernel-density plot of the maximum∆ CFF values encountered in each run,
both for runs with CO
2
solubility and those without CO
2
solubility. Note that the soluble CO
2
histogram
overlays the insoluble CO
2
histogram and a third color appears where the two overlap.
CO
2
, and then gaseous CO
2
can be represented by a 1 if it is reflective of the phase transitions in the fault
of a given realization or 0 if it is not.
60
1500 1000 500 0 500 1000 1500
Fault Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Fault Z-Coordinate (m)
Soluble CO
2
Insoluble CO
2
Shallow aquifer
Deep aquifer
Caprock
Locations of max( ΔCFF) values
Figure 3.12: The locations of the points in the fault with the maximum∆ CFF change for each run. Note
that while the runs including solubility and not including solubility were plotted separately, the points
overlap.
0.992 0.994 0.996 0.998 1.000
Maximum permeable fault area fraction with CFF > 0.1 bar
0
200
400
600
800
1000
1200
1400
Probability Density Function
Soluble CO
2
Insoluble CO
2
Insoluble CO
2
Soluble CO
2
Figure 3.13: The fraction of the permeable area in the fault with∆ CFF> 0.1 bar.
61
(d)
(a) (b)
(c)
1500 1000 500 0 500 1000 1500
Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Z-Coordinate (m)
Average Increase in Temperature in Fault (End of Simulation)
0.3
0.2
0.1
0.0
0.1
Change in Temperature (
o
C)
1500 1000 500 0 500 1000 1500
Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Z-Coordinate (m)
Average Increase in Temperature in Fault (Year of Max( CFF))
0.00
0.02
0.04
0.06
0.08
0.10
0.12
Change in Temperature (
o
C)
1500 1000 500 0 500 1000 1500
Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Z-Coordinate (m)
Standard Deviation of Increase in Temperature in Fault (End of Simulation)
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
Std. dev. (
o
C)
1500 1000 500 0 500 1000 1500
Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Z-Coordinate (m)
Standard Deviation of Increase in Temperature in Fault (Year of Max( CFF))
0.1
0.2
0.3
0.4
0.5
Std. dev. (
o
C)
Figure 3.14: All changes shown are deviations from initial conditions. (a) Ensemble average of theaverage
temperature adjacent to the fault. Temperatures taken at the time of the max(∆ CFF). (b) Standard deviation
corresponding to (a). (c) Ensemble average temperature increase in the fault at the end of the simulation.
(d) Standard deviation of the ensemble temperatures in (c).
62
(b)
1500 1000 500 0 500 1000 1500
Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Z-Coordinate (m)
Average Thermal Stress Change in Fault Max( CFF) Year
0.04
0.03
0.02
0.01
0.00
Change (MPa)
(a)
1500 1000 500 0 500 1000 1500
Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Z-Coordinate (m)
Average Thermal Stress Change in Fault at End of Simulation
0.04
0.02
0.00
0.02
0.04
0.06
0.08
Change (MPa)
(c)
1500 1000 500 0 500 1000 1500
Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Z-Coordinate (m)
Standard Deviation of Thermal Stress Change in Fault at End of Simulation
0.1
0.2
0.3
0.4
0.5
Std. dev. (MPa)
1500 1000 500 0 500 1000 1500
Y-Coordinate (m)
500
1000
1500
2000
2500
3000
3500
4000
Z-Coordinate (m)
Standard Deviation of Thermal Stress Change in Fault Max( CFF) Year
0.025
0.050
0.075
0.100
0.125
0.150
0.175
0.200
Std. dev. (MPa)
(d)
Figure 3.15: All changes shown are deviations from initial conditions. (a) Ensemble average thermal
stresses corresponding to Fig. 3.14(a). (b) Standard deviation of the ensemble thermal stresses in (a). (c)
Ensemble average thermal stresses corresponding to Fig. 3.14(c). (d) Standard deviation corresponding to
(c).
63
(b)
(a)
0 5 10 15 20 25 30 35
Thermal Stress or Overpressure (MPa)
0
50
100
150
200
Count
Bottom of Fault (Y = 0.0 m, Z = 500.0 m)
0 1 2 3 4 5
Thermal Stress or Overpressure (MPa)
0
20
40
60
80
100
120
Count
Top of Fault (Y = 0.0 m, Z = 4050.0 m)
Thermal Stress, End of Simulation
Overpressure, Year of Max( CFF)
Figure 3.16: Comparison of thermal stress and overpressure at the top and bottom of the fault. (a) The
thermal stress (positive in tension) and the overpressure in the fault relative to the initial conditions at the
top of the fault. (b) The thermal stress and the overpressure in the fault relative to the initial conditions at
the bottom of the fault but slightly above the storage (deep) aquifer.
64
3.4 Discussion
3.4.1 SolubilityintheFault-LeakageContext
Solubility in GCS has been frequently studied within the context of a CO
2
storage reservoir, but this
study’s novelty includes examining the effects of solubility in the context of fault leakage. Extensive work
has been done on solutal convection and density-driven fingering (Amooie et al., 2018; Ershadnia et al.,
2020; Singh et al., 2021). The model problem studied in these papers has a CO
2
plume positioned verti-
cally above a brine column and migrating laterally beneath the bottom boundary of a caprock layer while
gradually dissolving in brine. CO
2
dissolution causes density instabilities to form, causing fingering and
convection due to the sinking of denser, CO
2
-laden brine into the relatively less dense brine below. The
consequence of such convection is increased solubility trapping in the storage reservoir. In contrast, solu-
bility during fault-leakage is primarily a matter of CO
2
leakage-rate effects: as the fault itself would not be
expected to store much CO
2
, the effect of solubility manifests in the magnitude of free-phase CO
2
leakage.
Figure 3.9 shows that brine leakage rates do not change much with solubility, although they do exhibit a
slight shift toward lower leakage rates while free-phase CO
2
leakage rates decline substantially. Here, the
flow in the fault is driven by the pressure gradient between the two aquifers and buoyant force due to the
density difference between the CO
2
and brine phases. Nevertheless, CO
2
dissolution does occur within
the fault and the storage reservoir. Unlike previous studies, which focused on the effects of solubility in
the storage reservoir, this study shows an initial estimate of the effect of CO
2
solubility on the magnitude
of free-phase CO
2
and brine leakage.
3.4.2 GeneralizabilityandLimitationsforRealGCSSites
This study definitively shows that thermal stress changes can rival pore pressure increases in faults
during fault leakage if there are phase change-induced temperature changes in the fault, a finding which
impacts site evaluation for operators. Recent studies have incorporated thermal effects in the case of fault
65
leakage behavior (N.-H. Jung et al., 2015) and noted the importance of stress fields in affecting failure, which
then affects leakage pathways (Miocic et al., 2020). This work connects these two lines of investigation
by showing that thermal stresses resulting from phase change-induced temperature changes during fault
leakage during GCS can affect stress fields enough (on the order of 1 MPa - see Figure 3.16 and compare
to a threshold of 0.1 bar (Reasenberg & Simpson, 1992)) to induce fault failure and so potentially affect the
leakage behavior. As the changes in thermal stresses would tend to be concentrated near the phase-change
boundary from liquid to gaseous CO
2
, operators should estimate the phase-transition depth and use it to
enhance their fault-seal analysis.
Given the simplified fault model, there are limitations to this study’s generalizability, especially with
respect to its modeling of a fault seal. Fault zones at sites considered for GCS can be complex (Dockrill &
Shipton, 2010; Mulrooney et al., 2020), and the sealing capacity of a fault can vary widely depending on
factors such as the fault’s throw (Wu et al., 2021), fault-rock strength parameters (Rahman et al., 2021), and
in-situ stresses (Rahman et al., 2021). Thermal stresses induced by CO
2
phase change can also vary with
the rock properties, such as the thermal expansion coefficient or drained bulk modulus (Eqn. 3.2). This
study’s model also does not take into account differences in leakage due to dip angle, which have been
found to be important (Meguerdijian & Jha, 2021). While CO
2
migration parallel to the fault is taken into
account in this study, the use of a fault seal analysis algorithm (e.g., the shale-gouge ratio (SGR) (Yielding
et al., 1997)) and the incorporation of dip angle would be helpful in a site study.
3.5 Conclusions
Thermal stresses and CO
2
solubility are phenomena that can affect the fault stability and leakage rates
during CO
2
leakage from a deep saline aquifer. In this study, we find that:
1. Solubility effects can be neglected insofar as estimating poroelastic fault reactivation or brine leakage
rates are concerned, but solubility cannot be neglected for estimating free-phase CO
2
leakage rates.
66
2. Thermal stresses can be more significant at shallow depths than at deep depths in a fault when CO
2
phase change can take place, unlike pore pressure changes, which are larger at deep depths in the
fault.
3. Increases in tensile thermal stress post-injection can be sufficient to cause failure in the fault when
phase change is present.
4. Given the possibility of CO
2
phase change during fault leakage and the associated thermal stresses,
it may be advantageous to estimate the phase change path of potentially leaking CO
2
and use it as
a variable or parameter in fault-leakage ROMs.
We show that pore pressure and thermal stresses exhibit distinct effects in different locations but that
both can potentially cause fault reactivation. Future work involves investigation of dip-angle effects on
leakage rates and the creation of fault-leakage ROMs.
67
Chapter4
Physics-InformedMachineLearningforFault-LeakageReduced-Order
Modeling
Abstract
Geologic carbon storage (GCS) is a promising technology for mitigating CO
2
emissions. The overall
success of GCS depends on safe operations that are informed by risk assessment and have proper mitiga-
tion plans in place. Performing quantitative probabilistic risk assessment for a GCS site using traditional
reservoir simulators can be challenging due to the high computational costs. To overcome this challenge,
the US Department of Energy’s National Risk Assessment Partnership (NRAP) project has developed an
integrated assessment modeling approach that utilizes computationally efficient reduced-order models
(ROM) for simulating various parts of a GCS storage site to quantify uncertainty. In this study, we de-
velop a reduced-order model for fault leakage risk assessment. We use a deep learning approach to build
the reduced-order model. We perform a sensitivity analysis and find that the deep learning model yields
high accuracy with a much smaller computational cost than full-physics simulation. We also evaluate the
performance of the model in scenarios where simulations are not possible to run, providing analysis not
previously performed in fault-leakage ROM analyses. Based on a sensitivity analysis of the model, we
suggest a simplified conceptual model for fault leakage and site monitoring.
68
Keywords: geologic carbon storage, fault leakage, reduced-order modeling, deep learning, sensitivity
analysis
4.1 Introduction
Climate change due to rapidly increasing anthropogenic greenhouse gas emissions has become a ma-
jor environmental issue over past decades. Geologic Carbon Storage (GCS) has been proposed as a ma-
jor climate change mitigation strategy (IPCC, 2005). In the US, interest in deploying GCS technology is
growing in part to take advantage of recently-expanded §45Q tax credits (26 USC §45Q, 2019). While a
well-managed GCS site has minimal risks from injection, there are risks associated with GCS, including
seismicity (Dana et al., 2022; Jagalur-Mohan et al., 2018; Jha & Juanes, 2014; White & Foxall, 2016), ground-
water impacts (Keating et al., 2010), leakage (Newell & Ilgen, 2019), and their interactions (Guglielmi et al.,
2017; Guglielmi et al., 2021; Meguerdijian & Jha, 2021; Meguerdijian et al., 2022; Newell & Ilgen, 2019).
Consequently, comprehensive risk assessment and risk management/mitigation are important to increase
stakeholder confidence in GCS technology and make its large-scale deployment viable.
Risk assessment requires taking into account a multitude of phenomena and couplings. These phenom-
ena include: fluid flow, thermal effects, chemical effects, geomechanics, and their various couplings (Newell
& Ilgen, 2019). Evaluating these phenomena and couplings is often performed through full-physics nu-
merical simulations because numerical simulators are able to incorporate the various sources of data and
physical governing equations into one comprehensive model. However, these numerical simulations can
be computationally expensive, can suffer from strong nonlinearities preventing simulation completion, and
can require large numbers of runs (i.e. ensembles) to quantify the uncertainty of risks. Consequently, many
studies have turned to reduced-order models (ROM), also known as surrogate modeling in this context, to
reduce the computational cost. These include: low-fidelity modeling with linking functions for modeling
69
groundwater impacts (Bianchi et al., 2016); a ROM for CO
2
impacts on shallow, unconfined aquifers (Keat-
ing et al., 2016); a wellbore leakage ROM (D. R. Harp et al., 2016); and integrated assessment models linking
various ROMs to rapidly assess the complexity of an entire GCS operation for use in screening potential
GCS sites (Middleton et al., 2020; Pawar et al., 2016; Stauffer et al., 2020; Vasylkivska et al., 2021). These
types of ROMs have been applied to field cases (Onishi et al., 2019).
Leakage of brine and CO
2
through faults is one of the important risks to characterize and assess but is
also challenging (Guglielmi et al., 2017; Guglielmi et al., 2021; Meguerdijian & Jha, 2021; Miocic et al., 2019;
Rinaldi et al., 2014; Snippe et al., 2022; F. Zheng et al., 2021) and consequently has had attention for ROM
development. Multiple studies have discussed challenges regarding fault zones and their characterization,
including with respect to leakage (Johri, 2012; Scholz, 2018; White & Foxall, 2016). The fault-leakage ROMs
seeking to quantify these risks have included a ROM based on a 2D model (Lu et al., 2012), fault slip with
geomechanics (White et al., 2017), and a reduced physics model of the fault (Lindner, 2020; Vasylkivska
et al., 2022). Some field trials of fault leakage are currently being performed and could potentially be used
for validation (cf. “The CSIRO," n.d.) and others have yielded results (Guglielmi et al., 2017; Guglielmi et
al., 2021). The state-of-the-art has advanced so as to allow for improved and rapid risk assessment using
fault-leakage ROMs.
Several challenges exist with respect to building and using fault-leakage ROMs. These challenges fall
into three categories: (1) variability of potential field scenarios, (2) estimating parameters of interest for
these scenarios, and (3) lack of field data. Regarding the variability of potential leakage scenarios, there
are many different types of fault settings. There can be site-level and basin-scale interactions between
injection sites, and nearby wells and faults can also potentially interact with one another, complicating
the resulting leakage behavior. There are also numerous possible arrangements of storage, groundwater,
and drinking-water aquifers with depth, with leakage into one or more shallower aquifers from a storage
70
aquifer potentially being possible, with CO
2
leakage downward being an unusual but also possible out-
come (Meguerdijian & Jha, 2021). Different types of physics, such as phase change or CO
2
dissolution, may
occur. The fault itself may have a fault core and damage zone with different properties (Scholz, 2018). In
addition to all these, each part of the leakage system - the storage aquifer, shallow formations or ground-
water aquifers into which leakage can occur, and the faults themselves can have heterogeneity in their
properties. Even parts under operational control, such as the wells, can have potentially varying rates,
pressures, or other settings over time (F. Zheng et al., 2021). Taking into consideration these types of un-
certainties and variability requires exploring a large number of potential scenarios. Regarding estimating
parameters of interest for these potential scenarios, the parameter space from which to sample becomes
very large due to the number of variables and their possible values and combinations (cf. Razavi et al. 2012
for a discussion of the challenges of high-dimensional problems). This makes it necessary to reduce the
computational cost even for ROM generation, let alone uncertainty quantification. Simulators may expend
large computational cost for scenarios with strong nonlinearities, and in many scenarios, simulators may
simply fail to simulate the chosen scenario because the physics are too nonlinear for the simulator’s solver
to handle, with no available solver on hand being able to solve the problem. This leads to gaps in the
sample of the parameter space. However, simulators may also be restricted in the physics they can simu-
late. These all limit the ability to numerically study fault leakage for ROM development. Lack of field data
for fault leakage is and will be one of the major challenges since the goal of GCS operations is to ensure
no fault leakage and to-date few targeted field experiments inducing fault leakage have been performed.
The limited amount of applicable data for fault leakage shows large variations in fault permeability that
are non-monotonic in time (Guglielmi et al., 2021). The above-mentioned limitations imply that it will be
challenging to bound the uncertainties in identifying important parameters and potential scenarios.
In spite of these limitations and challenges, it is possible to develop a fault-leakage ROM that can be
used for risk assessment by taking into consideration the following:
71
1. The fault-leakage ROM is based on a conceptual model that concisely captures fault-leakage com-
plexities
2. The fault-leakage ROM captures the complexity of the underlying physical processes by reproducing
a full-physics simulator’s results exactly (cf. the comment in Razavi et al., 2012 about the ideal ROM
or surrogate model for deterministic simulators being an interpolating model)
3. The fault-leakage ROM has a minuscule computational burden compared to the full-physics simu-
lator
4. The fault-leakage ROM yields quantifiably trustworthy results where neither full-physics simulator
results nor field data is available
5. The fault-leakage ROM can be easily integrated into an integrated assessment model
In this paper, we build and evaluate a fault-leakage ROM and extract insights for field operations. We
use a conceptual model which accounts for variations in the properties of the fault, a storage aquifer, and
an offset aquifer. We focus on deep saline aquifers. We use a deep learning approach to develop a fault-
leakage ROM to obtain both accuracy and reduced computational cost. We also consider the performance
of the ROM over missing data and show the workflow for doing so, where missing data is data at time
points that were unable to be simulated. A sensitivity analysis of the ROM is performed, and a conceptual
model and monitoring recommendations based on the ROM’s key parameters are suggested.
72
4.2 Methods
4.2.1 DescriptionoftheMeshandNumericalModel
We perform 3D multiphase simulations of fluid flow and heat transfer of brine and CO
2
leakage through
a fault using the FEHM simulator (Zyvoloski et al., 2015; Zyvoloski et al., 1999). We consider fluid me-
chanics and thermal effects. Geomechanics or CO
2
solubility is not included in this study. We follow an
observed case of CO
2
leakage (see Figure 4.1 and Miocic et al., 2019). The domain consists of two aquifers
(a storage reservoir and a shallower leakage aquifer) connected by a fault penetrating through a layer of
impermeable caprock (see Figure 4.2). The mesh is unstructured and consists of triangular-prism elements.
We generated the meshes for varying fault angles using LaGriT (Gable et al., 2021). The fault consists of
hexahedral elements inserted into the unstructured mesh, resulting in a fixed fault thickness of 3 m. Higher
refinement is included in the fault to more accurately model flow. Additional refinement is also added at
the tops of the aquifers to improve the accuracy of plume movement near the fault. See Figure 4.3 for a
sample mesh. We use an unstructured mesh with smooth faults to ensure continuous fluid flow though the
fault as observed in the field. Our earlier attempts with a structured mesh resulted in artificial blocking of
upward CO
2
migration due to the right angle elements in the structured mesh.
To generate a suitable dataset for training the fault-leakage ROM, we run an ensemble of simulations.
The inputs to the simulations are varied using Latin Hypercube Sampling (Helton & Davis, 2002; McKay et
al., 1979) according to the parameters in Table 4.1. The fixed parameters - those that do not change between
simulations in the ensemble - are shown in Table 4.2. The dip angle varies between 40
◦ and 140
◦ (Table 4.1),
with seven different meshes being used to account for this variation. This is done because a previous study
showed the effect of dip angle on leakage (Meguerdijian & Jha, 2021). The caprock is adjusted in all dip-
angle cases to proceed 100 m from the centerline of the fault, extending outwards perpendicularly to the
fault. The well locations in Table 4.1 correspond to distances from the intersection of the fault with the top
73
of the deep aquifer to the well in the deep aquifer. The well locations correspond to distances from this
intersection roughly as follows: Location 0 = 200 m, Location 1 = 400 m, Location 2 = 600 m. We sample
certain parameters discretely and others continuously, as shown in Table 4.1. We use the Model Analysis
Toolkit (MATK) (D. Harp et al., 2013) to generate the ensemble and run the simulations using FEHM.
We invoke symmetry conditions at thexz-origin plane. Injection is performed within this plane. Con-
sequently, the simulation yields results for an aquifer twice the size of that which is simulated. The injec-
tion rate used in this study corresponds to one-half the injection rate into this doubly-sized aquifer.
St. Johns Dome
CWF
Travertines Travertines
GWC today
Paleo
GWC
Buttes
Fault
Reservoir
CO
2
E
CO
2
W
Figure 4.1: Physical setting of fault leakage. This figure is a modified version of Figure 1(b) from (Miocic
et al., 2019). Used under the Creative Commons Attribution 4.0 International License. To view a copy of
this license, please visit http://creativecommons.org/licenses/by/4.0/.
θ
Deep Aquifer
Shallow Aquifer
Fault (3 m)
Caprock
Caprock
200 m
200 m
100 m
4000 m
z
x
100 m
4000 m
1200 m in the y-direction (into the page)
Symmetry BC along x-z plane
Open side boundary
(constant at P
initial
)
P
top
= 8.15 MPa, T
top
= 50
o
C
CO
2
Injection Well
1750 m
z
Figure 4.2: Schematic of the model problem. “BC” denotes boundary condition.
74
Fault
Caprock
Figure 4.3: An example of the unstructured meshes used in these simulations. The dip angle (angle from
the horizontal) is 60
◦ . Note the additional refinement at the tops of the two aquifers near the fault to more
accurately capture CO
2
plume movement and the high refinement of the fault to model flow accurately.
Note also that the zoom-in figure of the fault has a different color scheme to differentiate the fault from
the surrounding caprock.
4.2.2 ROMDevelopmentWorkflow
We develop a ROM using a two-layer deep neural network (DNN). A schematic of the DNN is shown in
Figure 4.4. We use rectified linear unit (ReLU) functions as the activation functions because of their ability
to model nonlinear behaviors and known qualities for deep learning. We use a dense neural network with
a linear output layer. There are two separate DNNs: one for CO
2
and one for brine. These DNNs jointly
compose the ROM. We use the TensorFlow framework to build the ROMs (Martín Abadi et al., 2015) and
use L
1
kernel regularization with weights of 10
− 4
for each hidden layer. The loss function used is mean
absolute error, and the optimizer used is Adam (Kingma & Ba, 2017). All inputs to the network were scaled
to the range of [0,1] using Min-Max scaling available within the scikit-learn (Pedregosa et al., 2011) library.
We structure the ROM to predict brine and CO
2
leakage rates at individual points in time as opposed
to whole time sequences. As shown in Figure 4.5, leakage profiles of total brine and CO
2
leakage across
75
Figure 4.4: Schematic of the deep neural network used in this study.
the fault’s entrance to the shallow aquifer are obtained from FEHM simulation results. These leakage
profiles are subsequently split into training and test sets i.e. one set is used for training the DNN and the
second set is used for testing the resulting DNN. Each of these sets of leakage profiles are then broken
down into individual time points e.g. each leakage-rate time series of 101 points (for years 0-100 in the
simulation) results in 101 samples. The ROM sets the leakage rates at time t=0 years to 0.0 t/yr because no
leakage occurs before injection begins. The values from a given simulation’s leakage profile (time series)
are never used to predict other values in the same time series to avoid trivializing the ROM model into
a nearest-neighbor interpolation. Note that this method of predicting individual points in time contrasts
with methods involving recurrent neural networks, such as long short-term memory (LSTM) (Hochreiter
& Schmidhuber, 1997), which predict a sequence of outputs (e.g. a time series) based on an initial input
and changes in the network’s internal states. We adopt this atypical approach to build ROMs that function
like analytical solutions to a system of PDEs, i.e. the ROM should be able to output leakage values at any
point in time without the need to predict at other points in time. Furthermore, predicting in this manner
increases the number of training samples available for our ROM by orders of magnitude, allowing for the
training of larger and hopefully more accurate networks (consider thousands of simulations vs. hundreds
of thousands of time points).
76
Figure 4.5: Workflow for the generation of fault-leakage reduced-order models.
Variable Value
Damage zone permeability (log
10
m
2
) [-15, -12]
Damage zone porosity [0.001, 0.1]
Shallow aquifer permeability [-14, -12]
Shallow aquifer porosity [0.05, 0.5]
Deep aquifer permeability [-14, -12]
Deep aquifer porosity [0.05, 0.35]
Well location
∗ {0,1,2}
Injection rate
†
(kg/s) [0.25, 12.5]
Dip angle (
◦ ) {40, 60, 80, 90, 100, 120, 140}
Injection duration (years) [10,50]
Geothermal gradient (
◦ C/km) [8,44]
Table 4.1: The varying parameters of the model. Note that [] correspond to continuous ranges and {}
correspond to sets of discrete values.
∗ These correspond to distances from the well, with larger values
indicating larger distances. Because the meshes are unstructured, there is some variation in the distance
between dip angles. However, each well location value corresponds roughly to a distance of 200 m from
the intersection of the deep aquifer with the fault.
†
These injection rates are for the half-reservoir due to
the invocation of symmetry along the injection well.
77
Parameters Value
Deep aquifer thickness (m) 200
Caprock thickness (m) 100
Shallow aquifer thickness (m) 200
Density of the rock (solid phase) (kg/m
3
) 2650
Specific heat capacity of the rock (solid phase) (J/(kg · K)) 920
Surface temperature (at the top of the domain)(
◦ C) 50
Initial pressure at the top of the domain (MPa) 8.15
Porosity of the caprock (dimensionless) 0.1
Permeability of the caprock (m
2
) 10
− 18
Thermal conductivity of rock (W/(m· K)) 3
CO
2
injection temperature (
◦ C) 32
Height of the domain (m) 1750
y-width of the domain (m) 1200
Wellz-coordinate (m) 178-190
Total simulation time (years) 100
Relative permeability model linear
Table 4.2: Fixed parameters of the reservoir model.
4.3 Results
4.3.1 BaseCase
We begin the study by performing a base case simulation to demonstrate the physics inherent in this
fault-leakage problem. The varying parameters for this base case are shown in Table 4.3. The permeability
distribution is shown in Figure 4.6. The fault (labeled as the damage zone in Table 4.3) is moderately
permeable - about 80 millidarcy (mD). The deep aquifer is highly permeable, with a permeability of 1
darcy. The shallow aquifer has a permeability of about 25 mD. The aquifers are highly porous, while the
fault is somewhat porous (Figure 4.6). The injection rate of CO
2
is equivalent to 303,000 t/yr of CO
2
into
the aquifer (including the symmetry portion not simulated). The overpressure distribution is shown in
Figure 4.7, which shows a decline in pressure in the fault a few years after the end of injection. Figure 4.8
shows cooling around the well due to the injection of CO
2
cooler than the surrounding rock and fluid in
the deep aquifer. Figure 4.9(a) shows the CO
2
plumes in the deep and shallow aquifers. The CO
2
spreads
out in the deep aquifer and migrates both toward and away from the fault. The CO
2
that reaches the fault
78
Variable Value
Damage zone permeability (log
10
m
2
) -13.1
Damage zone porosity 0.09
Shallow aquifer permeability (log
10
m
2
) -13.6
Shallow aquifer porosity 0.3
Deep aquifer permeability (log
10
m
2
) -12
Deep aquifer porosity 0.23
Well location 0
Injection rate (kg/s) 4.8
Dip angle (
o
) 60
Injection duration (years) 29
Geothermal gradient (
o
C/km) 20.7
Table 4.3: The varying parameters of the base case.
Figure 4.6: The permeability and porosity distributions of the base case. Note that the fault appears as a
thin line.
then migrates upward and spreads across the shallow aquifer. The loss in resolution far from the shallow
aquifer is permissible because the emphasis in this study is on flow through and around the fault, and
the resolution of saturations further from the fault is consequently not important. Figure 4.9(b) shows
that water is drained from the locations where CO
2
enters the pore space. Note that the CO
2
and water
saturations appear more uniform in the shallow aquifer, potentially due to lower permeabilities in the
shallow aquifer promoting more even drainage than otherwise, and the open pressure boundary promoting
plume movement. These figures show the primary physics of the model problem, that is, multiphase fluid
movement in and near the fault.
79
Figure 4.7: The overpressure distribution of the base case at year 32.
Figure 4.8: The temperature distribution of the base case at year 32.
80
(a)
(b)
Figure 4.9: The base case’s (a) supercritical CO
2
saturation, and (b) water saturation at year 32 in the
simulation.
81
4.3.2 ROMPerformance
We train the ROM to nearly memorize the data and yet generalize well and with substantial speedups
over the simulator. We calculate brine and CO
2
leakage rates by aggregating leakage fluxes across the
fault interface with the leakage aquifer. The training accuracies of the leakage-rate predictions are near-
perfect (Figure 4.10(a) and (c)). The test accuracies also show excellent generalization, with the predicted
brine-leakage rates being close to the R
2
=1 line (Figure 4.10(b) and (d)), where R
2
is the coefficient of
determination (Scikit-learn User Guide 3.3.4.1, Pedregosa et al., 2011) and is defined as:
R
2
(y,ˆ y)=1− P
n
i=1
(y
i
− ˆ y
i
)
2
P
n
i=1
(y
i
− ¯y
i
)
2
(4.1)
where y is the true sample value, ˆ y is the predicted sample value, y
i
is the true i-th sample value, ˆ y
i
is
the predictedi-th sample value,n is the number of samples, and ¯y is the arithmetic mean of true sample
values (
1
n
Σ n
i=1
y
i
).
The CO
2
test set predictions show greater spread than the brine test set predictions. Figure 4.10(b)
shows more spread than Figure 4.10(d). Figure 4.11 confirms this: while both brine and CO
2
leakage
have prediction errors clustered around zero (as shown by the locations of the peaks), the extremes of
CO
2
leakage prediction errors are substantially further from zero than those of brine leakage prediction.
Still, the errors are relatively small, and the speedup from using the ROM in place of the simulator is
considerable: the average time for a single simulation of a given scenario is 6148 seconds on average using
64 processors and 264 GB RAM to simulate the ensemble, while obtaining the leakage profiles for the same
case requires only 0.02 seconds on average, an over 99.99% decrease. Consequently, the ROM achieves its
intended purpose reasonably well: predicting instantaneous leakage rates quickly.
The ROM also performs well according to a number of other measures. Cumulative leakage is an
important measure of performance because the concern with CO
2
leakage is not just with the leakage
82
rate, but the leakage amount (which is what ultimately affects the CO
2
concentration in the atmosphere).
Cumulative leakage of brine can also affect the salinity of shallow aquifers into which leakage occurs.
Consequently, we look at cumulative leakage amounts. We see in Figure 4.12 shows that the cumulative
fits for both CO
2
leakage (Figure 4.12(a)) and brine leakage (Figure 4.12(b)) are close to R
2
=1. Another
measure, the onset time of CO
2
leakage, measures the year in which the CO
2
leakage rate into the shallow
aquifer exceeds a given threshold. We varied the threshold value to see how well the ROM could perform
with different thresholds. With a threshold of 275 t/yr (chosen because of its high accuracy compared to
other threshold values), the ROM is able to predict the onset time fairly well, although not as accurately
as leakage rates (Figure 4.13). We note that this measure’s R
2
value is sensitive to the chosen leakage-
rate threshold. The maximum leakage rate is a measure that has previously been used to quantify leakage
risk (D. R. Harp et al., 2016) and can be used as estimate of maximum leakage risk. We see that the R
2
values
of the maximum brine and CO
2
leakage rates are predicted with a high degree of accuracy (Figure 4.14).
The maximum rates are predicted by using the ROM to calculate the brine and CO
2
leakage rate at years
0-100 for each set of static parameters and taking the maximum. The maximum outflow rate of brine from
the shallow aquifer is calculated similarly and yields acceptable accuracy (Figure 4.15). Consequently, we
see that the ROM preserves time series statistics despite each point in the time series being predicted
separately.
We train the ROM until further improvement in the test accuracy is not obtained. We train for 1000
epochs regardless of the test error (i.e. the test set is not used as a validation set to tune hyperparameters)
(Figure 4.16). We note that the training error continues to go down while the test error stabilizes. The
use of regularization likely helps the test error stay nearly constant even as the training error continues
to decline. Because the test error stabilizes during training, the ROM can be considered fully trained.
83
(a)
(b)
(c)
(d)
Figure 4.10: The performance of the ROM for predicting (a) instantaneous CO
2
leakage rate on the training
set, (b) instantaneous CO
2
leakage rate on the test set, (c) instantaneous brine leakage rate on the training
set, and (d) instantaneous brine leakage rate on the test set.
4.3.3 SensitivitytoInputParameters
4.3.3.1 DeltaSensitivityAnalysis
With respect to the importance of the varying parameters in the ROM’s calculations vs. in the FEHM
simulations, we note that the ROM considers a few varying parameters important while FEHM considers all
the varying parameters important. We perform a delta moment-independent measure sensitivity analysis
84
6 4 2 0 2 4 6
Prediction Error (Predicted-Actual, (t/yr)) x10
5
CO
2
Brine
Figure 4.11: The prediction error distributions of the instantaneous CO
2
leakage rate for the test set and
the instantaneous brine leakage rate for the test set. The prediction error is the predicted value of the
leakage rate minus the actual value from the simulation.
using the SALib package (Borgonovo, 2007; Herman & Usher, 2017; Plischke et al., 2013). The delta measure
is defined as (Borgonovo, 2007; Plischke et al., 2013):
δ i
=
1
2
E
X
i
[s(X
i
)] (4.2)
E
X
i
[s(X
i
)]=
Z
x
i
f
X
i
(x
i
)s(x
i
)dx
i
(4.3)
wheres(x
i
)=
R
y
i
|f
Y
(y)− f
Y|X
i
=x
i
(y)|dy,X
i
an input parameter which is assumed to be at a fixed value,
f
Y
(y) is probability density function of the output (in this case, brine or CO
2
leakage rates),f
X
i
(x
i
) is the
marginal density ofx
i
,f
Y|X
i
(y) is the distribution of outputs assuming that inputX
i
is fixed. Essentially,
δ i
measures the shift in the output distribution when parameter X
i
is held constant. The value of δ i
is
not dependent on any particular moment of the distribution and is invariant to monotonic transformation.
85
0 2 4 6
Actual Leakage Amount (t) 1e7
0
2
4
6
Predicted Leakage Amount (t)
1e7
Test Set - Cumulative CO
2
Leakage at Year 100
R
2
= 0.966
(a)
0 2 4
Actual Leakage Amount (t) 1e7
0
1
2
3
4
5
Predicted Leakage Amount (t)
1e7
Test Set - Cumulative Brine Leakage at Year 100
R
2
= 0.975
(b)
Figure 4.12: The performance of the ROM for predicting the following at 100 years: (a) cumulative CO
2
leakage on the test set, and (b) cumulative brine leakage on the test set.
We calculate the δ i
values for the varying parameters for both the ROM and the FEHM simulations, for
both brine and CO
2
leakage. We also calculate theδ i
value for a dummy variable that is random and has
no effect on the leakage amounts to estimate the δ i
value attainable via approximation error. δ i
values
would need to exceed these dummy values to be considered significant. We use bootstrapping to estimate
confidence intervals for δ i
.
86
0 20 40 60 80 100
Actual Onset Time (yr)
0
20
40
60
80
100
Predicted Onset Time (yr)
Test Set - Onset Time of CO
2
Leakage
R
2
= 0.645
Figure 4.13: The performance of the ROM for predicting the onset time of CO
2
leakage. Note that the
threshold for determining the onset of CO
2
leakage is 275 t/yr.
0.0 0.5 1.0 1.5 2.0 2.5
Actual Leakage Rate (t/yr) 1e6
0.0
0.5
1.0
1.5
2.0
2.5
Predicted Leakage Rate (t/yr)
1e6
Test Set - Maximum Rate of Brine Leakage
R
2
= 0.996
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
Actual Leakage Rate (t/yr) 1e6
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
Predicted Leakage Rate (t/yr)
1e6
Test Set - Maximum Rate of CO
2
Leakage
R
2
= 0.972
(a) (b)
Figure 4.14: The performance of the ROM for predicting (a) the maximum CO
2
leakage rate and (b) the
maximum brine leakage rate.
The sensitivity analysis indicates that the ROM achieves high accuracies by focusing on a few key
variables. In Figure 4.17(a), the sensitivity analysis for the ROM considers only the following variables sig-
nificant for CO
2
leakage: damage zone permeability, deep aquifer permeability (though barely), injection
rate, dip angle, injection duration, and time. In contrast, the sensitivity analysis for the FEHM simula-
tions shows that all the varying input parameters are significant. The sensitivity analysis for brine leakage
87
0 1 2 3 4 5 6
Actual Outflow Rate (t/yr) 1e5
0
1
2
3
4
5
6
Predicted Outflow Rate (t/yr)
1e5
T est Set - Maximum Brine Outflow Rate From the Shallow Aquifer
into the Fault
R
2
= 0.946
Figure 4.15: The performance of the ROM for predicting the maximum brine outflow rate from the shallow
aquifer.
shows similar results in Figure 4.17(b): the ROM’s significant parameters are determined to be damage
zone permeability, deep aquifer permeability, injection rate, dip angle, injection duration, and time, but all
the FEHM simulations’ varying parameters are considered significant. Unlike the FEHM simulations, the
ROM achieves high accuracies by primarily considering only a few variables that are key to the physics
of leakage. The reason why the ROM focuses on a few key parameters while FEHM is sensitive to all the
parameters may be that the ROM’s network weights are L
1
-regularized, promoting sparsity in the network
weights and, possibly as a consequence, sparsity in its sensitivity to inputs.
4.3.3.2 Sobol’SensitivityAnalysis
We use other global sensitivity analysis methods, such as the Sobol’ sensitivity analysis, to confirm
the delta sensitivity analysis’s determination of the ROM’s significant parameters. The Sobol’ sensitivity
estimates are based on the proportion of variance explainable by a particular parameter directly or with
its interactions with other parameters (Homma & Saltelli, 1996; Saltelli, 2002a, 2002b; Saltelli et al., 2010;
88
0 200 400 600 800 1000
Epochs
0
10000
20000
30000
40000
50000
60000
Mean Absolute Error (t/yr)
CO
2
Train and Test Error vs. Number of Training Epochs
Train
T est
(a)
0 200 400 600 800 1000
Epochs
0
10000
20000
30000
40000
50000
60000
70000
Mean Absolute Error (t/yr)
Brine Train and T est Error vs. Number of Training Epochs
Train
T est
(b)
Figure 4.16: The train and test errors of the ROM during training for the (a) CO
2
model, and (b) brine model.
Note that this training was performed excluding samples with time att=0 years because the leakage value
att=0 is set to 0.0 t/yr.
Sobol’, 2001). To illustrate the meaning of the Sobol’ analysis, lety =f(x
1
,x
2
,...,x
N
) be a model where
y is the output andx
i
is thei
th
input parameter. Let the functionf be decomposed as follows:
y =f(X)=f
0
+
X
i
f
i
+
X
i
X
j>i
f
ij
+...+f
12...N
(4.4)
89
Damage Zone Perm.
Damage Zone Porosity
Shallow Aquifer Perm.
Deep Aquifer Perm.
Shallow Aquifer Porosity
Deep Aquifer Porosity
Well Location
Injection Rate
Dip Angle
Injection Duration
Geothermal Gradient
Time
0.00
0.05
0.10
0.15
0.20
FEHM Dummy Line
ROM Dummy Line
ROM
FEHM
(a)
Damage Zone Perm.
Damage Zone Porosity
Shallow Aquifer Perm.
Deep Aquifer Perm.
Shallow Aquifer Porosity
Deep Aquifer Porosity
Well Location
Injection Rate
Dip Angle
Injection Duration
Geothermal Gradient
Time
0.00
0.05
0.10
0.15
0.20
0.25
0.30
FEHM Dummy Line
ROM Dummy Line
ROM
FEHM
(b)
Figure 4.17: Delta moment-independent measure values for (a) CO
2
-leakage DNN ROM and FEHM sim-
ulations and (b) brine-leakage DNN ROM and FEHM simulations. The significant variables for the ROM
and FEHM simulations are those above their respective dummy lines. The small black lines at the tops of
the columns denote the 95% confidence intervals for the δ values obtained from bootstrapping.
90
where each individual termf
ij...k
is a function only of the parametersx
i
,x
j
,...,x
k
. The interaction terms
are those that are a function of more than one variable. If the input parameters are independent and the
individual termsf
ij...k
are square integrable over the domain, then a unique decomposition of the variance
exists, which is:
V(y)=
X
i
V
i
+
X
i
X
j>i
V
ij
+...+V
12...N
(4.5)
whereV
i
=V(E(y|x
i
)),V
ij
=V(E(y|x
i
,x
j
))− V
i
− V
j
. The additional variance terms follow a similar
pattern and reflect the variance due to each input parameter x
i
both individually and via interactions with
other parameters. We then define the Sobol’ sensitivity parameters as follows:
S1
i
=
V
i
V(y)
(4.6)
ST
i
=
V(y)− V(E(y|x
− i
))
V(y)
(4.7)
where S1
i
is the first-order Sobol’ sensitivity index indicating the first-order effect of parameter x
i
on
the variance,ST
i
is the total-order Sobol’ sensitivity index indicating the effect of parameter x
i
when all
interactions are included, and x
− i
represents all variables exceptx
i
. V(E(y|x
− i
)) reflects the change in
variance due to all parameters and interactions not includingx
i
. Note that
P
i
ST
i
may exceed the value
of 1 due to the repetition of interaction terms in multipleST
i
values. Turning to the results from the Sobol’
sensitivity analysis, we find that the ROM’s first-order and total-order Sobol’ sensitivity indices for both
brine and CO
2
leakage rates yield the same five most important parameters (Figure 4.18): time, damage
zone permeability, injection rate, injection duration, and dip angle. Note that the black bars in these figures
denote bootstrapped confidence intervals. These results agree with the results from the delta sensitivity
analysis.
91
Dip Angle
Deep Aquifer Perm.
Well Location
Damage Zone Perm.
Injection Rate
Deep Aquifer Porosity
Time
Injection Duration
Geothermal Gradient
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Fraction of Variance
Sobol' Sensitivities
Brine - First-Order
Brine - T otal-Order
CO
2
- First-Order
CO
2
- Total-Order
Figure 4.18: A Sobol’ sensitivity analysis of the ROM.
4.3.3.3 PAWNSensitivityAnalysis
We further perform the PAWN sensitivity analysis, which is a distribution-based sensitivity analysis
based on cumulative distribution functions (CDF), unlike the delta sensitivity analysis method, which is
based on probability distribution functions (PDF) (Pianosi & Wagener, 2015, 2018). The PAWN sensitivity
index for inputx
i
is as follows:
PAWN
i
= stat
x
i
max
y
|F
y
(y)− F
y|x
i
(y|x
i
)| (4.8)
where stat is a statistic (such as the mean),F
y
(y) is the unconditional CDF of the outputy, andF
y|x
i
(y|x
i
)
is the conditional CDF of y. Note that max
y
|F
y
(y)− F
y|x
i
(y|x
i
)| is the Kolmogorov-Smirnov statistic,
which is used to compare the difference between CDFs (Kolmogorov, 1933; Smirnov, 1939). In this study,
we use the median as the statistic. We use the procedure as described in SALib (Herman & Usher, 2017).
92
We use the dummy-parameter method for determining significance, like in the delta sensitivity analysis.
The results are shown in Figure 4.19. We see that the PAWN values for CO
2
leakage rates are greatest
for: time, injection rate, injection duration, dip angle, damage zone permeability, and deep aquifer per-
meability. Likewise, the PAWN values for brine leakage rates are greatest for: damage zone permeability,
time, injection duration, injection rate, dip angle, and deep aquifer permeability. Thus all four sensitivity
analyses (delta, Sobol’ first-order, Sobol’ total-order, and PAWN) agree on the same five most important
parameters for the ROM for both brine and CO
2
leakage-rate prediction, with deep aquifer permeability
as a potential sixth most important parameter.
4.3.4 PredictionofFailedRuns
We use the ROM to predict the results for a batch of early-terminating simulations and show that the
ROM is able to indicate physics leading to simulation failure. Previous work (D. R. Harp et al., 2016) showed
that incomplete simulations tend to occur due to numerical instabilities at higher values of the wellbore
permeability, which is analogous to damage zone permeability in the fault-leakage scenario. Note that Lu et
al., 2012 also showed a significant proportion of unsuccessful simulations. In this study, we also encounter
early-terminating simulations, but we now seek to expand the concept of the ROM to not only act as a
surrogate for simulation, but also as an intelligent extrapolationbeyond the feasible space of the simulator.
We propose using the ROM to explore parameter spaces that the simulator was unable to access due to
early termination. This naturally involves extrapolation, and consequently, there is a need to ensure that
the results in the extrapolated space are physically realistic. We propose assessing the accuracy of the ROM
over the time points available from the failed simulations. The goal is to evaluate if the early-terminating
simulations can be considered to come from the same distribution as the completed runs. The result of
this evaluation is shown in Figure 4.20. We see that the ROM is able to predict the available data from the
93
Damage Zone Perm.
Damage Zone Porosity
Shallow Aquifer Perm.
Deep Aquifer Perm.
Shallow Aquifer Porosity
Deep Aquifer Porosity
Well Location
Injection Rate
Dip Angle
Injection Duration
Geothermal Gradient
Time
0.000
0.025
0.050
0.075
0.100
0.125
0.150
0.175
0.200
PAWN
FEHM Dummy Line
ROM Dummy Line
ROM
FEHM
(a)
Damage Zone Perm.
Damage Zone Porosity
Shallow Aquifer Perm.
Deep Aquifer Perm.
Shallow Aquifer Porosity
Deep Aquifer Porosity
Well Location
Injection Rate
Dip Angle
Injection Duration
Geothermal Gradient
Time
0.00
0.05
0.10
0.15
0.20
0.25
PAWN
FEHM Dummy Line
ROM Dummy Line
ROM
FEHM
(b)
Figure 4.19: PAWN values for (a) CO
2
-leakage ROM and FEHM simulations and (b) brine-leakage ROM
and FEHM simulations. The significant variables for the ROM and FEHM simulations are those above their
respective dummy lines.
94
early-terminating runs as well as test data from complete runs, suggesting that they are from the same
distribution, or at least that it cannot be stated that they are from different distributions.
We also consider the ability of the ROM to predict data that is completely missing. The motivation
for performing this analysis is to better assess the ROM’s generalization. Following work done on miss-
ing data imputation (Buuren, 2018), we note that the typical approach of eliminating data entries that are
incomplete, such as the incomplete leakage time series from failed simulations, assume that data is miss-
ing completely at random (MCAR) if evaluating ROM performance using measures such as R
2
with only
complete data. This is a strategy known as listwise deletion. Multiple imputation is generally considered
a superior method to listwise deletion and consists of imputing multiple values for the missing data based
on a missing data model to form multiple completed data sets, performing the desired analysis (e.g. R
2
computation) on each completed data set, and then aggregating the analysis results. Listwise deletion is
superior to multiple imputation if: the goal is to determine regression weights (e.g. training a ROM) when
the missing data is only in the output (e.g. leakage rates), the probability of data being missing does not
depend on the output value, or the model consists of logistic regression with the probability of data being
missing depending only on the output value (Chapter 2.7, Buuren, 2018). However, multiple imputation
is a more general method than listwise deletion, and can account for relationships between inputs and
outputs as well among inputs in the presence of missing-at-random (MAR) and missing-not-at-random
(MNAR) missing data models, allowing us to estimate ROM performance beyond the data available.
We perform stochastic imputation regression to form multiple datasets of missing data estimates and
follow the workflow of multiple imputation to assess the ROM’s performance over the missing data.
Stochastic imputation regression consists of fitting a model to available data and forming a normal distri-
bution of mean zero and variance equal to the variance of the residuals. Several complete datasets can then
be generated by imputing values for each missing data point equal to the sum of the predicted regression
value and a randomly sampled value from the aforementioned normal distribution. Each of these data sets
95
are then used to calculate a measure of ROM performance (e.g. R
2
), and so a distribution of R
2
values
is obtained. In this study, we generate two imputation deep neural networks, each based on 70% of the
available test data for brine and CO
2
, both from complete and incomplete simulations. These imputation
deep neural networks have the same neural network structure as the ROM but were trained using only
100 epochs. The residuals of the imputation networks are obtained from their predictions on the remain-
ing 30% of the test data. The variances of these residuals are used to generate two normal distributions
from which to sample residuals, one for brine and another for CO
2
. Completed sets of missing data are
then generated via imputation using the imputation deep neural networks and residual sampling from the
normal distributions. This assumes that the missing data model is ignorable, which implies that the data is
missing-at-random (MAR) (i.e. the values of the missing data depend on the values of the observed data).
After generating these completed data sets, missing-not-at-random (MNAR) scaling is applied to account
for the non-randomness of the missing values. We perform MNAR scaling because the distribution of
missing data suggests that the missing data points have faster fluid flow physics (e.g. higher damage zone
permeabilities) than the completed-simulation data, which likely implies higher leakage rates in general.
MNAR scaling is performed by generating a random number between 1.0 and 2.0 (known as the MNAR
constant) and then multiplying an entire given imputed dataset realization by the MNAR constant. This
amounts to an adjustment of +0% to +100% in the leakage rates for the MNAR-adjusted cases. This follows
the recommendation in literature to account for the existence of MNAR by adding or multiplying by a con-
stant in the absence of being able to specify models that are obviously more realistic than others (Chapters
3.8.5-3.8.7 Buuren, 2018). The values of CO
2
leakage are floored to 0.0 t/yr after MNAR scaling is included
because the CO
2
leakage rates’ residual distribution showed that 0.0 t/yr was the floor for leakage rates.
We find that the ROM’s performance for brine leakage is expected to be reasonably accurate over the
range of completely missing data while the performance for CO
2
leakage is expected to degrade substan-
tially. Figure 4.21 shows the performance of the ROM over the range of missing data with MNAR scaling
96
varying from +0% to +100%. Figure 4.21(a) shows that the ROM’s performance for CO
2
is expected to
degrade from R
2
≈ 0.95 to R
2
= 0.5-0.8. Toward the lower end of R
2
, near 0.50 - which has a significant
probability of occurring - the ROM’s performance becomes questionable. In contrast, Figure 4.21(b) shows
that the ROM’s performance for brine leakage tends to be concentrated at higher values, near R
2
=0.875.
The lowest value for R
2
is near 0.70. Consequently, while the ROM is expected to experience a drop in
performance for both brine and CO
2
in the range of missing data, the ROM’s performance for brine is
expected to be acceptable while that of CO
2
may not be acceptable.
We note the limitations of stochastic regression imputation in estimating the ROM’s performance over
missing data. First, stochastic regression imputation does not take into account the uncertainty in the R
2
values due to parameter uncertainty (Chapter 3.2, Buuren, 2018). In this study, this would mean taking
into account the effect of varying the neural network parameters. We also assume a normally distributed
residual with constant residual variance over all predicted leakage ranges. This might not be the case
for brine and is known not to be the case for CO
2
because 0.0 t/yr is seen to be a floor for CO
2
leakage
rates. While this latter leakage-rate floor is taken into account, it is taken into account after calculating the
variance, which may affect the final R
2
distribution by overstating the variance. We also assume that the
non-randomness can be accounted for using a constant scaling factor varying uniformly between 1.0 and
2.0 (+0% and +100% to the leakage rate values). While this type of constant scaling is recommended by the
literature due to clarity of interpretation, it might not be representative of the variety of complex behavior
possible in the subsurface. With these limitations noted, the results still show generally how the ROM per-
formance would vary over the missing data. This work is the first time a stochastic regression imputation
workflow has been attempted for CO
2
storage. Field data for CO
2
storage technology is extremely limited,
and this work shows the potential of a stochastic regression imputation or multiple imputation approach
for addressing the lack of available field data and also the promise of transfer learning.
97
(a)
(b)
Figure 4.20: Predicted vs. actual leakage rates for the available data from incomplete simulations for (a)
CO
2
and (b) brine.
98
0.50 0.55 0.60 0.65 0.70 0.75 0.80
R
2
0
5
10
15
20
25
30
35
Frequency
(a)
0.700 0.725 0.750 0.775 0.800 0.825 0.850 0.875
R
2
0
10
20
30
40
50
60
Frequency
(b)
Figure 4.21: Stochastic regression-imputation estimates of R
2
of the ROM’s performance over missing data
for (a) CO
2
and (b) brine.
99
4.4 Discussion
4.4.1 DistributionofParametersandLeakageRatesAmongCompletevs. Incomplete
ReservoirSimulations
We analyze the parameters to elaborate on our previous comment that the incomplete or failed reser-
voir simulations typically exhibit faster fluid flow physics than the complete reservoir simulations. Fig-
ure 4.22 shows the distributions for the reservoir simulations used for the CO
2
leakage rates’ test set and
the incomplete reservoir simulations. Note that because the reservoir simulations used for the CO
2
leak-
age rates’ test set are chosen randomly from the set of completed simulations, the CDFs calculated from
this test set are representative of the entire set of completed simulations. We observe that for a given pa-
rameter, the CDF lying lower than the other exhibits a shift toward higher values. We note that the CDFs
of the incomplete runs display shifts to higher values for the following parameters: injection duration,
injection rate, and damage zone permeability. The CDFs of the incomplete runs show shifts to lower val-
ues for the following parameters: deep aquifer porosity, shallow aquifer porosity, and well location. The
dip angle exhibits a non-monotonic shift, with dip angles closer to 90
◦ being associated with completed
simulations and those further away being associated with incomplete simulations. This can be seen from
the increasing gap between the two curves until 60
◦ , then a decreasing gap and overtaking of the incom-
plete simulations’ curve by the completed simulations’ curve until 120
◦ . At 140
◦ , the gap between the
CDFs of the complete and incomplete simulations close, indicating a greater concentration of incomplete
simulations’ at 140
◦ . Analyzing these for physical insights, we note that high injection rates and injection
durations correspond to higher overpressures in the storage aquifer, with all other factors remaining the
same. Lower well location indices correspond to injection closer to the fault. Lower deep aquifer porosities
promote higher Darcy velocities. Higher damage zone permeabilities promote faster leakage through the
fault. Low shallow aquifer porosities promote faster propagation of leaked fluid in the shallow aquifer,
100
which could cause challenges at the open boundary condition. All these changes, except for that of the
dip angle, promote faster fluid flow physics and are associated with incomplete simulations, the connec-
tion being that faster fluid flow physics typically leads to nonlinear behavior that is harder to solve using
numerical approaches. In the case of dip angle, angles further from 90
o
have longer faults, which could
result in stronger nonlinearities during brine backflow after injection. However, this effect is not clear,
and further investigation is needed. Notably, this kind of analysis would aid in other ROM development
studies to identify the characteristics of incomplete simulations. These characteristics can then be used to
constrain the uncertainties regarding ROM performance over the range of incomplete simulations.
With respect to the leakage rates, we note that ROM performance can be evaluated beyond the range
of input data. The typical way of dealing with incomplete simulations (which implies missing output
data) is simply to eliminate the incomplete runs from the analysis. While this may be valid for building
a ROM (Buuren, 2018), it is not necessarily valid for testing a ROM: incomplete simulations still carry
useful information which can be used in ROM evaluation. However, because there is data missing in the
incomplete simulations, there is uncertainty regarding the data’s values. But this uncertainty is not infinite:
the missing data values are constrained by physical laws, and so the ROM’s performance over these missing
values can be quantified. By performing uncertainty quantification over missing values using a multiple
imputation workflow, the ROM’s generalizability can be better evaluated. This information can then be
given to operators to give them greater confidence regarding the ROM’s range of validity.
4.4.2 PhysicalInsights
The sensitivity analysis shows that certain parameters are influential for the ROM: simulation time,
damage zone permeability, injection rate, injection duration, dip angle, and deep aquifer permeability (to
some extent). Notably, the parameters that are judged to be less influential on leakage rates are those
pertaining to the shallow aquifer (shallow aquifer permeability and porosity), well location, deep aquifer
101
40 60 80 100 120 140
Dip Angle (
o
)
0.0
0.2
0.4
0.6
0.8
1.0
Density
Complete CO 2 Runs
Incomplete Runs
0 2 4 6 8 10 12
Injection Rate (kg/s)
0.0
0.2
0.4
0.6
0.8
1.0
Density
Complete CO 2 Runs
Incomplete Runs
0 1 2
Well Location
0.0
0.2
0.4
0.6
0.8
1.0
0.05 0.10 0.15 0.20 0.25 0.30 0.35
Deep Aquifer Porosity
0.0
0.2
0.4
0.6
0.8
1.0
0.1 0.2 0.3 0.4 0.5
Shallow Aquifer Porosity
0.0
0.2
0.4
0.6
0.8
1.0
Density
14 13 12
Deep Aquifer Perm. (log 10 m
2
)
0.0
0.2
0.4
0.6
0.8
1.0
14 13 12
Shallow Aquifer Perm. (log 10 m
2
)
0.0
0.2
0.4
0.6
0.8
1.0
Density
0.00 0.02 0.04 0.06 0.08 0.10
Damage Zone Porosity
0.0
0.2
0.4
0.6
0.8
1.0
10 15 20 25 30 35 40 45 50
Injection Duration (yr)
0.0
0.2
0.4
0.6
0.8
1.0
Density
Complete CO 2 Runs
Incomplete Runs
10 15 20 25 30 35 40 45
Geothermal Gradient (
o
C/km)
0.0
0.2
0.4
0.6
0.8
1.0
15 14 13 12
Damage Zone Perm. (log 10 m
2
)
0.0
0.2
0.4
0.6
0.8
1.0
40 60 80 100 120 140
Dip Angle (
o
)
0.0
0.2
0.4
0.6
0.8
1.0
0 2 4 6 8 10 12
Injection Rate (kg/s)
0.0
0.2
0.4
0.6
0.8
1.0
Density
Figure 4.22: The distributional differences between the complete and incomplete simulations’ parameters
using the test set of CO
2
and the set of incomplete runs. Note that due to the use of a random selection
of the CO
2
simulations’ test set, these parameter distributions hold for the entire completed simulations’
dataset.
102
porosity, damage zone porosity, and geothermal gradient. These suggest that the ROM considers the leak-
age problem as a “balloon with a straw” problem, where pressure is propagated quickly through the reser-
voir and then relieved gradually through fault leakage, such that the fault is the primary constraint on
leakage. Notably, the FEHM simulator considers all parameters important, so we cannot take the param-
eters considered important by the ROM as the only ones that are important. However, given the good
performance of the ROM, we can seek to approximate the leakage ROM as consisting of two parts: a stor-
age aquifer modeled as a storage tank of fluid; and a fault modeled as a homogeneous, isotropic porous
medium. The implication for storage monitoring is that characterizing the deep aquifer and adjacent faults
may be sufficient to constrain the leakage rate without needing to identify specific leakage destinations.
This kind of targeted characterization and monitoring could reduce the cost of de-risking CO
2
storage sites
and assuring CO
2
containment, thereby moving GCS toward deployment.
With respect to the ROM’s performance over time series, we note that the preservation of time-series
statistics improves the ROM’s applicability to field cases. Figures 4.12, 4.13, 4.14, and 4.15 show the ROM’s
performance in preserving statistics over the entire duration of the simulation, which the ROM is able to do
despite not being trained to predict the entire time series but rather individual points in time. Time-series
statistics are important to preserve because analyses of uncertainty can often resort to summary statistics
to quantify leakage risk e.g., maximum leakage rates (D. R. Harp et al., 2016). By preserving these statistics,
our ROM preserves key elements of the distributions of outputs. This is especially notable in the case of
brine backflow (Figure 4.15): the preservation of brine backflow rates allows potential salinity changes in
the deep aquifer due to an influx of fluid from the shallow aquifer to be better estimated and summarized. If
the ROM had a high R
2
but poor time-series statistic preservation, we would have concluded that the ROM
does not capture the key dynamics of leakage and so would have not used it for leakage risk assessment.
103
4.4.3 ROMLimitationsandApplicability
There are limitations to this ROM’s applicability. The ROM understates the leakage values for cases
with “fast physics”, as described above. The geologic representation of the fault is a homogeneous, isotropic
medium. This is a simplified representation (Scholz, 2018), especially given the complexities of fault seal
analysis (Dockrill & Shipton, 2010; Mulrooney et al., 2020; Yielding et al., 1997). The fault is only three
meters thick; however, the variation in damage zone permeability results in varying fault transmissibilities,
which to an extent accounts for faults of differing thickness. The ROM also does not take into account
geomechanics or phase change, which can be important (Meguerdijian et al., 2022). The applicability of
the ROM is limited to deep saline aquifers in the supercritical CO
2
regime and does not extend to cases
with the presence of substantial amounts of hydrocarbon. The deep and shallow aquifer sizes are fixed.
The minimum sampled permeability is 10
− 15
m
2
(about 1 md), so the fault is never fully impermeable in
the simulations. This is a consequence of considering the fault as an upscaled fractured medium. The
range of values of fault parameters must be further clarified, and doing so is an area of future work. Fault
anisotropy is not taken into consideration though fault anisotropy is possible. Future work is needed to
understand the degree to which fault anisotropy can potentially cause changes in fault leakage behavior.
4.5 Conclusion
In this work, we build a fault-leakage ROM with good generalization for brine and CO
2
leakage rates.
The ROM uses two deep neural networks to predict brine and CO
2
leakage rates, respectively. The ROM
predicts individual time points as opposed to whole time series and contains no memory of previous pre-
dictions, unlike recurrent neural networks. The following achievements are noted:
1. We achieve a large speedup by using the ROM instead of the simulator, and time-series statistics are
preserved.
104
2. Based on sensitivity analyses, we find that the simulator considers all inputs significant while the
ROM considers the following inputs significant: damage zone permeability, simulation time, deep
aquifer permeability, injection rate, injection duration, dip angle.
3. The ROM exhibits good generalizability over the set of complete simulations and the available por-
tion of incomplete simulations.
4. Our analysis further includes an assessment of the ROM’s performance over the set of missing data,
the first time this is done for a fault-leakage ROM to the best of our knowledge.
We note that the ROM is expected to perform robustly over the missing data, which tends to exhibit
fast fluid flow physics, in the case of brine leakage but not for CO
2
leakage. We note that the inputs
considered significant by the ROM do not involve properties of the shallow aquifer, and this suggests a
new ROM conceptual model of leakage consisting only of the injection aquifer and a fault. This simplifies
the conceptual model by eliminating the shallow aquifer. However, we note that this conceptual model
would only apply to the prediction of brine and CO
2
leakage rates and not to the prediction of other risk
metrics.
Future work can expand the ROM’s validity and test the suggested conceptual model. Fault anisotropy
is not considered in this study but could be investigated to determine the impact of having horizon-
tal/vertical anisotropy and having a fault core within a damage zone where their properties differ. Other
risk metrics should be investigated, such as the CO
2
plume size in the shallow aquifer. This could yield
results showing that shallow aquifer properties are important for some risk metrics but not others. The
ROM can be expanded to additional cases involving different values for fault properties or aquifer sizes.
However, these values can be highly uncertain and need to be informed from field studies regarding faults
(e.g. fault apertures, extent of fracturing in a surrounding damage zone, fault seal analysis). Additional
105
ROMs can be generated to model the effect of fault leakage (as opposed to wellbore leakage, which has
been previously studied) on the aquifers into which leakage occurs.
106
Chapter5
ConclusionsandFutureWork
This dissertation has sought to address the knowledge gaps outlined in the Introduction. The first
study in Chapter 2 addressed the gap regarding the effect of fault dip angle on fault stability and leakage
dynamics, finding that the fault dip angle significantly impacts leakage, potentially even causing down-
ward leakage due to the combination of fault permeability increase due to slip and the buoyancy of the
leaking gas. Surface deformation effects are apparent from the fault dip angle, and this phenomenon can
be used for monitoring a GCS site. A new measure of leakage, the leakage magnitude, was also proposed
to complement the leakage mass. Future work in this area consists of combining fracturing and failure
mechanisms with fault leakage and dip angle effects in a combined study to identify potentially novel
phenomena arising from their combination, similar to the downward gas leakage in this study. There is
also an opportunity to study the use of surface deformation effects for inverting the fault dip angle (e.g.
using Bayesian inversion) during a GCS operation. This is currently the subject of an ongoing study by
the author and holds potential for improving site characterization during GCS operations.
The second study in Chapter 3 investigated the effects of thermal stresses and CO
2
solubility on CO
2
fault leakage. The study found that solubility effects are important for estimating CO
2
leakage rates and
thermal stresses are significant for estimating fault destabilization when CO
2
phase change is present.
107
The study further found that the spatial distribution of thermal stresses is different from that of poroelas-
tic stresses and that thermal stresses are not necessarily overwhelmed by poroelastic stresses. The future
work for this study was the development of fault-leakage ROMs, which was pursued in the third study
in Chapter 4. This third study involved running an ensemble of simulations and then building a ROM
using DNNs. The ROM’s performance and sensitivities were analyzed. The ROM’s performance was also
assessed over the parameter space of failed simulations, a novel addition that can help build confidence
among operators for this and future ROMs. Based on the ROM’s sensitivity analyses, a conceptual model
for leakage rates consisting of a reservoir and fault was suggested. Future work for this study involves
investigating risk metrics other than leakage rates (e.g. plume sizes), taking into account the presence of
a fault consisting of both a damage zone and fault core with anisotropic properties, varying the properties
of the fault along the fault (as a consequence of fault seal analysis), and the generation of additional ROMs.
Among these ROMs can be: fault-leakage ROMs that also consider geomechanics and CO
2
solubility ef-
fects, fault-leakage aquifer ROMs to estimate the impact of fault leakage on aquifer properties (e.g. pH,
Total Dissolved Solids, metal ion concentration changes), and ROMs with expanded parameter spaces that
include different aquifer sizes and longer distances from the well to the fault.
On the whole, data integration and derisking emerge as enabling areas of future work. The study
of the fault-leakage dynamics of CO
2
and its risk assessment would benefit from having increased data
available. For example, the incorporation of field data into ROM development could help constrain the
values of varying parameters. Efforts toward this end have been pursued by the DOE via the Energy Data
eXchange (EDX) initiative (National Energy Technology Laboratory, n.d.), but further efforts are needed
which integrate academia, national laboratories, industry, and other stakeholders. The data which exists
should be made easily available, and data sharing and curation should be promoted. Finally, fault-leakage
risk assessment should be incorporated into derisking efforts. Large amounts of pore space are anticipated
to be needed for GCS to affect climate change mitigation at scale, and faults can have different leakage
108
patterns than other leakage sources, such as wellbores. Efforts to advance the certification of pore space
at scale can generate data regarding fault-leakage dynamics, which can be analyzed and used to improve
the risk assessment of future GCS sites.
109
Bibliography
26 USC 45Q: Credit for carbon oxide sequestration. (2019). Retrieved October 20, 2019, from
https://uscode.house.gov/view.xhtml?req=(title:26%5C%20section:45Q%5C%20edition:prelim)
Alexander, W. R., Reijonen, H. M., & McKinley, I. G. (2015). Natural analogues: studies of geological
processes relevant to radioactive waste disposal in deep geological repositories. Swiss J. Geosci.,
108, 75–100.
Amooie, M. A., Soltanian, M. R., & Moortgat, J. (2018). Solutal convection in porous media: Comparison
between boundary conditions of constant concentration and constant flux. Physical Review E,
98(3). https://doi.org/10.1103/PhysRevE.98.033118
Bach, W. (1979). Impact of increasing atmospheric CO
2
concentrations on global climate: Potential
consequences and corrective measures. Environment International, 2(4), 215–228.
https://doi.org/10.1016/0160-4120(79)90004-7
Bachu, S. (2000). Sequestration of CO
2
in geological media: Criteria and approach for site selection in
response to climate change. Energy Conversion Management, 41, 953–970.
https://doi.org/10.1016/S0196-8904(99)00149-1
Bianchi, M., Zheng, L., & Birkholzer, J. T. (2016). Combining multiple lower-fidelity models for emulating
complex model responses for CCS environmental risk assessment. International Journal of
Greenhouse Gas Control, 46, 248–258. https://doi.org/10.1016/j.ijggc.2016.01.009
Biot, M. A., & Willis, D. G. (1957). The Elastic Coefficients of the Theory of Consolidation. Journal of
Applied Mechanics, 594–601. http://www.pmi.ou.edu/Biot2005/papers/FILES/068.PDF
Birkholzer, J. T., Oldenburg, C. M., & Zhou, Q. (2015). CO
2
migration and pressure evolution in deep
saline aquifers. International Journal of Greenhouse Gas Control, 40, 203–220.
Birkholzer, J. T., & Zhou, Q. (2009). Basin–scale hydrogeologic impacts of CO
2
storage: Capacity and
regulatory implications. International Journal of Greenhouse Gas Control, 3, 745–756.
Borgonovo, E. (2007). A new uncertainty importance measure. Reliability Engineering & System Safety,
92(6), 771–784. https://doi.org/10.1016/j.ress.2006.04.015
110
Brooks, R. H., & Corey, A. T. (1964). Hydraulic properties of porous media. Hydrological Papers (Colorado
State University), 3, 892–898.
Burghardt, J., & Appriou, D. (2021). State of Stress Uncertainty Quantification and Geomechanical Risk
Analysis for Subsurface Engineering. Retrieved February 22, 2022, from
http://onepetro.org/ARMAUSRMS/proceedings/ARMA21/All-ARMA21/ARMA-2021-2129/468335
Buuren, S. v. (2018). Flexible Imputation of Missing Data (2nd ed.) [Also available online in HTML format
at: https://stefvanbuuren.name/fimd/]. Chapman; Hall/CRC.
https://doi.org/10.1201/9780429492259
Caine, J. S., & Forster, C. B. (1999). Fault Zone Architecture and Fluid Flow: Insights from Field Data and
Numerical Modeling. Faults and Subsurface Fluid Flow in the Shallow Crust (pp. 101–127).
American Geophysical Union (AGU). https://doi.org/10.1029/GM113p0101
Cappa, F., & Rutqvist, J. (2011). Impact of CO
2
geological sequestration on the nucleation of earthquakes.
Geophys. Res. Lett., 38, L17313.
CC. (2019). Creative Commons – Attribution 4.0 International – CC BY 4.0 [Accessed: 2019-05-16].
Cesca, S., Grigoli, F., Heimann, S., Gonzalez, A., Buforn, E., Maghsoudi, S., Blanch, E., & Dahm, T. (2014).
The 2013 September–October seismic sequence offshore Spain: a case of seismicity triggered by
gas injection? Geophys. J. Int., 198, 941–953. https://doi.org/10.1093/gji/ggu172
Chaves, E. J., Schwartz, S. Y., & Abercrombie, R. E. (2020). Repeating earthquakes record fault weakening
and healing in areas of megathrust postseismic slip. Sci. Adv., 6.
https://doi.org/10.1126/sciadv.aaz9317
Chester, F. M., Evans, J. P., & Biegel, R. L. (1993). Internal structure and weakening mechanisms of the San
Andreas Fault. Journal of Geophysical Research: Solid Earth, 98(B1), 771–786.
https://doi.org/10.1029/92JB01866
Coussy, O., Eymard, R., & Lassabatére, T. (1998). Constitutive modeling of unsaturated drying deformable
materials. J. Eng. Mech., 124(6), 658–557.
Coussy, O. (1995). Mechanics of porous continua. John Wiley; Sons.
Coussy, O. (2004). Poromechanics. John Wiley & Sons.
The CSIRO In-Situ Laboratory. (n.d.). Retrieved June 7, 2022, from https://research.csiro.au/in-situ/
Dana, S., Zhao, X., & Jha, B. (2022). A two-grid simulation framework for fast monitoring of fault stability
and ground deformation in multiphase geomechanics. J. Comput. Phys., 466, 111405.
https://doi.org/10.1016/j.jcp.2022.111405
Dempsey, D., Kelkar, S., Lewis, K., Hickman, S., Davatzes, N., Moos, D., & Zemach, E. (2013). Modeling
shear stimulation of the desert peak EGS well 27-15 using a coupled
thermal-hydrological-mechanical simulator. 47th US Rock Mechanics/Geomechanics Symposium,
(ARMA-2013-608).
111
Dockrill, B., & Shipton, Z. K. (2010). Structural controls on leakage from a natural CO
2
geologic storage
site: Central Utah, U.S.A. Journal of Structural Geology, 32(11), 1768–1782.
https://doi.org/10.1016/j.jsg.2010.01.007
Dor, O., Ben-Zion, Y., Rockwell, T. K., & Brune, J. (2006). Pulverized rocks in the Mojave section of the
San Andreas Fault Zone. Earth and Planetary Science Letters, 245(3), 642–654.
https://doi.org/10.1016/j.epsl.2006.03.034
Downey, M. W. (1984). Evaluating seals for hydrocarbon accumulations. Am. Assoc. Petrol. Geol., 68,
1752–1763.
Duan, Q., Chen, J., & Yang, X. (2020). A comparison of gas and water permeability in clay-bearing fault
and reservoir rocks: Experimental results and evolution mechanisms. J. Geophys. Res., 125.
https://doi.org/https://doi.org/10.1029/2019JB018278
Ershadnia, R., Wallace, C. D., & Soltanian, M. R. (2020). CO
2
geological sequestration in heterogeneous
binary media: Effects of geological and operational conditions. Advances in Geo-Energy Research,
4(4), 392–405. https://doi.org/10.46690/ager.2020.04.05
Fialko, Y., & Simons, M. (2000). Deformation and seismicity in the Coso geothermal area, Inyo County,
California: Observations and modeling using satellite radar interferometry. J. Geophys. Res., 105,
781–794.
Gable, C. W., Gammel, T., George, D., Kuprat, A., Miller, T. A., Ortega, F. A., Trease, H., Trease, L., &
Walker, R. B. (2021). Los Alamos Grid Toolbox (LaGriT). Retrieved March 12, 2022, from
https://lanl.github.io/LaGriT/
Ghanbarzadeh, S., Hesse, M. A., Prodanovic, M., & Gardner, J. E. (2015). Deformation-assisted fluid
percolation in rock salt. Science, 350, 1069–1072.
Goebel, T., Hosseini, S., Cappa, F., Hauksson, E., Ampuero, J.-P., Aminzadeh, F., & Saleeby, J. (2016).
Wastewater disposal and earthquake swarm activity at the southern end of the Central Valley,
California. Geophys. Res. Lett., 45, 1092–1099. https://doi.org/10.1002/2015GL066948
Guglielmi, Y., Cappa, F., Avouac, J.-P., Henry, P., & Elsworth, D. (2015). Seismicity triggered by fluid
injection–induced aseismic slip. Science, 348, 1224–1226.
Guglielmi, Y., Birkholzer, J., Rutqvist, J., Jeanne, P., & Nussbaum, C. (2017). Can Fault Leakage Occur
Before or Without Reactivation? Results from an in Situ Fault Reactivation Experiment at Mont
Terri. Energy Procedia, 114, 3167–3174. https://doi.org/10.1016/j.egypro.2017.03.1445
Guglielmi, Y., Nussbaum, C., Cappa, F., De Barros, L., Rutqvist, J., & Birkholzer, J. (2021). Field-scale fault
reactivation experiments by fluid injection highlight aseismic leakage in caprock analogs:
Implications for CO
2
sequestration. International Journal of Greenhouse Gas Control, 111, 103471.
https://doi.org/10.1016/j.ijggc.2021.103471
Hanks, T. C., & Kanamori, H. (1979). A moment magnitude scale. J. Geophys. Res., 84, 2348–2350.
112
Hansen, O., Gilding, D., Nazarian, B., Osdal, B., Ringrose, P., & Kristoffersen, J.-B. (2013). Snøhvit: The
history of injecting and storing 1 Mt CO
2
in the fluvial Tubaen Fm. Energy Procedia, 37,
3565–3573.
Harbert, W., Daley, T. M., Bromhal, G., Sullivan, C., & Huang, L. (2016). Progress in monitoring strategies
for risk reduction in geologic CO
2
storage. International Journal of Greenhouse Gas Control, 51,
260–275.
Harp, D., Vasylkivska, V., Porter, M., Chen, B., Dempsey, D., & King, S. (2013). Model Analysis ToolKit
(MATK) - Python toolkit for model analysis — MATK 0 documentation. Retrieved March 13,
2022, from http://dharp.github.io/matk/
Harp, D. R., Pawar, R., Carey, J. W., & Gable, C. W. (2016). Reduced order models of transient CO
2
and
brine leakage along abandoned wellbores from geologic carbon sequestration reservoirs.
International Journal of Greenhouse Gas Control, 45, 150–162.
https://doi.org/10.1016/j.ijggc.2015.12.001
Helton, J. C., & Davis, F. J. (2002). Latin Hypercube Sampling and the Propagation of Uncertainty in
Analyses of Complex Systems (Sandia Report SAND2001-0417). Sandia National Laboratories.
United States of America. Retrieved April 24, 2021, from https://core.ac.uk/reader/194436110
Hepburn, C., Adlen, E., Beddington, J., Carter, E. A., Fuss, S., Mac Dowell, N., Minx, J. C., Smith, P., &
Williams, C. K. (2019). The technological and economic prospects for CO
2
utilization and
removal. Nature, 575(7781), 87–97. https://doi.org/10.1038/s41586-019-1681-6
Hepple, R., & Benson, S. M. (2002). Implications of surface seepage on the effectiveness of geologic
storage of carbon dioxide as a climate change mitigation strategy. In J. Gale & Y. Kaya (Eds.),
Greenhouse Gas Control Technologies (pp. 261–266). Elsevier.
Hepple, R. P., & Benson, S. M. (2005). Geologic storage of carbon dioxide as a climate change mitigation
strategy: Performance requirements and the implications of surface seepage. Environmental
Geology, 47(4), 576–585. https://doi.org/10.1007/s00254-004-1181-2
Herman, J., & Usher, W. (2017). SALib: An open-source python library for sensitivity analysis. The
Journal of Open Source Software, 2(9). https://doi.org/10.21105/joss.00097
Hochreiter, S., & Schmidhuber, J. (1997). Long Short-Term Memory. Neural Computation, 9(8), 1735–1780.
https://doi.org/10.1162/neco.1997.9.8.1735
Homma, T., & Saltelli, A. (1996). Importance measures in global sensitivity analysis of nonlinear models.
Reliability Engineering & System Safety, 52(1), 1–17.
https://doi.org/10.1016/0951-8320(96)00002-6
Horton, S. (2012). Disposal of hydrofracking waste fluid by injection into subsurface aquifers triggers
earthquake swarm in central Arkansas with potential for damaging earthquake.Seismo.Res.Lett.,
83, 250–260.
113
Hosseini, S. M., Goebel, T. H. W., Jha, B., & Aminzadeh, F. (2018). A probabilistic approach to
injection-induced seismicity assessment in the presence and absence of flow boundaries.
Geophys. Res. Lett., 45, 8182–8189. https://doi.org/10.1029/2018GL077552
Ide, S., & Takeo, M. (1997). Determination of constitutive relations of fault slip based on seismic wave
analysis. J. Geophys. Res., 102, 27379–27391.
IPCC. (2005). IPCC special report on carbon dioxide capture and storage (B. Metz, O. Davidson,
H. De Coninck, M. Loos, & L. Meyer, Eds.; tech. rep.). Intergovernmental Panel on Climate
Change. Cambridge University Press, UK.
Jacquey, A. B., Cacace, M., Blöcher, G., & Scheck-Wenderoth, M. (2015). Numerical Investigation of
Thermoelastic Effects on Fault Slip Tendency during Injection and Production of Geothermal
Fluids. Energy Procedia, 76, 311–320. https://doi.org/10.1016/j.egypro.2015.07.868
Jaeger, J. C., & Cook, N. G. W. (1979). Fundamentals of rock mechanics. Chapman; Hall.
Jagalur-Mohan, J., Jha, B., Wang, Z., Juanes, R., & Marzouk, Y. (2018). Inferring fault frictional and
reservoir hydraulic properties from injection-induced seismicity. Geophys. Res. Lett., 45,
1313–1320. https://doi.org/10.1002/2017GL075925
Jha, B., Bottazzi, F., Wojcik, R., Coccia, M., Bechor, N., McLaughlin, D., Herring, T., Hager, B. H.,
Mantica, S., & Juanes, R. (2015). Reservoir characterization in an underground gas storage field
using joint inversion of flow and geodetic data. Int. J. Numer. Anal. Meth. Geomech., 39,
1619–1638. https://doi.org/10.1002/nag.2427
Jha, B., & Juanes, R. (2014). Coupled multiphase flow and poromechanics: A computational model of pore
pressure effects on fault slip and earthquake triggering. Wat. Resour. Res., 50(5), 3776–3808.
https://doi.org/10.1002/2013WR015175
Johri, M. (2012). Fault Damage Zones - Observations, Dynamic Modeling, and Implications of Fluid Flow
(Doctoral dissertation). Stanford University. Stanford, California. Retrieved June 29, 2020, from
https://pangea.stanford.edu/departments/geophysics/dropbox/SRB/public/docs/theses/SRB_131_
OCT12_Johri.pdf
Jones, D. G., Barkwith, A. K. A. P., Hannis, S., Lister, T. R., Gal, F., Graziani, S., Beaubien, S. E., &
Widory, D. (2014). Monitoring of near surface gas seepage from a shallow injection experiment at
the CO
2
Field Lab, Norway. International Journal of Greenhouse Gas Control, 28, 300–317.
https://doi.org/10.1016/j.ijggc.2014.06.021
Jung, N.-H., Han, W. S., Han, K., & Park, E. (2015). Regional-scale advective, diffusive, and eruptive
dynamics of CO
2
and brine leakage through faults and wellbores. Journal of Geophysical
Research: Solid Earth, 120(5), 3003–3025. https://doi.org/10.1002/2014JB011722
Jung, Y., Zhou, Q., & Birkholzer, J. T. (2013). Early detection of brine and CO
2
leakage through abandoned
wells using pressure and surface-deformation monitoring data: Concept and demonstration. Adv.
Wat. Resour., 62, 555–569.
114
Karimi-Fard, M., Durlofsky, L., & Aziz, K. (2004). An Efficient Discrete-Fracture Model Applicable for
General-Purpose Reservoir Simulators. SPE Journal, 9, 227–236.
https://doi.org/10.2118/88812-PA
Karimi-Fard, M., & Firoozabadi, A. (2003). Numerical Simulation of Water Injection in Fractured Media
Using the Discrete-Fracture Model and the Galerkin Method. SPE Reservoir Evaluation &
Engineering, 6, 117–126. https://doi.org/10.2118/83633-PA
Keating, E. H., Fessenden, J., Kanjorski, N., Koning, D. J., & Pawar, R. (2010). The impact of CO
2
on
shallow groundwater chemistry: Observations at a natural analog site and implications for
carbon sequestration. Environmental Earth Sciences, 60(3), 521–536.
https://doi.org/10.1007/s12665-009-0192-4
Keating, E. H., Harp, D. H., Dai, Z., & Pawar, R. J. (2016). Reduced order models for assessing CO
2
impacts
in shallow unconfined aquifers. International Journal of Greenhouse Gas Control, 46, 187–196.
https://doi.org/10.1016/j.ijggc.2016.01.008
Keranen, K. M., Savage, H. M., Abers, G. A., & Cochran, E. S. (2013). Potentially induced earthquakes in
Oklahoma, USA: Links between wastewater injection and the 2011 M
w
5.7 earthquake sequence.
Geology, 41, 699–702.
Khilyuk, L. F., Chilingar, G. V., Robertson, J. O., Jr., & Endres, B. (2000). Gas Migration : Events Preceding
Earthquakes. Elsevier. https://doi.org/10.1016/B978-0-88415-430-3.X5000-9
Kingma, D. P., & Ba, J. (2017). Adam: A Method for Stochastic Optimization [arXiv: 1412.6980].
arXiv:1412.6980 [cs]. https://doi.org/https://doi.org/10.48550/arXiv.1412.6980
Klusman, R. W. (2003). Evaluation of leakage potential from a carbon dioxide EOR/sequestration project.
Energy Conversion Management, 44, 1921–1940.
Kolmogorov, A. (1933). Sulla determinazione empirica di una lgge di distribuzione. Inst. Ital. Attuari,
Giorn., 4, 83–91.
Lee, H. S., & Cho, T. F. (2002). Hydraulic Characteristics of Rough Fractures in Linear Flow under Normal
and Shear Load. Rock Mechanics and Rock Engineering, 35(4), 299–318.
https://doi.org/10.1007/s00603-002-0028-y
Lewicki, J. L., Evans, W. C., Hilley, G. E., Sorey, M. L., Rogie, J. D., & Brantley, S. L. (2003). Shallow soil
CO
2
flow along the San Andreas and Calaveras faults, CA. J. Geophys. Res., 108, 2187.
https://doi.org/10.1029/2002JB002141
Lewicki, J. L., Hilley, G. E., & Oldenburg, C. M. (2005). An improved strategy to detect CO
2
leakage for
verification of geologic carbon sequestration. Geophys. Res. Lett., 32, L19403.
Lewis, R. W., & Schrefler, B. A. (1998). The Finite Element Method in the Static and Dynamic deformation
and Consolidation of Porous Media (Second Edition). Wiley.
Lewis, R. W., & Sukirman, Y. (1993). Finite element modelling of three-phase flow in deforming saturated
oil reservoirs. Int. J. Numer. Anal. Methods Geomech., 17, 577–598.
115
Lindner, E. (2020). NRAP - Wednesday Morning Part 4 : NRAP fault
modeling: A new ROM for fault leakage. Retrieved February 25, 2021, from
https://www.youtube.com/watch?v=5QmfXMqDIUc&feature=youtu.be
Llewellyn, G. T., Dorman, F., Westland, J. L., Yoxtheimer, D., Grieve, P., Sowers, T., Humston-Fulmer, E., &
Brantley, S. L. (2015). Evaluating a groundwater supply contamination incident attributed to
Marcellus Shale gas development. Proc. Natl. Acad. Sci. U.S.A., 112, 6325–6330.
https://doi.org/10.1073/pnas.1420279112
Lowell, R., Kolandaivelu, K., & Rona, P. (2014). Hydrothermal activity. Reference module in earth systems
and environmental sciences. Elsevier.
https://doi.org/https://doi.org/10.1016/B978-0-12-409548-9.09132-6
Lu, C., Sun, Y., Buscheck, T. A., Hao, Y., White, J. A., & Chiaramonte, L. (2012). Uncertainty quantification
of CO
2
leakage through a fault with multiphase and nonisothermal effects. Greenhouse Gases:
Science and Technology, 2(6), 445–459. https://doi.org/https://doi.org/10.1002/ghg.1309
Luth, S., Ivanova, A., Ivandic, M., & Gotz, J. (2015). 4d seismic monitoring at the Ketzin pilot site during
five years of storage–results and quantitative assessment. Energy Procedia, 76, 536–542.
Ma, X., & Elbanna, A. E. (2015). Effect of off-fault low-velocity elastic inclusions on supershear rupture
dynamics. Geophys. J. Int., 203, 664–677. https://doi.org/10.1093/gji/ggv302
Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro,
Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow,
Andrew Harp, Geoffrey Irving, Michael Isard, Jia, Y., Rafal Jozefowicz, Lukasz Kaiser,
Manjunath Kudlur, ... Xiaoqiang Zheng. (2015). TensorFlow: Large-scale machine learning on
heterogeneous systems [Software available from tensorflow.org]. https://www.tensorflow.org/
Mazzoldi, A., Rinaldi, A. P., Borgia, A., & Rutqvist, J. (2012). Induced seismicity within geological carbon
sequestration projects: Maximum earthquake magnitude and leakage potential from undetected
faults. International Journal of Greenhouse Gas Control, 10, 434–442.
https://doi.org/10.1016/j.ijggc.2012.07.012
McKay, M. D., Conover, W. J., & Beckman, R. J. (1979). A comparison of three methods for selection values
of input variables in the analysis of output from a computer code. Technometrics, 22(2), 239–245.
Meguerdijian, S., & Jha, B. (2021). Quantification of fault leakage dynamics based on leakage magnitude
and dip angle. Int. J. Numer. Anal. Methods Geomech., 45(16), 2303–2320.
https://doi.org/10.1002/nag.3267
Meguerdijian, S., Pawar, R. J., Harp, D. R., & Jha, B. (2022). Thermal and solubility effects on fault leakage
during geologic carbon storage. International Journal of Greenhouse Gas Control, 116, 103633.
https://doi.org/https://doi.org/10.1016/j.ijggc.2022.103633
116
Middleton, R. S., Chen, B., Harp, D. R., Kammer, R. M., Ogland-Hand, J. D., Bielicki, J. M., Clarens, A. F.,
Currier, R. P., Ellett, K. M., Hoover, B. A., McFarlane, D. N., Pawar, R. J., Stauffer, P. H.,
Viswanathan, H. S., & Yaw, S. P. (2020). Great SCO
2
T! Rapid tool for carbon sequestration
science, engineering, and economics. Applied Computing and Geosciences, 7, 100035.
https://doi.org/10.1016/j.acags.2020.100035
Mikumo, T., Olsen, K. B., Fukuyama, E., & Yagi, Y. (2003). Stress-breakdown time and slip-weakening
distance inferred from slip-velocity functions on earthquake faults. Bull. Seismol. Soc. Amer., 93,
264–282.
Miocic, J. M., Gilfillan, S. M., Frank, N., Schroeder-Ritzrau, A., Burnside, N. M., & Haszeldine, R. S. (2019).
420,000 year assessment of fault leakage rates shows geological carbon storage is secure.
Scientific Reports , 9(1), 1–9. https://doi.org/10.1038/s41598-018-36974-0
Miocic, J. M., Johnson, G., & Gilfillan, S. M. V. (2020). Stress field orientation controls on fault leakage at a
natural CO
2
reservoir. Solid Earth, 11(4), 1361–1374. https://doi.org/10.5194/se-11-1361-2020
Mitchell, T. M., Ben-Zion, Y., & Shimamoto, T. (2011). Pulverized fault rocks and damage asymmetry
along the Arima-Takatsuki Tectonic Line, Japan. Earth and Planetary Science Letters, 308(3),
284–297. https://doi.org/10.1016/j.epsl.2011.04.023
Morris, J. P., Hao, Y., Foxall, W., & McNab, W. (2011). A study of injection-induced mechanical
deformation at the In Salah CO
2
storage project. International Journal of Greenhouse Gas Control,
5(7), 270–280.
Mulrooney, M. J., Osmond, J. L., Skurtveit, E., Faleide, J. I., & Braathen, A. (2020). Structural analysis of
the Smeaheia fault block, a potential CO
2
storage site, northern Horda Platform, North Sea.
Marine and Petroleum Geology, 121, 104598. https://doi.org/10.1016/j.marpetgeo.2020.104598
National Energy Technology Laboratory. (n.d.). EDX - NETL’s Energy Data eXchange. Retrieved August
6, 2022, from https://edx.netl.doe.gov/
Newell, P., & Ilgen, A. G. (Eds.). (2019). Science of Carbon Storage in Deep Saline Formations. Elsevier.
https://doi.org/10.1016/C2016-0-03237-0
Nicot, J.-P., Oldenburg, C. M., Houseworth, J. E., & Choi, J.-W. (2013). Analysis of potential leakage
pathways at the Cranfield MS, USA, CO
2
sequestration site. International Journal of Greenhouse
Gas Control, 18, 388–400. https://doi.org/10.1016/j.ijggc.2012.10.011
Okazaki, K., Katayama, I., & Noda, H. (2013). Shear-induced permeability anisotropy of simulated
serpentinite gouge produced by triaxial deformation experiments. Geophys. Res. Lett., 40,
1290–1294. https://doi.org/https://doi.org/10.1002/grl.50302
Olsen, K. B., Madariaga, R., & Archuleta, R. J. (1997). Three-dimensional dynamic simulation of the 1992
Landers Earthquake. Science, 278, 834–838.
117
Onishi, T., Nguyen, M. C., Carey, J. W., Will, B., Zaluski, W., Bowen, D. W., Devault, B. C., Duguid, A.,
Zhou, Q., Fairweather, S. H., Spangler, L. H., & Stauffer, P. H. (2019). Potential CO
2
and brine
leakage through wellbore pathways for geologic CO
2
sequestration using the National Risk
Assessment Partnership tools: Application to the Big Sky Regional Partnership. International
Journal of Greenhouse Gas Control, 81, 44–65. https://doi.org/10.1016/j.ijggc.2018.12.002
Osborn, S. G., Vengosh, A., Warner, N. R., & Jackson, R. B. (2011). Methane contamination of drinking
water accompanying gas-well drilling and hydraulic fracturing. Proc. Natl. Acad. Sci. U.S.A.,
108(20), 8172–8176.
Panda, D., Kundu, B., Gahalaut, V. K., Burgmann, R., Jha, B., Asaithambi, R., Yadav, R. K., Vissa, N. K., &
Bansal, A. K. (2018). Seasonal modulation of deep slow-slip and earthquakes on the Main
Himalayan Thrust. Nature Comm., 9(4140). https://doi.org/10.1038/s41467-018-06371-2
Pawar, R. J., Bromhal, G. S., Chu, S., Dilmore, R. M., Oldenburg, C. M., Stauffer, P. H., Zhang, Y., &
Guthrie, G. D. (2016). The National Risk Assessment Partnership’s integrated assessment model
for carbon storage: A tool to support decision making amidst uncertainty. International Journal
of Greenhouse Gas Control, 52, 175–189. https://doi.org/10.1016/j.ijggc.2016.06.015
Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P.,
Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., &
Duchesnay, E. (2011). Scikit-learn: Machine learning in Python. Journal of Machine Learning
Research, 12, 2825–2830. https://jmlr.csail.mit.edu/papers/v12/pedregosa11a.html,https:
//scikit-learn.org/stable/modules/model_evaluation.html#r2-score
Pereira, L. C., Guimaraes, L. J. N., Horowitz, B., & Sanchez, M. (2014). Coupled hydro-mechanical fault
reactivation analysis incorporating evidence theory for uncertainty quantification. Comput.
Geotech., 56, 202–215.
Pianosi, F., & Wagener, T. (2015). A simple and efficient method for global sensitivity analysis based
on cumulative distribution functions. Environmental Modelling & Software, 67, 1–11.
https://doi.org/10.1016/j.envsoft.2015.01.004
Pianosi, F., & Wagener, T. (2018). Distribution-based sensitivity analysis from a generic input-output
sample. Environmental Modelling & Software, 108, 197–207.
https://doi.org/10.1016/j.envsoft.2018.07.019
Pickup, G., Kiatsakulphan, M., & Mills, J. (2010). Analysis of grid resolution for simulations of CO
2
storage in deep saline aquifers. ECMOR XII - 12th European Conference on the Mathematics of Oil
Recovery, Article cp-163-00021. https://doi.org/https://doi.org/10.3997/2214-4609.20144939
Plischke, E., Borgonovo, E., & Smith, C. L. (2013). Global sensitivity measures from given data. European
Journal of Operational Research, 226(3), 536–550. https://doi.org/10.1016/j.ejor.2012.11.047
Pruess, K. (2005). Numerical studies of fluid leakage from a geologic disposal reservoir for CO
2
show
self-limiting feedback between fluid flow and heat transfer. Geophys. Res. Lett., 32, L14404.
https://doi.org/10.1029/2005GL023250
118
Rahman, M. J., Choi, J. C., Fawad, M., & Mondol, N. H. (2021). Probabilistic analysis of Vette fault stability
in potential CO
2
storage site Smeaheia, offshore Norway. International Journal of Greenhouse Gas
Control, 108, 103315. https://doi.org/10.1016/j.ijggc.2021.103315
Razavi, S., Tolson, B. A., & Burn, D. H. (2012). Review of surrogate modeling in water resources. Wat.
Resour. Res., 48(7). https://doi.org/10.1029/2011WR011527
Reasenberg, P. A., & Simpson, R. W. (1992). Response of Regional Seismicity to the Static Stress Change
Produced by the Loma Prieta Earthquake. Science, 255(5052), 1687–1690.
https://doi.org/10.1126/science.255.5052.1687
Riding, J. B., & Rochelle, C. A. (2009). Subsurface characterisation and geological monitoring of the CO
2
injection operation at Weyburn, Saskatchewan, Canada. In D. J. Evans & R. A. Chadwick (Eds.),
Underground Gas Storage: Worldwide Experiences and Future Development in the UK and Europe:
London (pp. 227–256). Spec. Publ. Geol. Soc.
Rinaldi, A. P., Vilarrasa, V., Rutqvist, J., & Cappa, F. (2015). Fault reactivation during CO
2
sequestration:
Effects of well orientation on seismicity and leakage. Greenhouse Gases:Science and Technology,5,
645–656. https://doi.org/10.1002/ghg.1511
Rinaldi, A. P., & Rutqvist, J. (2013). Modeling of deep fracture zone opening and transient ground surface
uplift at KB-502 CO
2
injection well, In Salah, Algeria. International Journal of Greenhouse Gas
Control, 12, 155–167. https://doi.org/10.1016/j.ijggc.2012.10.017
Rinaldi, A. P., Rutqvist, J., & Cappa, F. (2014). Geomechanical effects on CO
2
leakage through fault zones
during large-scale underground injection. International Journal of Greenhouse Gas Control, 20,
117–131. https://doi.org/10.1016/j.ijggc.2013.11.001
Robertson, E. (1988). Thermal properties of rocks (USGS Numbered Series No. 88-441). U.S. Geological
Survey. Retrieved February 21, 2020, from http://pubs.er.usgs.gov/publication/ofr88441
Romanak, K. D., Bennett, P. C., Yang, C., & Hovorka, S. D. (2012). Process-based approach to CO
2
leakage
detection by vadose zone gas monitoring at geologic CO
2
storage sites. Geophys. Res. Lett., 39,
L15405. https://doi.org/10.1029/2012GL052426
Rucci, A., Vasco, D. W., & Novali, F. (2013). Monitoring the geologic storage of carbon dioxide using
multicomponent SAR interferometry. Geophys. J. Int., 193, 197–208.
https://doi.org/10.1093/gji/ggs112
Rutqvist, J., Birkholzer, J. T., & Tsang, C. F. (2008). Coupled reservoir–geomechanical analysis of the
potential for tensile and shear failure associated with CO
2
injection in multilayered
reservoir–caprock systems. Int. J. Rock Mech. Min. Sci., 45, 132–143.
Saltelli, A. (2002a). Making best use of model evaluations to compute sensitivity indices. Computer
Physics Communications, 145(2), 280–297. https://doi.org/10.1016/S0010-4655(02)00280-1
Saltelli, A. (2002b). Sensitivity Analysis for Importance Assessment. Risk Analysis, 22(3), 579–590.
https://doi.org/10.1111/0272-4332.00040
119
Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., & Tarantola, S. (2010). Variance based
sensitivity analysis of model output. Design and estimator for the total sensitivity index.
Computer Physics Communications, 181(2), 259–270. https://doi.org/10.1016/j.cpc.2009.09.018
Scholz, C. H. (1988). The critical slip distance for seismic faulting. Nature, 336, 761–763.
Scholz, C. H. (2018). The mechanics of earthquakes and faulting. Cambridge University Press. Retrieved
September 22, 2020, from https://doi.org/10.1017/9781316681473
Shakiba, M., & Hosseini, S. A. (2016). Detection and characterization of CO
2
leakage by multi-well pulse
testing and diffusivity tomography maps. International Journal of Greenhouse Gas Control, 54,
15–28. https://doi.org/10.1016/j.ijggc.2016.08.015
Sibson, R. H. (1977). Fault rocks and fault mechanisms. Journal of the Geological Society, 133(3), 191–213.
https://doi.org/10.1144/gsjgs.133.3.0191
Sibson, R. H. (1994). Crustal stress, faulting and fluid flow. Geol. Soc. London, Special Publications, 78,
69–84.
Singh, M., Chaudhuri, A., Soltanian, M. R., & Stauffer, P. H. (2021). Coupled multiphase flow and
transport simulation to model CO
2
dissolution and local capillary trapping in permeability and
capillary heterogeneous reservoir. International Journal of Greenhouse Gas Control, 108, 103329.
https://doi.org/10.1016/j.ijggc.2021.103329
Sinn, H.-W. (2012). The Green Paradox: A Supply-Side Approach to Global Warming. MIT Press. Retrieved
December 13, 2019, from
http://ebookcentral.proquest.com/lib/socal/detail.action?docID=3339373
Smirnov, N. V. (1939). On the estimation of the discrepancy between empirical curves of distribution for
two independent samples. Bull. Math. Univ. Moscou, 2(2), 3–14.
Snippe, J., Kampman, N., Bisdom, K., Tambach, T., March, R., Maier, C., Phillips, T., Inskip, N. F.,
Doster, F., & Busch, A. (2022). Modelling of long-term along-fault flow of CO
2
from a natural
reservoir. International Journal of Greenhouse Gas Control, 118, 103666.
https://doi.org/10.1016/j.ijggc.2022.103666
Sobol’, I. M. (2001). Global sensitivity indices for nonlinear mathematical models and their Monte Carlo
estimates. Mathematics and Computers in Simulation, 55(1-3), 271–280.
https://doi.org/10.1016/S0378-4754(00)00270-6
Stauffer, P., Chu, S., Tauxe, C., & Pawar, R. (2020). NRAP Integrated Assessment Model-Carbon Storage
(NRAP-IAM-CS) Tool Users Manual Version: 2016.11-1.1. https://doi.org/10.18141/1592692
Strazisar, B. R., Wells, A. W., Diehl, J. R., Hammack, R. W., & Veloski, G. A. (2009). Near-surface
monitoring for the ZERT shallow CO
2
injection project. International Journal of Greenhouse Gas
Control, 3(6), 736–744. https://doi.org/10.1016/j.ijggc.2009.07.005
120
Streit, J. E., & Hillis, R. R. (2004). Estimating fault stability and sustainable fluid pressures for
underground storage of CO
2
in porous rock. Energy, 29, 1445–1456.
https://doi.org/10.1016/j.energy.2004.03.078
Sun, A. Y., & Nicot, J. P. (2012). Inversion of pressure anomaly data for detecting leakage at geologic
carbon sequestration sites. Adv. Wat. Resour., 44, 20–29.
https://doi.org/10.1016/j.advwatres.2012.04.006
Teatini, P., Castelletto, N., Ferronato, M., Gambolati, G., Janna, C., Cairo, E., Marzorati, D., Colombo, D.,
Ferretti, A., Bagliani, A., & Bottazzi, F. (2011). Geomechanical response to seasonal gas storage in
depleted reservoirs: A case study in the Po River basin, Italy. J. Geophys. Res., 116, F02002.
Tiwari, D. K., Jha, B., Kundu, B., Gahalaut, V. K., & Vissa, N. K. (2021). Groundwater extraction-induced
seismicity around Delhi region, India. Scientific Reports , 11.
Tran, M., & Jha, B. (2021). Effect of poroelastic coupling and fracture dynamics on solute transport and
geomechanical stability. Wat. Resour. Res., 57(10), e2021WR029584.
https://doi.org/https://doi.org/10.1029/2021WR029584
van Genuchten, M. T. (1980). A closed-form equation for predicting the hydraulic conductivity of
unsaturated soils. Soil Sci. Soc. Am., 44, 892–898.
Vasco, D. W., Rucci, A., Ferretti, A., Novali, F., Bissell, R. C., Ringrose, P. S., Mathieson, A. S., &
Wright, I. W. (2010). Satellite-based measurements of surface deformation reveal fluid flow
associated with the geological storage of carbon dioxide. Geophys. Res. Lett., 37, L03303.
https://doi.org/10.1029/2009GL041544
Vasco, D. W., Bissell, R. C., Bohloli, B., Daley, T. M., Ferretti, A., Foxall, W., Goertz-Allmann, B. P.,
Korneev, V., Morris, J. P., Oye, V., Ramirez, A., Rinaldi, A. P., Rucci, A., Rutqvist, J., White, J., &
Zhang, R. (2018). Monitoring and Modeling Caprock Integrity at the In Salah Carbon Dioxide
Storage Site, Algeria. Geological Carbon Storage (pp. 243–269). American Geophysical Union
(AGU). https://doi.org/10.1002/9781119118657.ch12
Vasylkivska, V., King, S., Bacon, D., Harp, D., Lindner, E., Baek, S., Keating, E., Yang, Y., Zhang, Y.,
Chen, B., & Mansoor, K. (2022). NRAP-Open-IAM User’s Guide [Release alpha 2.5.0-22.03.10].
https://gitlab.com/NRAP/OpenIAM/-/raw/master/User_Guide.pdf
Vasylkivska, V., Dilmore, R., Lackey, G., Zhang, Y., King, S., Bacon, D., Chen, B., Mansoor, K., & Harp, D.
(2021). NRAP-open-IAM: A flexible open-source integrated-assessment-model for geologic
carbon storage risk assessment and management. Environmental Modelling & Software, 143,
105114. https://doi.org/10.1016/j.envsoft.2021.105114
Verdon, J. P., Kendall, J.-M., Stork, A. L., Chadwick, R. A., White, D. J., & Bissell, R. C. (2013). Comparison
of geomechanical deformation induced by megatonne-scale CO
2
storage at Sleipner, Weyburn,
and In Salah. Proc. Natl. Acad. Sci. U.S.A., 110, E2762–E2771.
Vilarrasa, V., Makhnenko, R., & Gheibi, S. (2016). Geomechanical analysis of the influence of CO
2
injection location on fault stability. J Rock Mech. Geotech. Eng., 8, 805–818.
https://doi.org/10.1029/2002JB002141
121
Vilarrasa, V., & Rutqvist, J. (2017). Thermal effects on geologic carbon storage. Earth-Science Reviews, 165,
245–256. https://doi.org/10.1016/j.earscirev.2016.12.011
White, J. A., Castelletto, N., Klevtsov, S., Bui, Q. M., Osei-Kuffuor, D., & Tchelepi, H. A. (2019). A two-stage
preconditioner for multiphase poromechanics in reservoir simulation. Computer Methods in
Applied Mechanics and Engineering, 357, 112575. https://doi.org/10.1016/j.cma.2019.112575
White, J. A., & Foxall, W. (2016). Assessing induced seismicity risk at CO
2
storage projects: Recent
progress and remaining challenges. International Journal of Greenhouse Gas Control, 49, 413–424.
https://doi.org/10.1016/j.ijggc.2016.03.021
White, J. A., Lu, C., & Chiaramonte, L. (2017). A Reduced-Order Model of Fault Leakage for
Second-Generation Toolset [Report Number: NRAP-TRS-III-012-2017].
https://doi.org/10.18141/1433155
Wu, L., Thorsen, R., Ottesen, S., Meneguolo, R., Hartvedt, K., Ringrose, P., & Nazarian, B. (2021).
Significance of fault seal in assessing CO
2
storage capacity and containment risks – an example
from the Horda Platform, northern North Sea. Petroleum Geoscience, 27(3).
https://doi.org/10.1144/petgeo2020-102
Yang, X., Lassen, R. N., Jensen, K. H., & Looms, M. C. (2015). Monitoring CO
2
migration in a shallow sand
aquifer using 3D crosshole electrical resistivity tomography. International Journal of Greenhouse
Gas Control, 42, 534–544. https://doi.org/10.1016/j.ijggc.2015.09.005
Yang, Z., & Juanes, R. (2018). Two sides of a fault: Grain-scale analysis of pore pressure control on fault
slip. 97, 022906. https://doi.org/10.1103/PhysRevE.97.022906
Yielding, G., Freeman, B., & Needham, D. T. (1997). Quantitative Fault Seal Prediction. AAPG Bulletin,
81(6), 897–917. https://doi.org/10.1306/522B498D-1727-11D7-8645000102C1865D
Zeidouni, M., & Pooladi-Darvish, M. (2012). Leakage characterization through above-zone pressure
monitoring: 1-inversion approach. J. Petrol. Sci. Eng., 98, 95–106.
https://doi.org/10.1016/j.petrol.2012.09.006
Zhao, X., & Jha, B. (2019). Role of well operations and multiphase geomechanics in controlling fault
stability during CO
2
storage and enhanced oil recovery. J. Geophys. Res., 124.
https://doi.org/10.1029/2019JB017298
Zhao, X., & Jha, B. (2021). A new coupled multiphase flow–finite strain deformation–fault slip framework
for induced seismicity. J. Comput. Phys., 433, 110178. https://doi.org/10.1016/j.jcp.2021.110178
Zheng, F., Jahandideh, A., Jha, B., & Jafarpour, B. (2021). Geologic CO
2
Storage Optimization under
Geomechanical Risk Using Coupled-Physics Models. International Journal of Greenhouse Gas
Control, 110, 103385. https://doi.org/10.1016/j.ijggc.2021.103385
122
Zheng, L., Apps, J. A., Spycher, N., Birkholzer, J. T., Kharaka, Y. K., Thordsen, J., Beers, S. R.,
Herkelrath, W. N., Kakouros, E., & Trautz, R. C. (2012). Geochemical modeling of changes in
shallow groundwater chemistry observed during the MSU-ZERT CO
2
injection experiment.
International Journal of Greenhouse Gas Control, 7, 202–217.
https://doi.org/10.1016/j.ijggc.2011.10.003
Zoback, M. D., & Gorelick, S. M. (2012). Earthquake triggering and large-scale geologic storage of carbon
dioxide. Proc. Natl. Acad. Sci. U.S.A., 109, 10164–10168.
Zyvoloski, G. A., Robinson, B. A., Dash, Z. V., Kelkar, S., Viswanathan, H. S., Pawar, R. J., Stauffer, P. H.,
Miller, T. A., & Chu, S. (2015). Software Users Manual (UM) for the FEHM Application Version 3.1
- 3.X [LA-UR-12-24493]. Retrieved June 6, 2020, from
https://www.lanl.gov/orgs/ees/fehm/docs/FEHM_UM_V3.3.0.pdf
Zyvoloski, G. A., Robinson, B. A., Dash, Z. V., & Trease, L. L. (1999). Models and Methods Summary for
the FEHM Application. Retrieved June 6, 2020, from
https://www.lanl.gov/orgs/ees/fehm/pdfs/fehm_mms.pdf
123
AppendixA
SupplementalInformationforChapter3
We provide these supplementary figures to show that the 217 runs chosen from the ensemble for further
analysis are unbiased in their input parameter distributions. Of these runs, 112 runs include solubility
effects and 105 do not include solubility effects. Runs with damage zone permeabilities above 10
− 13.5
m
2
are excluded to ensure a similar distribution in damage zone permeability among runs with and without
solubility.
124
18 17 16 15 14 13 12
Fault Core Perm. (log 10 m
2
)
0
1
CDF
No Solubility
Solubility
0.00 0.02 0.04 0.06 0.08 0.10
Fault Core Porosity
0
1
CDF
15.00 14.75 14.50 14.25 14.00 13.75 13.50
Damage Zone Perm. (log 10 m
2
)
0
1
CDF
0.00 0.02 0.04 0.06 0.08 0.10
Damage Zone Porosity
0
1
CDF
14.0 13.5 13.0 12.5 12.0
Shallow Aquifer Perm. (log 10 m
2
)
0
1
CDF
14.0 13.5 13.0 12.5 12.0
Deep Aquifer Perm. (log 10 m
2
)
0
1
CDF
0.1 0.2 0.3 0.4 0.5
Shallow Aquifer Porosity
0
1
CDF
0.05 0.10 0.15 0.20 0.25 0.30 0.35
Deep Aquifer Porosity
0
1
CDF
0.0 0.1 0.2 0.3 0.4 0.5
Fault Core Thickness (m)
0
1
CDF
750 1000 1250 1500 1750 2000 2250
Fault Permeable Length (m)
0
1
CDF
40 50 60 70
Damage Zone Thickness (m)
0
1
CDF
1 2 3 4 5 6
Rock Conductivity (W/m*K)
0
1
CDF
150 200 250 300 350 400
Well X-coordinate (m)
0
1
CDF
200 100 0 100 200
Well Y-coordinate (m)
0
1
CDF
7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0
Well Rate (kg/s)
0
1
CDF
1 2 3 4 5
Vol. Thermal Expan. Coeff. (
o
C
1
) 1e 5
0
1
CDF
0.200 0.225 0.250 0.275 0.300 0.325 0.350
Poisson's Ratio
0
1
CDF
5.0 7.5 10.0 12.5 15.0 17.5 20.0
Young's Modulus (GPa)
0
1
CDF
Figure S.1: Cumulative distribution functions of the input parameters for simulations with and without
solubility effects.
125
0.16
-0.01 -0.01
0.01 0.08 0.18
-0.01 -0.03 -0.06 0.15
-0.08 0.03 -0.07 0.04 -0.16
0.02 0.02 -0.07 -0.06 0.20 -0.11
0.09 -0.04 0.13 -0.14 -0.07 0.07 0.08
-0.09 -0.01 -0.13 -0.14 0.04 0.00 -0.04 -0.11
0.02 0.22 0.03 -0.10 0.06 -0.05 -0.01 -0.01 0.13
-0.01 0.14 -0.06 0.05 0.06 0.01 -0.10 -0.03 0.11 0.10
0.12 0.06 0.04 0.06 -0.03 -0.21 -0.16 -0.03 0.02 -0.05 0.15
-0.20 0.00 0.06 0.06 0.03 -0.16 0.14 -0.13 0.08 0.01 -0.10 0.01
0.03 -0.01 0.01 -0.07 -0.11 -0.02 0.05 -0.06 0.17 0.15 -0.11 -0.12 0.08
0.01 0.02 -0.08 0.08 -0.08 0.15 -0.10 0.11 0.02 -0.02 -0.04 0.03 0.01 -0.04
0.09 0.19 -0.09 0.01 0.05 0.04 -0.13 -0.09 0.17 -0.02 0.08 -0.03 0.00 -0.02 -0.02
-0.03 0.01 0.03 -0.06 -0.18 -0.04 0.06 0.03 -0.05 -0.03 -0.16 0.02 0.06 0.04 0.07 -0.17
0.08 -0.12 -0.08 -0.02 -0.19 -0.06 -0.14 -0.09 0.01 0.01 0.13 0.10 0.03 -0.03 0.01 -0.01 0.06
1.00
0.75
0.50
0.25
0.00
0.25
0.50
0.75
1.00
Fault Core Porosity
Damage Zone Permeability (log
10
m
2
)
Damage Zone Porosity
Shallow Aquifer Permeability (log
10
m
2
)
Deep Aquifer Permeability (log
10
m
2
)
Shallow Aquifer Porosity
Deep Aquifer Porosity
Fault Core Thickness (m)
Fault Permeable Length (m)
Rock Conductivity (W/(m*K))
Damage Zone Thickness (m)
Well Y-coordinate (m)
Well X-coordinate (m)
Well Rate (kg/s)
Vol. Ther. Expan. Coe fficient (
o
C
-1
)
Poisson's Ratio
Young's Modulus (GPa)
Fault Core Porosity
Damage Zone Permeability (log
10
m
2
)
Damage Zone Porosity
Shallow Aquifer Permeability (log
10
m
2
)
Deep Aquifer Permeability (log
10
m
2
)
Shallow Aquifer Porosity
Deep Aquifer Porosity
Fault Core Thickness (m)
Fault Permeable Length (m)
Rock Conductivity (W/(m*K))
Damage Zone Thickness (m)
Well Y-coordinate (m)
Well X-coordinate (m)
Well Rate (kg/s)
Vol. Ther. Expan. Coefficient (
o
C
-1
)
Poisson's Ratio
Fault Core Permeability (log
10
m
2
)
Figure S.2: Pearson correlation matrix of input parameters among simulations with solubility.
126
-0.06
-0.09 0.02
0.03 0.04 0.06
-0.04 -0.04 0.04 -0.04
0.17 -0.02 0.02 0.04 0.06
0.03 -0.07 -0.09 -0.10 0.05 -0.04
0.17 0.01 0.16 0.11 0.02 0.19 -0.04
0.07 0.11 0.04 -0.10 -0.05 0.02 0.05 0.11
-0.02 -0.11 0.21 0.01 0.26 0.03 -0.03 0.01 0.03
-0.21 -0.06 0.02 0.08 -0.06 -0.10 -0.07 -0.11 0.12 -0.06
-0.05 0.08 -0.15 0.09 0.03 0.09 0.02 0.10 -0.01 0.05 0.02
-0.18 0.15 -0.05 0.08 -0.14 0.09 0.03 -0.13 0.03 -0.15 0.06 0.17
0.02 0.03 0.02 -0.06 0.01 0.13 0.05 0.21 -0.04 -0.10 -0.18 -0.01 0.19
0.02 -0.10 0.00 -0.09 0.07 -0.02 0.02 0.02 -0.11 -0.02 0.06 -0.02 -0.04 0.03
0.06 -0.05 -0.08 0.01 0.11 -0.13 0.07 -0.14 -0.16 0.04 -0.12 -0.01 -0.03 -0.02 0.00
0.12 0.05 -0.01 0.11 -0.09 -0.09 0.06 0.02 0.13 -0.13 -0.05 -0.08 0.08 -0.01 -0.24 0.01
-0.05 -0.01 -0.04 0.13 0.18 0.11 0.07 -0.01 -0.01 0.13 0.02 -0.04 0.19 0.11 0.08 0.21 -0.03
1.00
0.75
0.50
0.25
0.00
0.25
0.50
0.75
1.00
Fault Core Porosity
Damage Zone Permeability (log
10
m
2
)
Damage Zone Porosity
Shallow Aquifer Permeability (log
10
m
2
)
Deep Aquifer Permeability (log
10
m
2
)
Shallow Aquifer Porosity
Deep Aquifer Porosity
Fault Core Thickness (m)
Fault Permeable Length (m)
Rock Conductivity (W/(m*K))
Damage Zone Thickness (m)
Well Y-coordinate (m)
Well X-coordinate (m)
Well Rate (kg/s)
Vol. Ther. Expan. Coefficient (
o
C
-1
)
Poisson's Ratio
Young's Modulus (GPa)
Fault Core Porosity
Damage Zone Permeability (log
10
m
2
)
Damage Zone Porosity
Shallow Aquifer Permeability (log
10
m
2
)
Deep Aquifer Permeability (log
10
m
2
)
Shallow Aquifer Porosity
Deep Aquifer Porosity
Fault Core Thickness (m)
Fault Permeable Length (m)
Rock Conductivity (W/(m*K))
Damage Zone Thickness (m)
Well Y-coordinate (m)
Well X-coordinate (m)
Well Rate (kg/s)
Vol. Ther. Expan. Coefficient (
o
C
-1
)
Poisson's Ratio
Fault Core Permeability (log
10
m
2
)
Figure S.3: Pearson correlation matrix of input parameters among simulations without solubility.
127
Fault Core Porosity
Damage Zone Permeability (log
10
m
2
)
Damage Zone Porosity
Shallow Aquifer Permeability (log
10
m
2
)
Deep Aquifer Permeability (log
10
m
2
)
Shallow Aquifer Porosity
Deep Aquifer Porosity
Fault Core Thickness (m)
Fault Permeable Length (m)
Rock Conductivity (W/(m*K))
Damage Zone Thickness (m)
Well Y-coordinate (m)
Well X-coordinate (m)
Well Rate (kg/s)
Vol. Ther. Expan. Coefficient (
o
C
-1
)
Poisson's Ratio
Young's Modulus (GPa)
Fault Core Porosity
Damage Zone Permeability (log
10
m
2
)
Damage Zone Porosity
Shallow Aquifer Permeability (log
10
m
2
)
Deep Aquifer Permeability (log
10
m
2
)
Shallow Aquifer Porosity
Deep Aquifer Porosity
Fault Core Thickness (m)
Fault Permeable Length (m)
Rock Conductivity (W/(m*K))
Damage Zone Thickness (m)
Well Y-coordinate (m)
Well X-coordinate (m)
Well Rate (kg/s)
Vol. Ther. Expan. Coefficient (
o
C
-1
)
Poisson's Ratio
Fault Core Permeability (log
10
m
2
)
0.15
0.00 -0.00
0.02 0.08 0.19
-0.02 -0.03 -0.06 0.16
-0.07 0.04 -0.07 0.04 -0.15
0.02 0.01 -0.07 -0.06 0.20 -0.13
0.08 -0.04 0.14 -0.13 -0.08 0.07 0.08
-0.11 -0.01 -0.13 -0.13 0.04 0.01 -0.05 -0.12
0.02 0.21 0.03 -0.10 0.06 -0.06 -0.01 -0.00 0.12
-0.01 0.14 -0.07 0.04 0.05 0.02 -0.10 -0.03 0.12 0.10
0.12 0.06 0.03 0.07 -0.03 -0.19 -0.16 -0.03 0.02 -0.05 0.16
-0.16 0.02 0.06 0.07 0.02 -0.17 0.14 -0.13 0.07 0.01 -0.09 0.02
0.03 -0.00 0.01 -0.08 -0.11 -0.03 0.04 -0.08 0.17 0.15 -0.11 -0.12 0.08
-0.00 0.01 -0.08 0.08 -0.08 0.17 -0.12 0.11 0.03 -0.03 -0.05 0.04 0.01 -0.04
0.08 0.18 -0.07 0.02 0.07 0.06 -0.14 -0.10 0.17 -0.01 0.09 -0.04 -0.01 -0.02 -0.02
-0.03 0.02 0.02 -0.07 -0.16 -0.07 0.07 0.02 -0.05 -0.03 -0.16 0.02 0.07 0.04 0.07 -0.17
0.08 -0.12 -0.08 -0.03 -0.18 -0.05 -0.12 -0.07 0.01 0.00 0.14 0.10 0.03 -0.01 0.01 -0.01 0.06
1.00
0.75
0.50
0.25
0.00
0.25
0.50
0.75
1.00
Figure S.4: Spearman correlation matrix of input parameters among simulations with solubility.
128
-0.06
-0.10 0.02
0.05 0.05 0.06
-0.05 -0.05 0.02 -0.03
0.19 -0.02 0.02 0.04 0.04
0.03 -0.06 -0.09 -0.12 0.06 -0.03
0.17 0.02 0.16 0.12 0.01 0.18 -0.04
0.07 0.10 0.02 -0.10 -0.05 0.04 0.05 0.11
-0.04 -0.12 0.19 0.01 0.26 0.02 -0.03 0.01 0.04
-0.21 -0.07 0.03 0.07 -0.06 -0.10 -0.07 -0.11 0.12 -0.05
-0.02 0.07 -0.14 0.10 0.03 0.10 0.02 0.09 -0.00 0.05 0.00
-0.15 0.15 -0.03 0.07 -0.14 0.11 0.03 -0.13 0.05 -0.14 0.06 0.16
0.01 0.04 0.02 -0.05 0.01 0.12 0.06 0.22 -0.04 -0.10 -0.18 -0.01 0.20
-0.00 -0.10 0.01 -0.11 0.07 -0.02 0.02 0.01 -0.10 -0.01 0.07 -0.03 -0.04 0.03
0.04 -0.06 -0.08 0.01 0.10 -0.15 0.07 -0.14 -0.17 0.05 -0.13 -0.02 -0.07 -0.03 -0.00
0.12 0.06 -0.01 0.11 -0.10 -0.09 0.07 0.02 0.14 -0.13 -0.04 -0.08 0.08 -0.02 -0.22 0.02
-0.05 -0.01 -0.03 0.16 0.19 0.10 0.07 -0.01 0.00 0.14 0.02 -0.04 0.20 0.11 0.08 0.20 -0.05
1.00
0.75
0.50
0.25
0.00
0.25
0.50
0.75
1.00
Fault Core Porosity
Damage Zone Permeability (log
10
m
2
)
Damage Zone Porosity
Shallow Aquifer Permeability (log
10
m
2
)
Deep Aquifer Permeability (log
10
m
2
)
Shallow Aquifer Porosity
Deep Aquifer Porosity
Fault Core Thickness (m)
Fault Permeable Length (m)
Rock Conductivity (W/(m*K))
Damage Zone Thickness (m)
Well Y-coordinate (m)
Well X-coordinate (m)
Well Rate (kg/s)
Vol. Ther. Expan. Coe fficient (
o
C
-1
)
Poisson's Ratio
Young's Modulus (GPa)
Fault Core Porosity
Damage Zone Permeability (log
10
m
2
)
Damage Zone Porosity
Shallow Aquifer Permeability (log
10
m
2
)
Deep Aquifer Permeability (log
10
m
2
)
Shallow Aquifer Porosity
Deep Aquifer Porosity
Fault Core Thickness (m)
Fault Permeable Length (m)
Rock Conductivity (W/(m*K))
Damage Zone Thickness (m)
Well Y-coordinate (m)
Well X-coordinate (m)
Well Rate (kg/s)
Vol. Ther. Expan. Coe fficient (
o
C
-1
)
Poisson's Ratio
Fault Core Permeability (log
10
m
2
)
Figure S.5: Spearman correlation matrix of input parameters among simulations without solubility.
129

Abstract (if available)

Abstract Climate change due to rising anthropogenic greenhouse gas emissions has become a major environmental issue over the past several decades. Geologic carbon storage is a promising technology to mitigate these emissions by sequestering CO₂ in the subsurface and away from the atmosphere for long periods of time. Geologic carbon storage requires the integrity of the subsurface storage system. Potential leakage pathways, such as faults, may compromise this integrity, allowing CO₂ to migrate out of its storage reservoir to other locations in the subsurface and possibly even the atmosphere. Other risks also exist, such as seismicity and groundwater impacts. Assessing the risk of a geologic carbon storage operation requires understanding the physical phenomena involved, the leakage pathways which may exist, and the additional risks such as induced seismicity. These risks should be assessed in an integrated manner, which is commonly done using full-physics simulations. However, these simulations can be computationally costly to run in an ensemble, which is needed to quantify risk. In this dissertation, we seek to extend the physical understanding of the dynamics of CO₂ leakage through faults and produce reduced-order models (ROM) to reduce the computational cost of risk assessment. Over the course of three studies, we find that: the fault dip angle is an important factor in the leakage behavior, yielding behavior such as downward gas leakage and surface deformation effects that are dip-angle sensitive; thermal stresses during CO₂ leakage can be sufficient to cause fault destabilization if CO₂ phase change is present; CO₂ solubility is important for estimating CO₂ leakage rates but does not appear to be significant for fault stability estimates or brine leakage-rate estimation; a deep learning approach yields a high-accuracy fault-leakage ROM; and ROM performance can be evaluated over parameter spaces that have neither data nor simulation results available. Future work includes extending the ROM with parameter ranges informed by field data and incorporating fracturing, failure mechanisms, and CO₂ solubility.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

$Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model$

PDF

Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs

PDF

Deep learning for subsurface characterization and forecasting

PDF

Integrated reservoir characterization for unconventional reservoirs using seismic, microseismic and well log data

PDF

Deep learning architectures for characterization and forecasting of fluid flow in subsurface systems

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

$Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing$

PDF

Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing

PDF

Subsurface model calibration for complex facies models

PDF

Optimization of coupled CO₂ sequestration and enhanced oil recovery

PDF

Stress and deformation analysis on fluid-exposed reservoir rocks

PDF

Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter

PDF

Efficient control optimization in subsurface flow systems with machine learning surrogate models

PDF

Stress-strain characterization of complex seismic sources

PDF

Stochastic oilfield optimization under uncertainty in future development plans

PDF

Application of data-driven modeling in basin-wide analysis of unconventional resources, including domain expertise

PDF

Sparse feature learning for dynamic subsurface imaging

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Feature learning for imaging and prior model selection

PDF

Uncertainty quantification and data assimilation via transform process for strongly nonlinear problems

Asset Metadata

Creator Meguerdijian, Saro (author)

Core Title Dynamics of CO₂ leakage through faults

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Degree Conferral Date 2022-12

Publication Date 09/07/2023

Defense Date 09/01/2022

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag CO₂,deep learning,fault leakage,fault reactivation,geologic carbon storage,induced seismicity,OAI-PMH Harvest,reduced-order model,sensitivity analysis,underground gas storage

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jha, Birendra (committee chair), Elsworth, Derek (committee member), Ershaghi, Iraj (committee member), Ghanem, Roger (committee member), Jafarpour, Behnam (committee member)

Creator Email saro.meguerd@gmail.com,smeguerd@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC111764651

Unique identifier UC111764651

Legacy Identifier etd-Meguerdiji-11186

Document Type Dissertation

Format application/pdf (imt)

Rights Meguerdijian, Saro

Type texts

Source 20220908-usctheses-batch-978 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu