University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Optimization of CO2 Storage Efficiency under Geomechanical Risks using
Coupled Flow-geomechanics-fracturing Model
by
Fangning Zheng
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(PETROLEUM ENGINEERING)
December 2023
Copyright 2024 Fangning Zheng

I dedicate this thesis to the
preservation of our planet’s sustainability.
ii

Acknowledgements
I would like to acknowledge the partial financial support of the United States Department
of Energy, National Energy Technology Laboratory through NETL-Penn State University
Coalition for Fossil Energy Research (UCFER, contract number DE- FE0026825). I would
also like to thank Energi Simulation for supporting the Advanced Reservoir Characterization and Forecasting Center at USC. The numerical simulations in this manuscript were
performed using the Computer Modeling Group software donated by CMG Limited. Additionally, I would like to express my gratitude for the Dr. Paul Fellowship, which has provided
invaluable support for my academic journey.
I would like to thank my advisors, Dr. Behnam Jafarpour and Dr. Birendra Jha, for
their invaluable support in shaping the direction of my research, inspiring the novelty of my
research, and providing critical insights into the content of my work.
I would also like to thank my dissertation committee, Dr. Behnam Jafarpour, Dr. Birendra Jha, and Dr. Roger Ghanem, for their invaluable guidance, expertise, and constructive
feedback throughout the process of completing this research.
Finally, I wish to express my heartfelt gratitude to my husband, Tianbo Xie, my parents,
Yongjun Fang and Dong Zheng, and my feline companions, Erya and Sansha, for their unwavering presence, support, encouragement, and affection during my academic journey. Their
faith in my abilities, constant support, and academic guidance have played an indispensable
role in helping me overcome challenges and attain my goals.
iii

Table of Contents
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x
Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxii
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxv
Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Introduction and Literature Review . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Topic 1: Well-placement Optimization of CO2 Storage under geomechanical risks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.2 Topic 2: Well-controlled Optimization under Caprock Fracturing and
CO2 Leakage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.3 Topic 3: Geological Uncertainty . . . . . . . . . . . . . . . . . . . . . 9
1.1.4 Topic 4: Surrogate Modeling of CO2 Coupled-physics Simulation . . . 10
1.2 Motivation Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3 Thesis Scope and Overview of Chapters . . . . . . . . . . . . . . . . . . . . . 15
Chapter 2: Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
iv

2.1 Fundamental Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.1 Coupled-Physics Governing Equations . . . . . . . . . . . . . . . . . 19
2.1.2 Mohr-Coulomb Failure Criterion . . . . . . . . . . . . . . . . . . . . . 21
2.1.3 Plasticity Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.1.4 Caprock Fracturing Mechanism . . . . . . . . . . . . . . . . . . . . . 25
2.2 Optimization Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1 Multi-objective Optimization . . . . . . . . . . . . . . . . . . . . . . 28
2.2.2 Controlled-optimization Based on Flow-only Simulation . . . . . . . . 31
2.2.3 Advanced Numerical Optimization Algorithms of CO2 Leakage Through
Caprocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2.4 Robust Optimization under Geologic Uncertainties . . . . . . . . . . 34
2.3 General Optimization Framework . . . . . . . . . . . . . . . . . . . . . . . . 35
2.4 Numerical Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.4.1 Model 1: Gaussian Type . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.4.2 Model 2: Multi-facies Type . . . . . . . . . . . . . . . . . . . . . . . 39
2.4.3 Model 3: Modeling with Caprock and Aquifer Layers . . . . . . . . . 41
2.5 CO2 Surrogate Modeling: U-Net Architecture . . . . . . . . . . . . . . . . . 44
2.6 General Workflow of Surrogate Modeling . . . . . . . . . . . . . . . . . . . . 47
2.6.1 Multi-model Learning . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.6.2 Multi-output Learning . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.6.3 General Workflow of Surrogate Modeling . . . . . . . . . . . . . . . . 50
Chapter 3: Well-placement Optimization of CO2 Storage under Geomechanical Risks 55
3.1 Case 3.1: Gaussian-type model . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.2 Case 3.2: Multi-facies model . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.3 Case 3.3: Robust Optimization for Multi-facies model . . . . . . . . . . . . . 61
3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
v

Chapter 4: Well-controlled Optimization and CO2 Leakage . . . . . . . . . . . . . . 67
4.1 Case 4.1: Optimization with Flow-only Simulation . . . . . . . . . . . . . . . 67
4.2 Case 4.2: Deterministic Optimization . . . . . . . . . . . . . . . . . . . . . . 70
4.3 Case 4.3: Robust Optimization under Geological Uncertainties . . . . . . . . 75
4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
Chapter 5: Surrogate Modeling of CO2 Coupled-physics Simulation . . . . . . . . . . 83
5.1 Multi-model Learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.1.1 U-Net with permeability . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.1.2 U-Net without permeability . . . . . . . . . . . . . . . . . . . . . . . 90
5.2 Multi-output Learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
Chapter 6: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Chapter 7: Ongoing and Future Works . . . . . . . . . . . . . . . . . . . . . . . . . . 103
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
A Evolutionary Multi-objective Optimization (EMO) Method . . . . . . . . . . 118
A.1 Elitist Non-dominated Sorting GA (NSGA-II) . . . . . . . . . . . . . 118
A.2 Controlled NSGA-II . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
B Preliminary Tests on Numerical Simulations and Objective Functions . . . . 123
B.1 Model 1: Gaussian-type . . . . . . . . . . . . . . . . . . . . . . . . . 123
B.2 Model 3: CO2 Storage with Caprock and Aquifer layers . . . . . . . 124
B.3 Model 3: Objective Function Behavior . . . . . . . . . . . . . . . . . 128
C Sensitivity Analysis of Numerical Simulation Model . . . . . . . . . . . . . . 133
D Sensitivity Analysis of Surrogate Models . . . . . . . . . . . . . . . . . . . . 140
D.1 Training Data Size of Multi-model Learning . . . . . . . . . . . . . . 140
vi

D.2 Training Data Size of Multi-output Learning . . . . . . . . . . . . . . 141
E Surrogate Modeling Data Generation . . . . . . . . . . . . . . . . . . . . . . 148
vii

List of Tables
1.1 CO2 injection geomechanical risk assessment. . . . . . . . . . . . . . . . . . . 4
1.2 Summary of CO2 leakage metrics used in the literature. . . . . . . . . . . . . 9
2.1 Simulation model settings for numerical model 1. . . . . . . . . . . . . . . . 41
2.2 Simulation model settings for numerical model 2. . . . . . . . . . . . . . . . 43
2.3 Geomechanical and flow properties for the multi-facies model. . . . . . . . . 43
2.4 Simulation model settings for numerical model 3. . . . . . . . . . . . . . . . 45
2.5 Modified Barton-Bandis Model Parameters. . . . . . . . . . . . . . . . . . . 46
2.6 Training hyperparameter settings. . . . . . . . . . . . . . . . . . . . . . . . . 49
4.1 Well flow rates for different initializations and the corresponding local solutions. 69
4.2 True value of the 95th percentile on CDF in Figure 4.11 for each optimization
iteration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.1 Testing errors across various training scenarios for effective plastic strain with
the same hyperparameter settings. . . . . . . . . . . . . . . . . . . . . . . . . 87
5.2 MSE losses of four training cases with different weight factor of well location
input. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
A.1 pseudocode of NSGA-II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
C.1 Aquifer grid resolution for model 1 in section 2.4.1 and model 2 in section 2.4.2.135
D.1 A summary of the Mean Squared Error (MSE) for four U-Net models including mapping to pressure, CO2 plume/saturation, vertical displacement,
and effective plastic strain, considering different training data sizes and the
presence (y) or absence (n) of permeability information. . . . . . . . . . . . . 141
viii

D.2 A summary of the Mean Squared Error (MSE) for pressure, CO2 saturation,
vertical displacement, and effective plastic strain of the multi-output U-Net
model with well location input weight of 10 (Case 2 in Table 5.2) under different training data sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
ix

List of Figures
1.1 Geomechanical Risks under different injection scenarios. . . . . . . . . . . . . 13
1.2 Simulation results of average reservoir pressure (a) and maximum reservoir
pressure (b) from flow-only simulations vs. CO2 leakage from coupled simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.1 Three types of failure mechanism, (a) tensile failure, (b) compaction failure,
(c) shear failure, (d) loads from all sides. . . . . . . . . . . . . . . . . . . . . 26
2.2 Modified Barton-Bandis model. . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.3 General Optimization Framework. . . . . . . . . . . . . . . . . . . . . . . . . 37
2.4 Dual-grid model for 3D coupled flow and geomechanics simulations. Left-hand
side is the geomechanics grid which encapsulates the reservoir domain. Righthand side is the reservoir domain grid. The two homogeneous permeability
regions on left and right sides of the heterogeneous injection region are the
aquifer regions for dissipating injection-induced pressure. . . . . . . . . . . . 40
2.5 Dual-grid model for 3D coupled flow and geomechanics simulations; (Left):
the geomechanics grid with the encapsulated reservoir grid. (Right): the
reservoir grid with the homogeneous aquifer surrounding the heterogeneous
injection region at the center. . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.6 Heterogeneous injection region reservoir properties for numerical model 2. . . 42
2.7 Numerical Simulation model for model 3. . . . . . . . . . . . . . . . . . . . . 44
2.8 U-Net architecture for mapping from 3-D input data including well locations
and permeability fields to 3-D output data including pressure, vertical displacement, CO2 plume and plastic strain fields. . . . . . . . . . . . . . . . . 47
x

2.9 A breakdown of the U-Net model’s layer specifications including output shapes
and the parameter counts for each layer. The output shape is represented in
four dimensions, capturing the number of filter layers and input dimensions
along the x-, y-, and z-directions, respectively. The kernel size is 3×3 and the
paddings are not shown in the table. . . . . . . . . . . . . . . . . . . . . . . 52
2.10 Schematic workflow to map from well location into plastic strain field through
two U-Net models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.11 General workflow of multi-model learning. . . . . . . . . . . . . . . . . . . . 54
3.1 The 3-D and 2-D views of Pareto-front (PF) for the three-well scenario. Figure (a), (b) and (c) show the solutions with non-zero effective plastic strain.
Figure (d), (e) and (f) show the solutions with zero effective plastic strain,
and the vertical axis is measured by negative of the safety factor defined in
section 2.2.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.2 For the three-well scenario, plots of well configuration on the injection layer’s
log10 permeability map, 3-D plots of gas plume, pressure (kPa), z-displacement
and effective plastic strain for the selected three cases. All plots are at year
15 (the end of injection year). Case 1 (a) and Case 2 (b) are having the same
CO2 storage amount of 5.97 mega-ton. Case 2 (b) and Case 3 (c) are having
the same level of safety factor and the effective plastic strain of both cases is
zero. The CO2 storage amount for Case 2 and 3 are 5.97 mega-ton and 5.4
mega-ton, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.3 Plot of the maximum storage vs. the number of wells with and without the
zero-failure constraint. The difference in the two curves highlights the role of
geomechanics in finding optimal well locations for maximum storage. . . . . 59
xi

3.4 The 3-D and 2-D views of Pareto-front (PF) for the four-well scenario using
model 2 in section 2.4.2. Plots (a), (b) and (c) show the solutions with nonzero effective plastic strain. Plots (d), (e) and (f) show the solutions with zero
effective plastic strain, and the vertical axis is measured by negative of the
safety factor defined in section 2.2.1. . . . . . . . . . . . . . . . . . . . . . . 60
3.5 Three simulations selected from the optimization results. Case 1 (a) and Case
2 (b) are having the same CO2 storage amount of 127 mega-ton. Case 2
(b) and Case 3 (c) are having the same level of safety factor of 0.15 and the
effective plastic strain of both cases is zero. The CO2 storage amount for Case
2 and 3 are 127 mega-ton and 77 mega-ton, respectively. . . . . . . . . . . . 61
3.6 Five realizations generated as examples using SISIM representing the uncertainty in permeability. The ensemble mean and standard deviation are also
plotted over the 100 realizations. . . . . . . . . . . . . . . . . . . . . . . . . 62
3.7 Simulation results for the ensemble mean (a) of gas plume, pressure, effective
plastic strain, and vertical displacement (top to bottom); simulation results
for Realization 91 (b) and Realization 73 (c); . . . . . . . . . . . . . . . . . . 63
3.8 Histograms of (a) CO2 storage, (b) maximum effective plastic strain, and (c)
maximum Z-displacement for the 100 individual simulations. . . . . . . . . . 64
4.1 (a) Distribution of well controls u discretized using du = 11000 representing
the solution space. (b) Total free CO2 (moles) for flow-only simulation using
the controls in (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.2 (a) Optimization objective function, (b) 3-D gas plume, and (c) the injection
rates corresponding to optimization solution with the flow-only simulation
and Initialization 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
xii

4.3 Solution space and optimization trajectory in log10 scale. The blue square,
blue curve and cross, and blue star are the initial solution, optimization trajectory, and optimal solution for initialization 1, respectively. The black square,
black curve and cross, and black star are the initial solution, optimization
trajectory, and optimal solution for initialization 2, respectively. The subplot
in the upper right corner illustrates an instance of a local minimum located
within the red circle on the response surface. . . . . . . . . . . . . . . . . . . 72
4.4 Optimization results of objective function vs. iteration number for initialization 1 (blue) and initialization 2 (black) in log10 scale. Solid lines and filled
circles are plotted against the left y-axis, while dashed lines and unfilled circles
are plotted against the right y-axis. . . . . . . . . . . . . . . . . . . . . . . . 73
4.5 Optimization results for the two initializations. Initial gas plume (plot a, d)
and injection rate (plot c, f) are compared with the optimal (local) gas plume
(plot b, e) and injection rate (plot c, f) for the two cases in year 15. The blue
gridblocks are caprock regions with fractures, visualized as regions with high
fracture permeability in the k-direction. The red blocks underneath represent
the gas plume in the injection reservoir. . . . . . . . . . . . . . . . . . . . . . 74
4.6 Plots of the initial (solid curves) and optimal (dashed curves) pressure (black
curves) and minimum effective stress (blue curves) over time for two gridblock
locations: (a) a location at well 1 with index={22, 8, 27}, and (b) a location
at well 3 with index={9, 18, 27}. . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.7 (a) Well locations on log10(permeability) map at the injection layer. (b) well
locations on the porosity map. (c) histogram of the log10 (permeability) map.
(a) Well locations on log10(permeability) map at the injection layer. (b) well
locations on the porosity map. (c) histogram of the log10 permeability map. 75
4.8 Permeability maps of ensemble mean (a) and standard deviation (b) over 100
realizations. (c) log10(permeability) maps for the first 5 realizations. . . . . . 76
xiii

4.9 (a) Robust optimization solution space and optimization trajectory in log10
scale. Linear interpolation is used to generate the smooth surface from the ∼50
filled circles. The black square, black curve and cross, and black star are the
initial solution, optimization trajectory, and optimal solution for initialization
1, respectively. (b) A zoomed-in view for the ‘plateau’ region. . . . . . . . . 78
4.10 Robust optimization results shown for initialization 1. The objective function
in log10 scale vs. the optimization iteration number. The black solid lines and
filled circles are plotted against the left y-axis, while the blue dashed lines and
unfilled circles are plotted against the right y-axis. . . . . . . . . . . . . . . . 78
4.11 (a) a box plot depicting the ensemble of objective function values at each
iteration, with black boxes plotted against the left y-axis and blue boxes
plotted against the right y-axis; (b) a CDF plot illustrating the ensemble
objective function values for each iteration. The light color box has a modified
x-axis for the ease of plotting; and (c) a magnified view of plot (b) focusing
on the x-axis range between -4.56 and -4.515. . . . . . . . . . . . . . . . . . 79
4.12 This figure showcases three plots, each displaying different variables against
time at the grid of well 3 with index={9, 18, 27}: (a) pressure, (b) minimum
effective stress, and (c) cumulative CO2 leakage amount. In all three plots, the
grey and cyan lines represent the ensemble values for initial and optimal cases,
respectively, while the black and blue lines illustrate the ensemble average
values for initial and optimal cases, respectively. . . . . . . . . . . . . . . . . 80
4.13 Robust optimization results for initialization 1. Averaged initial gas plume
(a) and injection rate (c) is compared with the averaged optimal (local) gas
plume (b) and injection rate (c) over ensemble at year 15. The blue blocks
represent caprock fractures in terms of their fracture permeability in the kdirection, and the lower red blocks represent the gas plume in the injection
reservoir. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
xiv

5.1 U-Net model prediction for pressure (with permeability as one of the input
data). The black vertical lines indicate wells on the maps. The first row plots
the wells on permeability map for each of the 11 cases. The second row shows
the 3-D pressure prediction. The third row shows the true observation of
pressure at the same layer, and the fourth row shows the difference between
predicted pressure and the true observation of pressure. . . . . . . . . . . . . 84
5.2 U-Net model prediction for vertical displacement, Dz (with permeability as
one of the input data). The black vertical lines indicate wells on the maps.
The first row plots the wells on permeability map for each of the 11 cases.
The second row shows the 3-D Dz prediction. The third row shows the true
observation of Dz at the same layer, and the fourth row shows the difference
between predicted Dz and the true observation of Dz. . . . . . . . . . . . . . 85
5.3 U-Net model prediction for CO2 Saturation (with permeability as one of the
input data). The black vertical lines indicate wells on the maps. The first
row plots the wells on permeability map for each of the 11 cases. The second
row shows the 3-D CO2 Saturation prediction. The third row shows the true
observation of CO2 Saturation at the same layer, and the fourth row shows
the difference between predicted and the true observation of CO2 Saturation. 85
5.4 Different cases on Table 5.1 of U-Net model prediction for effective plastic
strain. The green and black dots are well locations shown on the map. The
first row shows the well locations on permeability maps. The second row plots
the true observation for each of the 11 cases. The rest of the rows show the
prediction and mismatch at the first layer for each case. . . . . . . . . . . . . 88
xv

5.5 Feature maps at convolutional block no.3 with model input of well location
information concatenated with permeability maps. The first row plots the 16
layers on channel 1 for well configuration 1 with well locations plotted on the
last column. The second row shows the difference between the feature maps
at the same location between well configuration 1 and 0 (which is not shown
here). The third row shows the feature maps at the same location for well
configuration 2 and the fourth row shows the difference between the feature
maps at the same location between well configuration 1 and 2. . . . . . . . . 89
5.6 U-Net model prediction for vertical displacement, Dz (WITHOUT permeability as one of the input data). The green and black dots are well locations
shown on the map. The first row plots the wells on permeability map for
each of the 11 cases. The second row shows the Dz prediction at the injection
layer. The heterogeneous permeability region is labelled at the first plot inside
the dashed box. The third row shows the true observation of Dz at the same
layer, and the fourth row shows the difference between predicted Dz and the
true observation of Dz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.7 U-Net model prediction for CO2 Saturation (WITHOUT permeability as one
of the input data). The green and black dots are well locations shown on the
map. The first row plots the wells on permeability map for each of the 11
cases. The second row shows the CO2 Saturation prediction at the injection
layer. The heterogeneous permeability region is labelled at the first plot inside
the dashed box. The third row shows the true observation of CO2 Saturation
at the same layer, and the fourth row shows the difference between predicted
and the true observation of CO2 Saturation. . . . . . . . . . . . . . . . . . . 91
xvi

5.8 Seven examples of the 3-D pressure field mismatch for each of the four cases
in Table 5.2. The first row shows the well locations on permeability map, the
second row shows the true observations, and the third to the sixth rows shows
the pressure mismatch of the four cases. . . . . . . . . . . . . . . . . . . . . 94
5.9 Seven examples of the 3-D vertical displacement field mismatch for each of the
four cases in Table 5.2. The first row shows the well locations on permeability
map, the second row shows the true observations, and the third to the sixth
rows shows the vertical displacement mismatch of the four cases. . . . . . . . 95
5.10 Seven examples of the 3-D CO2 saturation field mismatch for each of the four
cases in Table 5.2. The first row shows the well locations on permeability
map, the second row shows the true observations, and the third to the sixth
rows shows the CO2 saturation mismatch of the four cases. . . . . . . . . . . 96
5.11 Seven examples of the 3-D effective plastic strain field mismatch for each of the
four cases in Table 5.2. The first row shows the well locations on permeability
map, the second row shows the true observations, and the third to the sixth
rows shows the effective plastic strain mismatch of the four cases. Row 3 is
empty because the model fails to predict the plastic strain output. . . . . . . 97
7.1 Physics behind the deep-learning U-Net neural network CO2 storage surrogate
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
A.1 Flowchart of multi-objective genetic algorithm. . . . . . . . . . . . . . . . . . 119
A.2 NSGA-II algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
B.1 Simulation results for a three-well CO2 injection example. Well locations
on the injection layer’s permeability map, gas plume shape after injecting a
certain mass, pressure, effective mean stress, Von Mises Stress, effective plastic
strain, and displacements in the x-, y-, and z-directions. . . . . . . . . . . . 125
xvii

B.2 Plots of (a) cumulative gas mass vs. injection time and (b) injection block
pressure vs. injection time for each of the three CO2 injection wells. The
three wells are injecting under BHP control of 24000 kPa for 15 years. . . . . 126
B.3 Plots of (a) 3D plot of effective plastic strain with a threshold > 0 (only
showing grid blocks that failed). Wells are labeled by the black stick and
numbers from ‘1’ to ‘3’. The blue stick labeled ‘a’ is used for comparison
later. (b) sliced 3D plot on x-direction. It shows that the failure region
is shaped like a bowl with larger and more widespread failure observed at
shallower depths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
B.4 Von Mises stress vs. mean effective normal stress plots showing the stress
trajectories of selected grid blocks. Marker locations correspond to time steps
in years, starting from Year 0 (labelled on the plot) to Year 15. (a) black
circle is the stress path for the reservoir top grid block of well 1. Red star
identifies the stress path of the reservoir top grid block of ‘a’ in Figure 5. Blue
square identifies the stress path of the reservoir top grid block of well 2. (b)
black square identifies the stress path of the reservoir top grid block of well
1. Black circle identifies the stress path of the reservoir injection grid block
of well 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
B.5 Preliminary test results on the injection layers for a one-well, 15-year injection
coupled flow-geomechanics-fracturing simulation. . . . . . . . . . . . . . . . . 128
xviii

B.6 CO2 Gas plumes (a,b,c,d) and kfz (e,f,g,h) at year 1, 5, 10, and 15, respectively. The 3D plots are in x − z view. In CO2 Gas plume plots (a,b,c,d),
the upper red blocks represent the gas plume in the upper aquifer, the middle
blue blocks represent the gas plume in the caprock fractures, and the lower
red blocks represent the gas plume in the injection reservoir. For plots of kfz
(e,f,g,h), the z-axis range is from 875 to 1000 meters in depth, which is the
caprock layers. The colored grid blocks represent the locations of caprock
fractures where the kfz has reached 100 mD. . . . . . . . . . . . . . . . . . . 129
B.7 Time evolution of kfz, minimum Effective Stress, and pressure for ‘CapPoint1’.130
B.8 Objective function f vs. injection rate u. (a) Before fracture, a near-linear
trend is observed. (b) After fracture opens at caprock, the trend of objective
function becomes exponential. (a) is a plot of the data points inside the
blue dashed-line box of (b) on a different y-scale. (c) Pressure at year 15
vs. injection rate u at the 27th layer (injection reservoir) for the blue dots
and at the 10th layer (the bottom layer of the caprocks) for the orange dots,
respectively. (d) Effective minimum stress at year 15 vs. injection rate u at
the 27th layer (blue dots) and the 10th layer (orange dots), respectively. (e)
Pressure vs. effective minimum stress at year 15 at the 27th layer (blue dots).
(f) Pressure vs. effective minimum stress at year 15 at the 10th layer (orange
dots). The y-axes are correlated with the same colors of points. . . . . . . . 131
B.9 Plots of pressure (a) and minimum effective stress (b) at grid {19, 15, 27} vs.
time. Different colors of lines are representing different simulation results with
the corresponding injection rate u specified with the legend on the right. . . 132
xix

C.1 (a) Plot of the pressure relative error and simulation run time for the sensitivity analysis of aquifer grid resolution for numerical model 1 (section 2.4.1).
The x-axis represents different cases of aquifer grid resolution listed in Table C.1. The left y-axis represents the percentage error relative to the finest
case (20×200). The right y-axis represents simulation CPU run time. (b) Plot
of the mean effective stress (I1) relative error from the sensitivity analysis of
aquifer grid resolution for numerical model 1. The x-axis represents different
cases of aquifer grid resolution listed in Table C.1. The y-axis represents the
percentage error relative to the finest case (20×200). . . . . . . . . . . . . . 136
C.2 Numerical Model 1 (section 2.4.1). Cross-sectional plots of (a) pressure and
(b) mean effective stress I1 at the injection layer of the reservoir for different
grid resolutions. Position of the cross-section A-B in 3D grid is shown in (c). 137
C.3 (a) Plot of pressure relative error and simulation run time for the sensitivity analysis of aquifer grid resolution for numerical model 2 (section 2.4.2).
The x-axis represents different cases of aquifer grid resolution listed in Table C.1. The left y-axis represents the percentage error relative to the finest
case (15×160). The right y-axis represents simulation CPU run time the cases.
(b) Plot of I1 relative errors for the sensitivity analysis of aquifer grid resolution for numerical model 2. The x-axis represents different cases of aquifer
grid resolution listed in Table C.1. The y-axis represents the percentage error
relative to the finest case (15×160). . . . . . . . . . . . . . . . . . . . . . . . 138
C.4 Numerical model 2 (section 2.4.2). Cross-sectional plots of (a) pressure and (b)
I1 at the injection layer of the reservoir for different grid resolutions. Position
of the cross-section A-B in the 3D grid is shown in (c) which is extended along
the x-axis from x=0 to x=15 km. Grid sizes are not to scale. . . . . . . . . . 139
xx

D.1 Multi-model learning model sensitivity analysis on training sample size. The
y-left axis is the MSE loss plotted for CO2 plume, vertical displacement, and
pressure while the y-right axis is the MSE error plotted for plastic strain. . . 142
D.2 Multi-output learning model sensitivity analysis on training sample size. The
y-left axis is the MSE loss plotted for pressure, vertical displacement, and
effective plastic strain while the y-right axis is the MSE error plotted for CO2
Saturation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
D.3 Sensitivity analysis training results of pressure. Seven examples are chosen for
3-D plotting. The first row plots the well locations on permeability map. The
second row plots the true observation of pressure. The third to the twelfth
rows plot the mismatch between true observation and prediction. . . . . . . . 144
D.4 Sensitivity analysis training results of vertical displacement. Seven examples
are chosen for 3-D plotting. The first row plots the well locations on permeability map. The second row plots the true observation. The third to the
twelfth rows plot the mismatch between true observation and prediction. . . 145
D.5 Sensitivity analysis training results of CO2 saturation. Seven examples are
chosen for 3-D plotting. The first row plots the well locations on permeability
map. The second row plots the true observation. The third to the twelfth
rows plot the mismatch between true observation and prediction. . . . . . . . 146
D.6 Sensitivity analysis training results of plastic strain. Seven examples are chosen for 3-D plotting. The first row plots the well locations on permeability
map. The second row plots the true observation. The third to the twelfth
rows plot the mismatch between true observation and prediction. . . . . . . . 147
E.1 Four quadrants on the heterogeneous permeability region. . . . . . . . . . . . 148
xxi

Nomenclature
α Biot’s coefficient
∆P Pressure difference
δ Kronecker delta
σ˙ Stress rate (with time)
ε˙e Elastic strain rate (with time)
λ Plastic multiplier
∥·∥2 Euclidean norm
σ
′
Effective stress tensor
σ Total stress tensor
ε Strain tensor
g Gravitational acceleration= 9.81 m/s2
m Model parameters
wβ Mass-flux of fluid phase β
ν Poisson Ratio
ϕ Porosity
ϕf Friction angle
ψ Dilatancy angle
ρb Bulk density
xxii

ρs Rock density
ρw Density of water
ρCO2 Density of CO2
σ1 Maximum principal stress
σ3 Minimum principal stress
σv Total volumetric stress
τ Shear strength
εp Effective plastic strain
c Cohesion
E Young’s Modulus
F body force
f Yield function
fβ Volumetric source term for phase β
G Shear Modulus
g Plastic potential function
Iinj Injectivity
k Absolute permeability
kabs Absolute permeability
kdr Drained bulk modulus
krg Relative permeability of CO2
krw Relative permeability of water
mβ Mass (per unit volume) of fluid phase β
MCO2 CO2 storage amount in mass.
xxiii

Mleak,g Size of a leakage event in analogy with the seismic moment M0
Nr Number of realizations in a subset of the ensemble
Nu Number of control variables
Nens Total number of ensemble realizations
p(x, t) Pore pressure at time t and location x
pw Water pressure
P0 Initial pressure
Pbh Bottom-hole pressure
pCO2 CO2 Pressure
Pg Gas phase pressure
Pr Reservoir pressure
q flow rate
Sg Gas saturation
Sw Water saturation
SCO2 CO2 saturation
U Displacement
u Well injection rate (m3/day)
vCO2 Velocity of CO2
Vof f Bulk volume of the offset aquifer
w Well location
z depth
xxiv

Abstract
Carbon Capture and Storage (CCS) has been recognized as one of the major mitigation
strategies for industrial CO2 emission and atmospheric CO2 accumulation. Numerical simulation of fluid flow, transport, and trapping mechanisms have been used to optimize CO2
storage in geologic formations by improving trapping efficiency through injection strategies.
Flow models, however, do not capture the geomechanical deformation that can occur during
CO2 injection, including reservoir expansion, ground surface uplift, caprock fracturing, and
induced seismicity. The coupled flow and geomechanical simulation models are increasingly
used to study the geomechanical effects during CO2 injection. An existing gap in the literature is related to the impact of geomechanical effects mentioned above on optimizing the
storage performance of geologic formations based on designing the well configuration and
control schedules.
I present a novel optimization framework for geologic CO2 storage under geomechanical risks, where coupled flow-geomechanics simulations are used to quantify the risks of
injection-induced ground surface deformation and rock failure in reservoir and caprock layers. A multi-objective optimization problem is formulated and solved to maximize CO2
storage while minimizing the two forms of geomechanical risks. Caprock failure could result
in creating a leakage pathway for CO2 to leak from the storage reservoir to shallower aquifers
and even to the atmosphere. To address this concern, I have developed a novel optimization
approach aimed at maintaining the caprock integrity during the storage of CO2 in geologic
formations under geological uncertainty. The developed workflow integrates advanced numerical optimization algorithms with coupled multiphase flow-geomechanics-fracturing models
xxv

for simulating the response of the storage reservoir to CO2 injection. Using the geomechanical response of the simulation, I define and quantify the potential caprock failure and CO2
leakage risks. An optimization formulation is used to minimize the risk of caprock fracturing
and CO2 leakage by finding the optimal distribution of dynamically changing CO2 injection
rates across several wells throughout the injection period. Multiple numerical experiments
with increasing complexity are presented to demonstrate the performance of the proposed
framework. The flow and geomechanical properties of deep saline aquifers are highly uncertain and can lead to significant variability in the geomechanical responses. To properly
account for the geomechanical risks of geologic CO2 storage, a stochastic description of rock
properties is incorporated into the optimization framework. Ensemble optimization and reduced sampling method are examined and compared for the stochastic optimization under
geologic uncertainties.
Simulating CO2 storage under geomechanical risks involves substantial computational
costs due to the coupling between multiphase flow and geomechanics. When identification
of optimal well locations is considered, the computational burden can escalate hundred-fold,
rendering industrial-scale applications of the technology impractical. In this work, we introduce a novel deep-learning framework designed to significantly reduce the computational
overhead associated with simulating and quantifying the geomechanical risks of CO2 storage.
Our approach leverages deep learning-based surrogate modeling to significantly enhance the
efficiency of coupled flow-geomechanics simulations for identifying suitable injection well locations for CO2 storage. We design a novel U-Net convolutional neural network to streamline
the process of mapping well location indices combined with permeability to the simulation
outputs. Subsequently, we train the U-Net to establish mappings from different well location
scenarios to corresponding pressure fields, CO2 saturation and geomechanical outputs, including vertical displacement and plastic strain. A coupled flow-geomechanics CO2 storage
simulation model is used to generate the training datasets. The U-Net model is subsequently
adopted as an efficient tool to replace the coupled flow-geomechanics simulation needed for
xxvi

identification of injection well locations to minimize geomechanical risks. The trained U-Net
model drastically reduces the computational cost by replacing coupled physics simulation
with a fast proxy model that can be used to approximately identify the geomechanical risk
associated with different well locations.
xxvii

Chapter 1
Introduction
1.1 Introduction and Literature Review
Carbon capture and sequestration (CCS) is proposed as a mitigation strategy for reducing
atmospheric CO2 accumulation. A deep saline aquifer is a practical environment for reliable,
long-term and large-scale CO2 storage [26]. The injected CO2 can be trapped and secured
in saline aquifers mainly through four trapping mechanisms: 1) geological or hydrodynamic
trapping, for example, due to impermeable rock layers and folded rock strata to prevent
CO2 plume from moving upwards, 2) residual trapping during displacement of resident fluids
by the injected supercritical CO2, which gets immobilized by capillary forces, 3) solubility
trapping due to dissolution of CO2 in saline water [102], and 4) ionic trapping when the
dissolved CO2 reacts with minerals and forms CO3 or HCO3 ions. When the concentration
of CO3 and HCO3 ions increases above their solubility limits, they precipitate to form CO3
minerals [26], which is referred to as mineral trapping of the injected CO2. Among the four
trapping mechanisms, geological trapping is the least secured mechanism over a long time
period e.g. a hundred year. Residual, solubility, and mineral trapping mechanisms have
increasing levels of security over ten to hundred-year timescale [10] and are considered as
permanent storing mechanisms [73]. Kumar et al. (2005) [73] conduct 3-D flow simulations
incorporating CO2 trapping mechanisms and observed that solubility and residual trapping
play major roles in immobilizing the injected CO2.
1

The problem of CO2 storage in saline aquifers has been widely studied over the past
two decades using both analytical modeling and numerical flow simulation. Pruess and
Garcia (2002) [97] develop a multiphase flow simulation model to simulate CO2 storage
in saline aquifers by ignoring chemical and mechanical effects. Nghiem et al. (2004) [89]
propose a fully coupled geochemical Equation-of-State compositional simulator, CMG GEMGHG, to simulate CO2 storage in saline aquifers. The analytical models for CO2 storage
problem also helped to develop a better understanding of specific phenomena during CO2
injection, for example, miscibility, non-Darcy flow, and reservoir heterogeneity, among others
[?]. The impact of relative permeability model and its hysteresis on CO2 storage has also
been highlighted in a number of studies, including [66] [59]. Bachu (2015) [7] summarizes the
storage capacity of a saline reservoir to be a function of reservoir geology and fluid properties
(e.g., reservoir thickness, pressure, temperature, water salinity), storage operation (e.g. well
control, well count and configuration, injection strategy), as well as regulatory constraints.
Operation parameters and well design strategies have been studied to maximize the storage
capacity of a storage formation, both in enhanced gas/oil recovery or aquifer storage projects
[18] [46] [28].
Optimization of CO2 injection has mainly focused on maximizing CO2 storage capacity
and efficiency, which depend on geologic properties and operational design. Shamshiri and
Jafarpour (2011) [112] show that a controlled CO2 injection which maximizes the sweep
efficiency of the CO2 flood helps to increase CO2 contact with brine, thereby improving
the solubility and residual trapping mechanisms and the overall aquifer storage potential.
Zhang and Agarwal (2012) [140] use Genetic Algorithm-based optimizer and TOUGH2 flow
simulator to optimize a water-alternating-gas (WAG) injection scheme and CO2 storage
for single-well injection in a homogeneous aquifer. They conclude that the WAG scheme
accelerates CO2 dissolution into water and results in less CO2 migration and the associated
risks. Heath et al. (2014) [53] examine the effect of injection well count and spacing on the
storage efficiency as a function of rock properties, boundary conditions, and brine production.
2

They observe that the storage pressure for closed boundary reservoir can be limited by
geomechanical considerations such as injection-triggered surface uplift and rock failure, and
that open boundary conditions as well as brine production can provide an order of magnitude
more space for storage. March et al. (2018) [81] study the possibility of CO2 storage in
naturally fractured reservoirs using a dual-porosity model and show that naturally fractured
reservoirs may be able to store CO2 with a properly designed injection rate.
During CO2 injection, the increase in reservoir pressure will induce stress and strain
changes in and around the injection zone and can trigger geomechanical hazards such as
fault reactivation [58] [56], ground deformation and uplift [105], as well as local fracturing
of the caprock and CO2 leakage. Table 1.1 summarizes a detailed list of geomechanical risk
factors that are associated with CO2 injection, the corresponding risk mechanisms, and the
locations where these risks happen. The suggested risk metrics that can be used to quantify
the associated risk factors are also listed.
1.1.1 Topic 1: Well-placement Optimization of CO2 Storage under
geomechanical risks
Wilson et al. (2003) [136] categorize ground displacement and induced seismicity as the
potential local risks of geologic CO2 storage. Celia et al. (2015) [23] propose injection
pressure build-up and induced microseismicity as the geomechanical effects related to CO2
storage. They state that even when the control of an injection well is regulated, the area
of injection pressure build-up could be larger than the area of CO2 plume and could cause
potential leakage far away from the injection well. Experiences from many field CO2 storage
projects suggest that significant geomechanical changes may indeed occur depending on
injection pressure as well as site geomechanical and flow conditions. Many researchers have
conducted investigations based on coupled flow-geomechanical analysis and interpretation of
microseismic events. InSAR measurements at the In Salah CO2 injection project site have
shown a 5 mm per year of surface uplift, which is correlated with the poroelastic expansion
3

Table 1.1: CO2 injection geomechanical risk assessment.
Location Risk
Factor Mechanism Risk Metric
Caprock/Injection
zone/Surrounding
formations
Fault
reactivation
Shear failure: injection induced pressure increase →
σef f decrease → fault shear
strength decrease → fault
slip
Caprock Wellbore
failure
• Borehole Instability:
◦ Rock yielding → rock
detachment from borehole wall → borehole
enlargement
◦ Fractured rock surrounding wellbore
• Casing Deformation and
Failure.
• CO2 Leakage:
◦ leakage rate
◦ CO2 injectivity
◦ ratio of CO2 escaped
◦ . . .
• Failure probability
• Slip Tendency, τ/σ′
n
• Effective plastic
strain (Mohr Coulomb
Failure Criterion)
• Permeability change
with time
Caprock/Injection
zone/Surrounding
formations
Rock
fracturing
1. Tensile fracturing
2. Expansion/contraction of reservoir during
injection/production
causes shear failure of
the formation.
Injection zone/
Surrounding formations
Ground
uplifting
Stress changes triggered by
injection induced pressure
change.
• Uplift rate
(mm/year)
• Total ground surface
uplift
Injection zone/
Surrounding formations
Induced
seismicity
Injection induced pressure
and stresses changes triggering rock’s plastic deformation.
• Effective plastic
strain (Mohr Coulomb
Failure Criterion)
• Seismic moment
4

from the injection zone and adjacent caprocks [105]. Thousands of microseismic events
were recorded during and after the injection period and the re-opening of the pre-existing
fractures were detected [117]. Micro-seismic events caused by pore pressure variations were
also induced at Weyburn oilfield in Canada during CO2 injection for enhanced oil recovery
[125].
These, and several other, field examples show the importance of modeling and monitoring
geomechanical effects during CO2 injection. The flow and geomechanical properties of a
reservoir play essential roles in designing the CO2 storage plan. Without proper management
of the injection pressure and stress state changes, engineering design and optimization of
geologic CO2 storage systems is unlikely to succeed in addressing the existing concerns. As
a result, most studies [78] [129] [130] [105] [49] [126] [58] [82] [56] [141] use coupled flow and
geomechanical simulation to assess the risks posed by CO2 sequestration projects. Goodarzi
et al. (2015) [48] investigates the thermal effects of CO2 injection using a coupled fluid flow,
heat transfer, and geomechanical simulation model. They propose optimization algorithms
to minimize the cost of CO2 injection while constraining the maximum fracture length by
either increasing the injection temperature or decreasing the injection rate. An existing
gap in the literature is related to the impact of geomechanical effects on optimizing the
storage performance of geologic formations based on designing the well configuration. A
holistic approach to simulation and optimization of geologic CO2 storage must incorporate
the underlying coupled processes to not only predict and enhance the trapping of CO2, but
also to represent potential geomechanical risks and account for geologic and other sources
of uncertainty. Past studies of coupled flow and geomechanics in the literature primarily
focus on analysis of geomechanical risks for existing wells in the field, without using the
coupled models to optimize future development plans. While geomechanical modeling and
risk assessment has been discussed in the literature, a principled approach to managing the
trade-off between storage efficiency maximization and risk minimization is not available.
Specifically, optimization of well configuration and other design aspects can fully exploit the
5

utility of coupled-physics modeling by not only maximizing the storage efficiency, but also
minimizing the geomechanical risks of the injection operation, which has become increasingly
important in recent years.
The geomechanical risk factors and mechanisms discussed above would occur in injection
layers as well as in the underburden and side-burden formations, activating the pre-existing
faults and fracturing reservoir rocks to cause induced seismicity. Ground surface displacement and regional uplifting can also be triggered by injection-induced pressurization of the
reservoir. It is considered as an enhancement of CO2 injectivity if the induced fractures are
contained within the injection reservoirs and do not grow toward the caprock [51]. However,
due to the critically stressed nature of the brittle crust [147], injection-induced pressurization
could change the stress state of the rock such that it enters a plastic regime at the small-scale
and fails along pre-existing weak planes at a larger scale, potentially causing (micro)seismic
events. It is more conservative to assess and model the plastic failure of the rocks within and
around the injection zone [144]. Management and optimization of the geologic CO2 storage
operations should take into account the geomechanical risk factors to arrive at efficient and
safe storage strategies.
In this work, we develop a novel optimization framework for geologic CO2 storage while
accounting for geomechanical risks. To this end, a 3D coupled fluid flow and geomechanics simulation model is developed to predict the CO2 storage performance in a large-scale
reservoir and to quantify the risks associated with rock deformation (surface uplifts) and
geomechanical failure. The objective is to use this model in an optimization problem formulation to identify well locations that maximize storage capacity while maintaining the
geomechanical risks below an acceptable level. Coupled flow and geomechanical simulation
is combined with a representative rock failure criterion to adopt a practical approach for
quantifying geomechanical risks. Multi-objective optimization based on Genetic Algorithm
is used to evaluate a Pareto-front for the decision space. The main contribution of the work
is to develop field development optimization workflows based on coupled-physics modeling to
6

identify well configurations that address the trade-off between storage efficiency and geomechanical risks of geologic CO2 storage. The proposed approach combines rigorous coupled
flow and geomechanical simulation with multi-objective optimization to present a range of
well configuration solutions for different levels of geomechanical risk. The developed framework also offers new insights about the effect of considering geomechanical risks in placing
CO2 injection wells in heterogeneous reservoirs. Specifically, well locations that may be
deemed as favorable without considering the geomechanical effects may result in geomechanical deformations that pose unfavorable risk levels. Hence, finding the optimal locations
for CO2 injection wells while accounting for geomechanical risks and formation heterogeneity
becomes a challenging problem that can be properly tackled only by optimization, which is
proposed in this paper as a principled approach to manage the trade-off between storage
efficiency and geomechanical risks.
1.1.2 Topic 2: Well-controlled Optimization under Caprock Fracturing
and CO2 Leakage
The geomechanical risks associated with caprock integrity could result in creation of CO2
leakage pathways from the storage reservoir to the shallower aquifers and even to the atmosphere. The related risk factors are fault reactivation, tensile fracturing, and wellbore failure
[14]. The mechanism of fault reactivation is associated with shear failure of the caprock.
When CO2 is injected into a deep underground formation, the increase in reservoir pressure
could lead to changes in effective normal stress and shear stress acting on the fault plane.
A previously existing fault may start slipping (fault reactivation) when the shear stress acting on the fault plane exceeds the shear strength of the fault [21] [58] [142] [83], leading to
induced seismicity [111] [19] and creation of permeable pathways in the caprock [37] [83]
[85] [?]. Tensile fracturing (hydraulic fracturing) happens when caprocks undergo tensile
failure due to injection-induced high pressures. The directions of fracture initiation and
7

propagation are highly dependent on the minimum principal stress directions and the heterogeneity of the reservoir [?] [15], and are important for assessing the geomechanical risks
and determining the caprock integrity. Fracture opening and propagation in the vertical
direction can potentially result in a leakage pathway toward shallower formations, including
groundwater aquifers [93]. Geomechanically-induced wellbore failure includes borehole instability and casing deformation and failure. Injection-induced pressure and stress changes
can cause rock deformation, yielding and even fracturing. This could lead to the detachment
and fracturing of the rock surrounding the borehole wall, causing borehole enlargement as
well as creating leakage pathways to shallower formations. The shear failure in the caprock
could cause shearing of the casing string, whereas caprock compaction that is caused by the
expansion of the injection layer could lead to buckling of the casing [51].
The existing literature on optimization of geologic CO2 storage does not include the
risks associated with caprock failure and the associated CO2 leakage through compromised
caprocks. In this work, we present a novel CO2 storage optimization workflow to minimize
CO2 leakage using coupled flow-geomechanics-fracturing simulation. A modified BartonBandis rock failure model [9] [123] is incorporated to model the initiation of rough fractures
and the fracture permeability change due to pressure build up from CO2 injection. The CO2
leakage metric is quantified based on the failure model. We propose a novel CO2 leakage
metric which consist of the total leakage amount from caprock fracturing and the fracturing
potential before caprock fracturing occurs.
Table 1.2 summarizes the CO2 leakage metrics used in the literature. Common leakage
metrics presented in literatures include leakage rate and injectivity. Meguerdijian and Jha
(2021a) [83] used leakage magnitude defined in Table 1.2 which is analogous to the moment
of earthquakes as the CO2 leakage metric. There is a missing relationship defined in the
literature to quantitatively describe the transition from pre-leakage event to leakage event.
In this study, we define a leakage metric that combines the caprock fracturing potential before
leakage happens as well as the leakage amount after leakage occurs. The first component
8

Table 1.2: Summary of CO2 leakage metrics used in the literature.
Risk Metric Equation References
CO2 leakage rate mtonCO2/year [38] [128] [17]
CO2 injectivity Iinj =
q
((Pbh − Pr)) =
q
∆P
[?] [115]
CO2 escape ratio
CO2 escaped from injection zone
total CO2 injected
× 100% [100] [101]
Leakage Magnitude Mleak,g(t) = Z
Voff
ϕ(Pg(x, z, t) − P0(x, z))Sg(x, z, t) dV [83]
is quantified as the stress differences between the minimum effective stress and the fracture
opening stress, which will be defined later in section 2.2.2. This measures the distance to the
fracturing threshold of the stress state for caprock. The second component is quantified as
the total CO2 leakage amount, which works as a barrier function added to the first component
due to its greater magnitude. The detailed formulation will also be described in section 2.2.2.
1.1.3 Topic 3: Geological Uncertainty
For CO2 sequestration sites such as saline aquifers, high uncertainty in the geological properties of storage reservoirs usually result from the lack of understanding and measurement
of the geologic representation at different scales [79]. Reservoir heterogeneity, anisotropy,
and lateral variation originating from uncertainties of rock composition, texture, pore structure, and rock type could also contribute to rock’s stress state and strength differences [86],
thus developing uncertainties in geomechanical properties to impact the CO2 storage capacities and costs in carbon capture and storage systems during pre-injection, injection, and
post-injection periods [6].
To properly account for the geological uncertainty, we developed an optimization framework where robust optimization under probabilistic approaches is incorporated into the proposed workflow. Sampling techniques and ensemble representations [139] are adopted with
hard data (well data at injector locations), soft data (porosity map), a log-normal permeability distribution, and the variogram information to sample a set of realizations of the
9

permeability map using sequential gaussian simulation by The Stanford Geostatistical Modeling Software [99].
1.1.4 Topic 4: Surrogate Modeling of CO2 Coupled-physics Simulation
Simulating CO2 storage under geomechanical risks frequently involves substantial computational costs attributed to the coupling between numerical simulations for flow and geomechanics. When history matching, optimization, and geological uncertainties are introduced,
the computational burden can escalate to hundreds or even thousands of times, rendering
industrial-scale applications of CO2 storage optimization using coupled models impractical.
There are mainly four surrogate modeling approaches [8] [55] [118] for alleviating computational costs in the application of reservoir simulation: 1) Statistics-based approach, this
technique relies on statistical methods to construct a response surface, effectively simplifying
the problem. It aims to alleviate both complexity and computational burdens, making it possible to predict reservoir performance under geological uncertainties and historical mapping
more efficiently [27] [36] [33] [40] [80] [35] [135] [67] [138] [68] [76] [114] [110]. This approach
can be considered reliable if the estimated values closely align with the actual values within
an acceptable margin of error [55]. 2) Reduced-physics approach. This approach simplifies
the physical models by different techniques, such as multi-fidelity grid resolution, to expedite computations [4] [41] [42] [106] [90] [16] [91]. 3) Reduced order modeling (ROM). ROM
involves projecting the high-dimensional system matrix onto a lower-dimensional space, enabling faster solving of equations while excluding irrelevant parameters [64] [50] [60] [52] [24]
[69] [45] [22] [124] [44] [114]. 4) AI-based models are surrogate models trained using machine
learning and pattern recognition techniques to map model’s input-output relationship. The
field of AI-based models can be broadly categorized into two major groups: data-driven proxy
models which are trained purely based on data acquired from physics-based simulation modeling or monitoring/historical field data [5] [2] [74], and physics-informed or physics-based
proxy models, which are proxy models trained using data while incorporating certain degrees
10

of physics into aspects such as loss functions and structures of the machine learning model
[113] [77] [84] [63] [137] [75].
Numerous studies have delved into the development of surrogate models for CO2 storage, primarily with the objective of predicting trapping mechanisms and storage outcomes
[71]. State-of-the-art neural network architectures, including convolutional neural networks,
U-Net, recurrent neural networks, LSTM, Residual neural networks, among others, have
been harnessed to forecast the spatial and temporal evolution of pressure and CO2 plume
behavior in 2-D or 3-D fields characterized by heterogeneous flow properties and fixed well
locations [134] [145]. Understanding the geomechanical responses resulting from CO2 injection is pivotal for ensuring the safety and effectiveness of injection operations within the
designated storage timeframe. Notably, researchers have made strides in developing CO2
storage surrogate models capable of predicting the spatial and temporal evolution of elastic
geomechanical responses, such as stress and strain fields [118] [119]. Furthermore, a substantial body of research has been dedicated to the development of surrogate models aimed
at facilitating the optimization of CO2 storage, whether online or offline [3][92] [107].
However, the challenge remains substantial, particularly when predicting the state variables under varying well locations. This challenge arises due to the significant spatial variability of pressure and other state variables in response to changing well locations within
a context of heterogeneous flow properties. Moreover, predicting the inelastic behavior of
geomechanical responses, such as plastic strain, is an intricate task due to the spatial and
temporal non-linearities involved [144] [143]. A critical gap in the existing literature is the
development of a surrogate model capable of predicting the spatial characteristics of state
variables, including non-linear geomechanical responses like plastic deformation, while accommodating well location changes. This model must achieve these predictions solely on the
basis of well configurations and flow property information as inputs.
In this work, we build an efficient CO2 storage surrogate model designed to facilitate
the transformation of well location data and permeability field information into key state
11

variables of interest. These state variables encompass 3-D pressure, CO2 plume, ground displacement, and plastic deformation of the storage field at the last time step. Our constructed
surrogate models offer substantial versatility and can be effectively employed in various applications, including data assimilation, field-scale optimization, and history matching. Their
primary benefit lies in their ability to significantly cut down on computational costs while
enabling the practical implementation of decision-making workflows that closely align with
real-world scenarios.
1.2 Motivation Examples
The In Salah pilot project [105] represented a pioneering endeavor to evaluate the feasibility
of injecting and storing carbon dioxide (CO2) underground within a depleted oil field. This
study delves into the valuable lessons derived from the project, where nearly 4 megatons of
CO2 were injected. The injection process induced a significant number of seismic events,
with over 9500 microseismic events recorded following the commencement of injection. Additionally, ground uplift of up to 25 mm was observed at the injection site over the course
of two years. These unexpected geomechanical responses underscore the critical importance
of geomechanical risk assessment and mitigation strategies in the context of CO2 injection
projects.
In order to comprehensively explore the influence of CO2 injection on geomechanical
responses, a numerical experiment was devised. This experiment involved the injection of
CO2 into a synthetic saline reservoir over a 15-year period, utilizing a configuration with five
injection wells operating at different injection scenarios. Figure 1.1 presents the relationship
between injection rates and geomechanical risks, where the red line corresponds to CO2
storage, while the blue and black lines depict the associated geomechanical risks, represented
by ground uplifting and the plastic strain of the rock. Our analysis reveals a notable trend:
as injection rates increase, the displacement and plastic deformation of the reservoir exhibit a
12

corresponding augmentation. This heightened geomechanical response signifies an increased
potential for the reservoir rocks to experience failure or fracturing, emphasizing the intricate
interplay between CO2 injection parameters and geomechanical risk. The goal of maximizing
CO2 storage in subsurface storage reservoirs while mitigating geomechanical risks presents
a complex and potentially conflicting set of objectives. This study aims to address this
challenge by identifying optimal combinations of injection schedules and injector locations
within a injection reservoir. In doing so, we seek to find a balance between the imperative
to enhance CO2 storage and the equally critical need to minimize geomechanical risks to the
reservoir.
Risks vs. Storage
Figure 1.1: Geomechanical Risks under different injection scenarios.
In Figure 1.2, we present a comparative analysis of our findings, focusing on the use of
flow-only simulation as a benchmark. A series of simulation runs were conducted, encompassing both flow-only and coupled simulations, with identical control schedules. The key
13

parameter under examination is the average reservoir pressure and maximum reservoir pressure, derived from the flow-only simulation results, plotted against the total CO2 leakage
amount from the coupled flow-geomechanics simulation results. Each blue dot within the
plot represents a unique simulation run.
Notably, our observations from the plot reveal an absence of any clear correlation between
the results obtained from flow-only pressure and the CO2 leakage outcomes derived from
coupled simulations. This observation strongly suggests that employing flow-only pressure
data as a surrogate to estimate geomechanical risks or CO2 leakage is a challenging endeavor.
Our study illustrates the limitations of relying solely on flow-based assessments when seeking
to predict and manage the geomechanical dynamics and potential CO2 leakage during CO2
injection and storage. Accurate prediction of geomechanical responses and risks during CO2
injection mandates the use of coupled flow-geomechanics simulations.
Figure 1.2: Simulation results of average reservoir pressure (a) and maximum reservoir
pressure (b) from flow-only simulations vs. CO2 leakage from coupled simulations.
14

1.3 Thesis Scope and Overview of Chapters
The primary objective of this thesis is to introduce an innovative and comprehensive optimization framework for enhancing the efficiency of CO2 storage while minimizing geomechanical risks and accounting for prediction uncertainty. The main research components and
workflow are outlined as follows:
1. Quantification of Geomechanical Risks:
• Assessment of geomechanical risks associated with CO2 injection, including ground
displacement, rock’s plastic deformation, and CO2 leakage.
• Development of a novel CO2 leakage risk metric integrated into a controlledoptimization workflow.
2. Building Coupled Flow-Geomechanics-Fracturing Simulation Model:
• Construction of a coupled flow-geomechanics-fracturing simulation model for CO2
storage under geomechanics.
• Implementation of a dual-grid system within the CMG-GEM simulator, facilitating separate gridding of flow and geomechanics, with grid interaction occurring
at designated locations.
• Adoption of the sequential coupling method to integrate the flow simulator with
the geomechanics simulator.
• Utilization of the Mohr-Coulomb’s Failure Criterion to represent the rock failure
model.
• Incorporation of the modified Barton-Bandis fracturing Model to simulate caprock
fracturing and CO2 leakage.
3. Development of two pivotal optimization frameworks for CO2 Storage Under Geomechanical Risks:
15

• Well-placement optimization framework to determine the optimal well locations,
aiming to maximize total CO2 storage while minimizing injection-induced ground
uplifting and rock plastic strain under a fixed injection bottom-hole pressure.
• Well-controlled optimization framework designed to identify the optimal well control schedules, with the objective of minimizing the CO2 leakage metric under a
fixed total injection volume constraint.
4. Incorporation of Geologic Uncertainty:
• Integration of geologic uncertainty considerations into the developed optimization
frameworks, resulting in the establishment of robust optimization workflows.
5. Development of a Fast Proxy Model for CO2 Storage:
• Creation of an efficient proxy model designed to assess geomechanical risks associated with different well locations, aiming to significantly reduce computational
costs associated with coupled-physics simulations.
• Implementation of a U-Net model for direct mapping from well location indices to
predictions of 3-D pressure, gas plume, vertical ground displacement, and effective
plastic strain fields.
This dissertation is organized as follows: Chapter 2 presents the methodology part of
the thesis. Chapter 2.1 shows the coupled flow-geomechanics governing equations and the
details of the failure and fracturing models. Chapter 2.2 represents the optimization formulations for a multi-objective optimization (Chapter 2.2.1), a well-controlled optimization with
flow-only simulation (Chapter 2.2.2), a well-controlled optimization with coupled simulation
and CO2 leakage (Chapter 2.2.3), and a robust optimization formulation (Chapter 2.2.4).
A detailed description of the algorithm and the pseudo-codes can be found in Appendix A.
Chapter 2.3 represents a comprehensive optimization framework of this study. Chapter 2.4
16

shows three distinct numerical simulation models built to be incorporated into the corresponding optimization workflows for a Gaussian’s type flow property model (Chapter 2.4.1),
a multi-facies model (Chapter 2.4.2), and a CO2 storage model with caprock and aquifer
layers (Chapter 2.4.3). The preliminary tests and model sensitivity analysis are also presented for each model in Appendix B and Appendix C, respectively. Chapter 2.5 shows the
U-Net architecture of a CO2 storage surrogate model and the corresponding workflow for
development of the surrogate modeling. The surrogate model sensitivity analysis and data
generation methods can be found in Appendix D and Appendix E, respectively.
Chapter 3 shows a detailed study of the numerical experiments for the well-placement
optimization using the Gaussian-type (Chapter 3.2) and multi-facies models (Chapter 3.3).
The robust optimization results (Chapter 3.4) are also presented for the multi-facies model.
Chapter 4 encompasses an extensive investigation of the well-controlled optimization with
caprock fracturing and CO2 leakage. It begins with a flow-only-based optimization (Chapter 4.1) to highlight its limitations. Following this, a deterministic optimization (Chapter 4.2)
employing a coupled-physics model is carried out, and the ensuing results are discussed. The
chapter is then concluded with the presentation of results from robust optimization under
geologic uncertainty (Chapter 4.3).
Chapter 5 provides an in-depth exploration of the surrogate modeling used to replace the
CO2 coupled-physics simulation across various scenarios of well locations. A comprehensive
summary and conclusion can be found in Chapter 6, while Chapter 7 delves into the ongoing
and future aspects of this research.
17

Chapter 2
Methodology
In this chapter, we describe the fundamental equations of coupled fluid flow and geomechanics and strategies to solve them using sequential two-way coupled approaches. The
Mohr’s Coulomb failure criterion and the plasticity model are also described. The caprock
fracturing mechanisms are illustrated and the modified Barton-Bandis model is used for the
development of the fracture permeability model. In the second part, we quantify geomechanical risks in the optimization formulation in terms of two scalar outputs of the coupled
flow-geomechanical simulations. We then formulate the multi-objective optimization problem with three objectives to be solved. We also construct the gradient-based well-controlled
optimization framework for minimizing CO2 leakage through caprock fracturing using a coupled flow-geomechanics-fracturing simulation. A new leakage metric is formulated as the
objective function in the optimization formulation. The robust optimization under geologic uncertainties are then formulated for both of the multi-objective optimization and the
controlled-optimization. In the third part, we describe the general optimization framework
of this study. Different types of numerical simulation models are described in the fourth
part. Finally, we conclude the chapter by presenting the architectural design and the general workflow of the CO2 storage surrogate model.
18

2.1 Fundamental Equations
2.1.1 Coupled-Physics Governing Equations
We begin the mathematical formulation of coupled flow and geomechanics simulations with
the linear momentum balance equations of quasi-static equilibrium under infinitesimal strains
[58] , which can be expressed as
∇ · σ + ρbg = 0 (2.1)
where σ is the Cauchy total stress tensor, g is the gravity vector, and the bulk density for
the two-phase (water-CO2) system is ρb = ϕ(ρCO2 SCO2 + ρwSw) + (1 − ϕ)ρs. Here, ϕ is the
true or Eulerian porosity, which is the ratio of the pore volume in the deformed configuration
to the bulk volume in the deformed configuration. ρCO2
, ρw and ρs are CO2, water and solid
grain densities, respectively, at reservoir conditions related to the respective fluid and solid
compressibilities, e.g. the CO2 compressibility is cCO2 = (1/ρCO2
)dρCO2 /dpCO2
. The fluid
mass-conservation equation for each fluid phase β (water or CO2) is expressed as
dmβ
dt + ∇ · wβ = ρβfβ, (2.2)
where the accumulation term dmβ/dt describes the time variation of fluid mass relative
to the motion of the solid skeleton, wβ is the mass-flux of fluid phase β relative to the solid
skeleton, and fβ is the volumetric source term for phase β due to wells. Mass (per unit bulk
volume in the reference configuration) of phase β is mβ = ρβSβϕ(1+εv), where the volumetric
strain εv is related to the change in effective volumetric stress δσ′
v
as εv = δσ′
v
/Kdr assuming
a linear elastic material (prior to failure) with the drained bulk modulus Kdr given in terms
of the drained Young’s modulus E and Poisson’s ratio ν. The fluid phase saturations satisfy
SCO2 + Sw ≡ 1. Flux wβ is related to the phase pressure pβ, density ρβ, viscosity µβ, and
relative permeability krβ through the multiphase extension of Darcy’s relation.
19

The effective stress equation is written as
σ
′ = σ − αpI, (2.3)
where σ
′
is the effective stress, σ is the total stress, α is the Biot’s coefficient of rock
saturated with water, p is the pore pressure (assumed to be saturation-weighted average of
phase pressures pβ), and I is the Identity tensor. The σ
′
is related to the elastic strain
ε
e = ε − ε
p
, where ε
p
is the plastic strain tensor, through a poroelastic constitutive model.
The constitutive equation of poroelasticity for fluid mass increment [142] can be expressed
as
dmβ
ρβ
= ( α
Kdr
dσv +
α
2
Kdr
dp + αdεp)Sβ + (Nwwdpw + NwCO2 dpCO2
) (2.4)
where Nww and NwCO2 are the water-water and water-CO2 inverse Biot moduli that
are functions of the rock and fluid compressibilities and the fluid saturations. The inverse
Biot moduli are components of the inverse Biot modulus tensor N = M−1 where M is
the multiphase Biot modulus tensor. For single-phase flow with a fluid compressibility cf ,
the inverse Biot modulus is expressed as 1
M = ϕcf +
(α−ϕ)(1−α)
Kdr
. Equation 2.4 relates the
increment in fluid mass (dmβ) to the change in total volumetric stress dσv and the change
in the fluid pressures dp. Substituting this equation in Eqn. 2.2 produces two pressure
equations for the two fluid phases. The coupled flow-geomechanics pressure equation is
different from the uncoupled (flow-only) pressure equation in two important ways: (1) a
∂σv/∂t term appears in the coupled equation that accounts for the effects of injectioninduced pore expansion near injection wells and matrix compression away from the wells on
the reservoir pressure, and (2) the rock compressibility term that multiplies ∂pβ/∂t in the
flow-only equation is now replaced by a term that combines α, Kdr, and N, which truly govern
the time-dependent compressibility of the rock-fluid system. The two pressure equations are
20

solved using an implicit finite volume method for the element-centered unknowns of nonwetting phase pressure pCO2 and the wetting phase saturation Sw for a given state of stress
and poroelastic parameters. The state of stress is obtained by solving Equation (1) using a
standard linear finite element method with nodal displacement vector unknowns and known
pressure solution from the flow module. We adopt the multiphase extension of Darcy’s law:
vβ = −
kk
r
β
µβ
(∇pβ − ρβg) (2.5)
where µβ and k
r
β
are the dynamic viscosity and the relative permeability, respectively, of
phase β in presence of other fluid phases.
The system of equations for flow and geomechanics can be solved iteratively to obtain the
pressure, saturation, displacement, stress and strain outputs. Iterative coupling at a given
time step is initialized by first solving the flow equations and passing the equivalent pressure
information to the geomechanics module to compute the stress, strain and displacement.
The solution from the geomechanics module is passed back to the flow module through
coupling variables (updated porosity due to geomechanical changes [121]). The flow module
recomputes the pressure and saturation distributions and sends the equivalent pressure to
the geomechanics module for another iteration. This iterative process continues until the
convergence criterion is met at which point the simulator marches to the next time step.
2.1.2 Mohr-Coulomb Failure Criterion
When a material behaves elastically, the relationship between stress and strain is reversible.
However, at large enough stresses, the material may enter an irreversible, plastic regime
where destressing the material will not remove all the strain in the material. The irreversible
part of the total strain is the plastic strain, which occurs after the material reaches a yield
state at a certain stress level and the rock fails. To model this, a failure criterion prescribing
the stress state at which plastic strain commences must be included along with a plastic
21

potential function to determine the amount of plastic strain and the new stress state post
failure. In this work, we assume that microseismicity is a direct manifestation of rock failure.
We use the Mohr-Coulomb failure criterion with a non-associated plastic potential function
to model injection-induced failure processes and quantify failure in terms of the effective
volumetric plastic strain. Below we provide a brief description of the model.
The Mohr-Coulomb criterion is a part of the Mohr-Coulomb theory to describe the shear
failure of soil and rock. The failure criterion represents the linear envelope obtained from
a plot of shear strength (τ ) of the rock versus the effective normal stress (σ
′
), expressed
as τ = c + σ
′
tan ϕf , where the quantity c is the cohesion of the rock and the angle ϕf is
called the angle of internal friction. The Mohr-Coulomb circle represents the stress state of
the rock using the maximum effective principal stress (σ
′
1
) and minimum effective principal
stress (σ
′
3
). Allowable stress states are those that do not intersect the Mohr–Coulomb linear
failure envelope. Stress states that describe a rock just at the failure point “touches” the
failure envelope. Stress states corresponding to Mohr circles which cross the failure line
are not allowed because failure of the rock would have occurred prior to the rock having
achieved such a stress state [146]. The Mohr Coulomb failure criterion depends on rock
properties, such as cohesion and friction angle, which explicitly appear in the definition of
failure criterion. Other geomechanical properties such as the Young’s Modulus, Poisson’s
ratio and flow properties, such as permeability and porosity which appear in coupled flow
and geomechanics equations, implicitly affect the stress and failure state of the rock. CO2
injection leads to an increase in pore pressure and reduction of effective principal stresses,
which moves the stress state (Mohr circle) closer to the Mohr-Coulomb failure envelope.
Theoretical details of the failure model and calculation of the effective volumetric plastic
strain is given in Vermeer (1984) [127]. The plastic failure model can be extended to large
deformation using the finite strain theory [142].
22

2.1.3 Plasticity Theory
To calculate the effective plastic strain (Vermeer and De Borst, 1984), we first assume from
plasticity theory that the total strain rate ( ˙ε) is the sum of an elastic ( ˙εe) and a plastic ( ˙εp)
contribution:
ε˙ = ˙εe + ˙εp (2.6)
From Hooke’s law, we have:
σ˙ = Dε˙e (2.7)
where D is the elastic constitutive matrix. Substitute Eqn. 2.6 into Eqn. 2.7, we obtain:
σ˙ = D(ε˙ − ε˙p) (2.8)
The Mohr-Coulomb criterion is written in terms of the stress component as:
τ
⋆ − σ
⋆
sin ϕf − c cos ϕf < 0 (2.9)
where the mean effective stress σ
⋆
,
σ
⋆ =
1
2
(σ
′
1 + σ
′
3
) (2.10)
is the center of the stress circle and the maximum shear stress τ
⋆
,
τ
⋆ =
1
2
(σ
′
1 + σ
′
3
) (2.11)
is the radius of the circle. In plasticity literature, the yield function (f) is commonly used
to distinguish plastic from elastic states. From Mohr-Coulomb criterion, the yield function
f is defined as:
f = τ
⋆ − σ
⋆
sin ϕf − cos ϕf (2.12)
23

where f is negative as long as the stress circle makes no contact with the Mohr-Coulomb
linear failure envelope, and it becomes zero when they touch. Hence, as long as f < 0 ,
the material is in elastic state, and when f = 0, the material enters the plastic state. In a
time-dependent simulation, where the increment in plastic strain over a time step becomes
the unknown, the condition implies,
ε˙p = 0 for f < 0 (2.13)
where the overdot symbol indicates time derivative. In plastic deformation, unlike elasticity
theory, there is no one-to-one correspondence between plastic strain and stresses. Instead,
the plastic strain rates are assumed to be derived from the plastic potential function g of
stresses as follows:
ε˙p = λ
∂g
∂σ
(2.14)
Here, λ is the plastic multiplier, which is positive if plastic deformation occurs (f = 0) and
is zero when f < 0. For granular materials, a suitable definition for g is (Radenkovic, 1961):
g = τ
⋆ − σ
⋆
sin ψ + constant (2.15)
where ψ is the dilatancy angle (which controls the amount of plastic volumetric strain developed during plastic shearing). The plastic potential g (Eqn. 2.15) resembles Mohr-Coulomb
yield function f (Eqn. 2.12), with the only difference being that the angle of friction ϕ in f
is replaced by dilatancy angle ψ. During a simulation, the calculated stresses at a time step
are used to compute the yield function f using Eqn. 2.12. If f ≥ 0, the Mohr circle touches
the failure envelope and plastic strain rate ( ˙εp) is calculated using Eqn. 2.14. The effective
plastic strain at time t + ∆t is then expressed as:
εp(t + ∆t) = εp(t) + ∆tε˙p (2.16)
24

To obtain the λ in Eqn. 2.14, we first substitute Eqn. 2.14 into Eqn. 2.17:
σ˙ = Dε˙ − λD
∂g
∂σ
(2.17)
We calculate the multiplier λ from a so-called consistency condition that a material remains
in plastic state when yielding happens, which can be expressed as
˙f =

∂f
∂σ
T
σ˙ = 0 (2.18)
By substituting Eqn. 2.17 into Eqn. 2.18, we obtain
˙f =

∂f
∂σ
T
Dε˙ − λD
∂g
∂σ

= 0 (2.19)
Solving this equation for λ, we obtain
λ =

∂f
∂σ
T Dε˙

∂f
∂σ
T D ∂g
∂σ
(2.20)
2.1.4 Caprock Fracturing Mechanism
During CO2 injection, an increase in pore pressure causes a decrease in effective stresses.
Rock will rupture when the minimum principle effective stress is reduced to be lower than
a threshold value which can be related to the rock’s fracture strength [104] and can be
measured from laboratory experiments. Three failure mechanisms [39], including tensile
failure, compaction failure, and shear failure, are illustrated in Figure 2.1. Tensile failure
happens when the block of rock is pulled apart until the effective tensile stress on the plane
of the load exceeds a critical value called the tensile strength. When rocks undergo tensile
failure, fracture usually develops in the direction normal to the tensile stress, as shown in
Figure 2.1 (a). Typical values of tensile strength can be ranged from 0 to a few megapascals. When a rock material with high porosity is undergoing compaction, grains will be

packed closer together to fill the pore spaces (Figure 2.1 (b)). A distributed shear failure
within the rock skeleton can be observed to cause pore collapse. Grain crushing will occur
if the compressive stresses become sufficiently high. In most sedimentary rocks, compressive
strength is much higher than tensile strength, which results in rocks fracturing at much lower
tensile stress values. Shear failure can develop when the shear stress reaches a critical value,
τmax (Figure 2.1 (c)). The Mohr–Coulomb failure criterion is usually used to describe the
shear failure mechanism. Figure 2.1 (d) shows a more general case where a rock containing
a failure plane oriented with the friction angle ϕf is loaded from all directions. The load
that is perpendicular to the fracture plane is called the fracture normal effective stress (σ
′
fn),
which is equivalent to the minimum principle effective stress. The stress parallel to the
fracture plane is called the fracture shear stress (τf ). In this paper, we assume the primary
fracture initiation criterion to be a tensile failure. In this case, the fracture will grow along
the direction normal to the minimum principle effective stress [131].
�!
�!"
#
�!
(a) Tensile failure (b) Compaction failure (c) Shear failure (d) Loads from all directions
Figure 2.1: Three types of failure mechanism, (a) tensile failure, (b) compaction failure, (c)
shear failure, (d) loads from all sides.
Rock fracturing tends to increase fracture permeability. In this work, the modified
Barton-Bandis model is adopted to compute the value of fracture permeability [9] [122]
[123]. The modified Barton-Bandis model [9] is considered to be more representative of
rough-walled fractures than the linear elastic fracture mechanics model (which is limited
to elliptical planar fracture shapes) and the Coulomb shear failure model (which does not
26

account for fracture spacing and length). The change of fracture permeability due to the
change of fracture normal effective stress is illustrated in Figure 2.2 During CO2 injection,
the fracture normal effective stress starts to decrease along path AB because of the pore pressure increase. The fracture opens when the fracture normal effective stress drops to a critical
value of fracture opening stress (σ
′
frs) at point B. Fracture permeability increases suddenly
to the maximum fracture permeability, kfmax, at point C. As injection continues, the fracture
normal effective stress keeps decreasing, and the fracture permeability is maintained at the
maximum permeability value along path CD. As injection stops and pore pressure starts
to dissipate, fracture normal effective stress begins to gradually increase following path CE,
and the fracture permeability (kf ) [29] can be computed as
kf = kccf
e
e0
4
≥ krcf , (2.21a)
where
e = e0 − Vj
, (2.21b)
Vj =
σ
′
fn
kni + σ
′
fn/Vm
, (2.21c)
Vm = e0
"
1 −

krcf
kccf 1/4
#
. (2.21d)
where kccf is the fracture closure permeability and krcf is the residual value of fracture
closure permeability. Due to the fracturing-induced permanent damage in the rock framework, fracture will not follow the original AB path towards complete closure. The fracture
permeability will approach to a residual value of the fracture closure permeability, krcf , which
is greater than zero. When injection starts again to cause a decrease in the fracture normal
effective stress, the fracture permeability will follow path CD back to the maximum fracture
permeability.
27

�!
�!"
#
A
B
D C
E
�!$%
#
�$&!
�!'()
Figure 2.2: Modified Barton-Bandis model.
2.2 Optimization Formulation
2.2.1 Multi-objective Optimization
This section describes the proposed CO2 storage optimization framework that includes both
storage efficiency and geomechanical risks as objective functions. We formulate and solve
a multi-objective optimization problem with the goal of maximizing CO2 storage capacity
while maintaining the geomechanical risks at an acceptable and safe level. The optimization
is performed to find the optimal well locations for CO2 injection with an allowable pressure
that result in maximizing CO2 mass injected while minimizing the geomechanical risks which
are two: the maximum vertical displacement (Uz) and the maximum plastic strain (εp),
both evaluated across spatial and temporal dimensions within the reservoir field. The multiobjective optimization problem can be formulated as
u
⋆ = argmin
u∈Θu
f(u) =



min f1(u) = −MCO2
,
min f2(u) = max (Uz),
min f3(u) = max (εp).
(2.22)
28

where u is the well locations and Θu is the feasible set for well locations. The amount
of CO2 stored (MCO2
) is defined as the total mass of CO2 injected into the reservoir. The
maximum vertical displacement measures the maximum change in the ground surface level
elevation (uplift) observed in the field due to injection. The maximum plastic strain measures
the maximum amount of irreversible rock failure in the field. The two risk measures in this
work are calculated by solving the coupled flow-geomechanical equations as discussed in the
section 2.1. The allowable values of the two risk measures are user-defined and depend on risk
tolerance and regulatory restrictions for specific injection sites [120] [57]. In this paper, we
assume the threshold for effective plastic strain is zero to avoid the initiation of plastic strain
and to reduce the potential for induced seismicity. We also define the maximum allowable
value for the change of vertical ground displacement as 5 mm/year. This value is based on
the In Salah CO2 injection project, where injection induced poroelastic expansion resulted in
a 5 mm/year ground surface uplifting [105]. The effective plastic strain magnitude measures
the plastic deformation after failure has happened. Before failure happens, i.e. in the elastic
deformation mode, we define the “safety factor” based on the proximity of the stress state of
a grid block to the Mohr-Coulomb failure criterion [29]. The safety factor is defined as the
shortest distance between the Mohr circle and the failure envelope normalized to the distance
between the center and the failure envelope. The safety factor will be used in plotting the
optimization results when the effective plastic strain is zero.
Well placement optimization for heterogeneous reservoirs is complicated and the objective functions typically tend to have multiple local optima [28] [116]. Global optimization
methods with gradient-free algorithms such as Genetic Algorithm (GA) and Particle Swarm
Optimization (PSO) have been proposed for solving the related optimization problems without considering the associated geomechanical risks. We formulate a multi-objective optimization problem to incorporate the geomechanical risks in the optimization problem and
to allow flexibility in adopting the risk level in the solution. Multi-objective optimization is
29

used to solve problems in which a number of (often conflicting) objectives must be met simultaneously. The problem is challenging since the optimal solution of an objective function
is different from those of the others. Naturally, the resulting solutions are expressed in terms
of a set of trade-off optimal solutions, popularly known as Pareto-optimal solutions [30]. The
advantage in using the multi-objective optimization is that the user is provided with a wide
range of Pareto-optimal solutions that represent the trade-off between risk measures and storage performance. The multiplicity of the solutions makes evolutionary (population-based)
algorithms, known as evolutionary multi-objective optimization (EMO) algorithms particularly suitable. In Appendix A, we discuss the approach to solve multi-objective optimization
problems (Eqn. 2.22) and provide details about the EMO method that is adopted in this
paper. We describe one of the state-of-the-art EMO solutions approaches using the genetic
algorithm, known as the Controlled Non-dominated Sorting GA or NSGA-II procedure [32].
In this paper, we implement the controlled NSGA-II to solve the multi-objective optimization problem. This choice is motivated by the derivative-free nature of the algorithm that
is suitable for discrete optimization problems (with well indices as decision variables), and
by the ease of implementation for multi-objective optimization problem. We use the built-in
multi-objective genetic algorithm ‘gamultiobj′
in MATLAB. Well locations are indexed as
decision variables. The population size is set to be 50 and the maximum generation is set to
be 300. We use the default selection (tournament selection), crossover (crossover intermediate), and mutation (adaptive feasible) functions to perform the reproduction procedure.
The elitism is controlled by pareto fraction and crowding distance options. We use a default
value of 0.35 as pareto fraction, i.e., the number of Pareto individuals will be limited in
the current population by 35% of the population size. The crowding distance function promotes the diversity along a Pareto-front by choosing individuals that have higher crowding
distances.
30

2.2.2 Controlled-optimization Based on Flow-only Simulation
To illustrate the importance of accounting for geomechanical effects, we compare optimization problems based on flow-only and coupled flow and geomechanics simulations. The
optimization problem with a flow-only model can be formulated using different performance
metrics. In our experiments, the total amount of injected CO2 is fixed; hence, a reasonable
choice for the objective function is to maximize the storage security by minimizing the total
free CO2 in the storage formation. The resulting problem formulation can be compactly
expressed as
u
⋆ = argmin
u∈Θu
f(u) = Mfree CO2
, (2.23a)
s.t. X
Nu
i=1
ui = Vtarget CO2 storage. (2.23b)
Where u is a vector of control variables, u = (uw1,t1
, . . . , uwW ,tS
), which represents surface
well flowing rates (m3/day) for well wi at time interval tj with i = 1 . . . W and j = 1 . . . S
where W is the total well number and S is the total number of the control schedule intervals,
and Θu is the feasible range of u. The number of control variables Nu depends on W and S
and is given by Nu = (W ×S). The term Mfree CO2
is defined as the total amount (in moles)
of free CO2 in supercritical state under reservoir conditions. Minimizing the free CO2 serves
the same purpose as maximizing the trapped (residual trapping) and dissolved (solubility
trapping) CO2 to enhance the storage security. We define a linear constraint for the control
variables to fix the total injection volume (Vtarget CO2 storage). Constraining the total amount
of CO2 injection is adopted to represent the target amount captured at a power plant for
injection into the storage formation.
31

2.2.3 Advanced Numerical Optimization Algorithms of CO2 Leakage
Through Caprocks
In this section, we formulate an optimization problem with the objective of reducing the
CO2 leakage metric resulting from caprock fracturing. This is achieved by manipulating
the well injection rates, which serve as the control variables. The well locations remain
constant throughout this process, and we employ a gradient-based optimizer to conduct the
optimization.
The deterministic form of the optimization problem can be formulated as
u
⋆ = argmin
u∈Θu
f(u, m) = −λ1∥σ
′
min(x, T|u,m) − σ
′
frs(x, T)∥2+λ2M2
leakage(T|u,m), (2.24a)
s.t. X
Nu
i=1
ui = Vtarget CO2 storage. (2.24b)
The symbol m represents the reservoir model parameters (reservoir flow and petrophysical
properties). The first term of Eqn. 2.24a, ∥σ
′
min(x, T|u,m)−σ
′
frs(x, T)∥2, is the Euclideannorm of stress difference between effective minimum stress (σ
′
min) and effective fracture
opening stress σ
′
frs defined in the modified Barton-Bandis model at the last time step T at
caprock layers given injection schedules u and reservoir permeability m, with a unit of kP a.
The fracture opening stress σ
′
frs is an inherent property of the rock and remains unaffected
by variations in injection schedules and flow properties. Before the injection induced caprock
fracturing happens, the larger the value of the first term is, the farther away the stress state
will be from fracturing. We put a minus sign in front of the first term so that minimizing the
minus stress difference equals maximizing the stress difference to prevent caprock fracturing.
The weight λ1 works as a parameter to scale the stress difference term. In our case, we use
λ1 = 1 for the unit of stresses to be kP a and the first term has a range of [−105
, 0]. The
second term, M2
leakage(T|u,m), is defined as the square of the total amount of CO2 leakage
in mole at the last time step T given injection schedules u and reservoir permeability m.
32

It is the cumulative sum of CO2 leaked into caprock fractures and upper aquifers due to
fracturing of the caprock. The weight λ2 serves to scale the leakage term and we use λ2 = 1
for the unit of leakage to be mole. Before caprock fracturing happens, the second term
remains zero. When caprock fracturing is triggered by CO2 injection, the second term starts
to accumulate to a large number ( 1021) so that the effect of the first term is diminished. In
this way, the second term works as a barrier function added to the first term. We not only
measure the extent of leakage caused by injection induced caprock fracturing, but we also
quantify the risk extent before fracturing occurs.
We define a linear constraint for the control variables so that the total injection volume
is set to a fixed value at each injection interval s. The amount of total injection can be
referred to as the target amount of CO2 released from a power plant and captured for longterm storage in a specific time frame. The injection duration is 15 years. In this paper, we
assign a constant injection rate for each well throughout the whole injection duration, and
u = (uw1
, uw2
, uw3
) with W = 3. In realistic scenarios, a valid injection schedule that consists
of a rate change every three or 12 months could be applied as the control variables in the
optimization problem. However, in our study, under the constraint that the total injection
volume is fixed for each time interval, and to save for computational cost on the gradient
calculation, we assume a constant injection rate for each well throughout the injection.
We use a gradient-based optimization algorithm to solve the problem. The forward finite
difference scheme is used to calculate the gradient vector ∇f(u) of the objective function as
∇f(u) = ∂f(u)
∂uw,t
=
f(uw,t + ∆) − f(uw,t)
∆
(2.25)
The control variable uw,t, which represents the well flowing rate for well w at time t is
normalized using min-max normalization. In this paper, the delta size is sensitive to the
simulation response and is set to be ∆ = 10−2
. The gradient is normalized as
∇f(u)norm =
∇f(u)
∥∇f(u)∥2
(2.26)
33

The Python Sci-py minimize package [108] [109] [72] [132] is adopted with Sequential Least
Squares Programming (SLSQP) algorithm to perform the proposed optimization with the
user defined (finite-difference method) gradient calculation and the linear constraints on the
optimization variables with an upper and a lower bound.
2.2.4 Robust Optimization under Geologic Uncertainties
The robust optimization formulation where geological uncertainties are incorporated can be
represented as the ensemble average over the model realizations, that is
u
⋆ = argmin
u∈Θu
F(u,m), (2.27a)
F(u,m) = 1
Nens
N
Xens
l=1

−λ1∥σ
′
min(x, T|u, ml) − σ
′
frs(x, T)∥2+λ2M2
leakage(T|u, ml)

,
(2.27b)
s.t. X
Nu
i=1
ui = Vtarget CO2 storage. (2.27c)
where ml
is the reservoir input parameter with l = 1 . . . Nens, where Nens is the total number
of the ensemble realizations. The control variables u which represents a vector of injection
rates is constrained by a total target CO2 storage volume.
The calculation of the expected value, F, as the performance measure presents a linear operation. The gradient of Eqn. 2.27 can then be calculated using the forward finite
difference method for each realization as
∇F(u,m) = 1
Nens
N
Xens
l=1
∂f(u, ml)
∂uw,t
=
1
Nens
N
Xens
l=1
f(uw,t + ∆, ml) − f(uw,t, ml)
∆
(2.28)
3

The gradient is also normalized using Eqn. 2.26. The computational cost of the ensemble
gradient F is given by (Nu × Nens), which can be significant if Nu and Nens are large. In
our case, we define Nens = 100 for the number of ensembles of permeability realizations and
Nu = 3 for the number of control variables to balance between computational cost and model
complexity. At each iteration of the robust optimization step, the number of objective function evaluations (i.e. number of simulations) is 100, and the number of gradient evaluations
is 300. The average simulation run time is one hour, and we use parallel computation with
40 nodes to perform the simulations.
2.3 General Optimization Framework
In this section, we describe a general optimization framework consists of well control and well
placement optimization under deterministic and robust cases to maximize CO2 storage efficiency under geomechanical risks and prediction uncertainties. Figure 2.3 shows a schematic
workflow based on the optimization problem we described. We first build a coupled flowgeomechanics-fracturing simulation model using CMG GEM. The modified Barton-Bandis
model is incorporated with Mohr-Coulomb Failure Criterion to describe the rock failure and
stress-dependent fracture permeability for accounting the CO2 storage amount and geomechanical risks induced by CO2 injection. Next, for control optimization, we formulate the
optimization problem to minimize the leakage metric defined in 2.2.3 through caprock fracturing by adjusting the well injection rates as the control variables. In this case, the well
locations are fixed, and a gradient-based optimizer is adopted to perform the optimization
and obtain multiple local solutions. Among these local solutions, we choose a desirable
solution satisfying higher-level criteria.
For well placement optimization, the problem can be defined as maximizing CO2 storage
while minimizing geomechanical risks such as rock failure and ground uplifting by finding
a preferable injector location. The plastic failure assumption is a conservative approach
35

to model the micro-seismicity potential. We adopt, without a loss of generality of our
overall optimization framework, the maximum effective plastic strain over the whole domain
as a measure of the geomechanical risks to be minimized. Our proposed framework allows
substitution of the risk measure with other site-specific metrics that may be more appropriate
for a given CO2 injection site with its unique geologic setting, geographic location, and
existing surface/subsurface infrastructure. Some alternative geomechanical risk measures
include the maximum differential uplift. Regardless of the risk measure used, a rigorous
coupled fluid flow and geomechanics simulation model should be constructed and used to
quantify the storage and geomechanical responses during CO2 injection, which constitute the
objective functions of our optimization problem. The optimization variables become well
locations (usually well indices), and we adopt a multi-objective optimization formulation
based on Genetic Algorithm to solve the optimization problem. Multiple trade-off solutions
will then be obtained to form a Pareto-front. And among the Pareto-front solutions, we
choose a desirable solution under higher-level criteria.
When geological uncertainties are presented in the model parameters, the robust optimization formulation should be employed under the control/well placement optimization
problem. Multiple realizations of the uncertain model parameters, such as porosity and
permeability fields, can be generated using appropriate geostatistical methods using hard
data (well data), soft data (seismic data), and variogram models of the field. The robust
optimization should be formulated based on a proper statistical performance measure (such
as the expected value, P10, P50, P90, etc.) of the model responses over the ensemble of the
model parameter realizations.
2.4 Numerical Simulation Model
In this section, we introduce three distinct numerical simulation models, each characterized
by unique reservoir geology features. The first model incorporates a heterogeneous Gaussian
36

Deterministic Optimization Formulation
Coupled Flow-GeomechanicsFracturing Simulation
Numerical Simulation Model (CMG)
Well Placement
Optimization
Multi-objective
Optimizer
Genetic
Algorithm
Multiple Trade-off Solutions
Choose a desirable solution
Higher Level Criteria
Gradient-based
Optimizer
Multiple local
solutions
Choose a desirable solution
Higher Level Criteria
Control Optimization
� � $
max �! � = CO" storage,
mi� �" � = ��ℎ ��
Gradient
Descent
� � = min ��
Hard Data Soft Data
Permeability Realizations �
Geostatistical
Technology
…
Variogram Model �" �, � = 1
�2
) �(�, �3)
4!
356
+
Robust Optimization Formulation
Figure 2.3: General Optimization Framework.
permeability field while maintaining homogeneous geomechanical properties. We conduct
multiple-well optimization scenarios within the framework of this model. The second example entails a more expansive, multi-facies model, marked by its heterogeneous flow and
geomechanical properties. For this model, the focus is on well-placement optimization. The
third example involves a reservoir domain segmented into three vertical regions: an upper
aquifer, a caprock layer, and an injection reservoir layer. This model is particularly configured to enable caprock fracturing, facilitating well-controlled optimization scenarios.
37

In this paper, we assume an isothermal environment for the CO2 injection and storage
processes. Geological trapping, residual trapping, and solubility trapping are considered in
the simulation model. Mineral trapping is effective in a longer time scale (over hundreds of
years) and is not considered in the model since we only model for 15-years of CO2 injection.
We use Brooks-Corey relative permeability model [13] with n = 2 to calculate the CO2-water
relative permeability. The maximum residual gas saturation (Sgrmax) is set to be 0.4. The
capillary pressure for the sand-CO2-brine system is referred to Plug and Bruining (2007) [95]
of the drainage capillary pressure curve for supercritical CO2 at 40 ◦C. The gas density is
calculated with the Peng-Robinson EOS [94]. The gas viscosity is calculated from the Jossi,
Stiel and Thodos correlation [96]. The aqueous phase density and viscosity are calculated
from the Rowe and Chou (1970) correlation [103] and the Kestin et al. (1981) correlation [70],
respectively. We also perform the preliminary tests on the three models using CMG-GEM
for the solution accuracy check in Appendix B. The purpose of the preliminary test is to
qualitatively validate results of the simulation model for the system of CO2 storage under
geomechanics.
2.4.1 Model 1: Gaussian Type
The 3D coupled fluid flow and geomechanics simulation model for this example is built using
CMG-GEM [29]. The synthetic simulation model is constructed with a dual grid system
(Figure 2.4), which allows for a geomechanics domain encapsulating a reservoir domain.
The reservoir domain is confined by over-burden, under-burden and side-burdens of the
geomechanics domain. The reservoir domain has two subregions: injection and aquifer.
The injection subregion is bounded by the two aquifer regions on the two sides along the
x-axis. The aquifer subregion permeability is constant at 20 mD to allow dissipation of
injection-induced pressure. The zero-flux boundary conditions are imposed on the reservoir
boundaries. For the geomechanics domain, the fixed displacement boundary conditions in
x− and y-directions are imposed on the four lateral boundaries and a fixed displacement
38

boundary condition is applied to the z-direction on the bottom boundary. The model settings
for both domains are summarized in Table 2.1.
The reservoir is at a depth of 1000 meters and the initial reservoir pressure follows
the pore pressure gradient and is approximately 9800 kPa. Each well is perforated in the
17th layer (injection layer depth is around 1043 m) out of 20 layers of the reservoir. The
heterogeneous permeability in the injection subregion ranges from 0.1 mD to 2000 mD. The
average permeability of the formation is 21 mD. The highest and lowest permeability values
in the injection layer are 138 mD and 0.13 mD, respectively, while the average permeability
in that layer is 7 mD. For the following numerical experiments, we assume multiple wells
are injecting CO2 into the reservoir for 15 years, each at a constant BHP of 24000 kPa.
The wells can be placed only within grid block numbers 4 through 27 in the x-direction to
minimize the interaction of wells with the aquifers on the two sides. Furthermore, we set a
distance constraint of at least 500 meters between the neighboring wells. During the 15-year
injection period, structural, residual and solubility trapping mechanisms are simulated and
quantified, and the total CO2 mass injected is quantified as the CO2 storage mass in the
first objective function in optimization formulation. For this model, the average simulation
run time on a 16 GB RAM, 2.2 GHz CPU machine is about 30 minutes. Parallel execution
is facilitated by a configuration comprising 40 nodes, with each node featuring 16 CPUs.
2.4.2 Model 2: Multi-facies Type
Figure 2.5 shows the simulation model with a geomechanics domain that encapsulates the
reservoir domain. The depth of the overburden starts from the ground surface (z = 0).
To mimic a constant pressure boundary condition at the edges of the storage volume, we
divide the reservoir region into two regions: an injection region where CO2 is injected and an
infinite-acting aquifer region along the boundary of the injection region. The aquifer region
has the grid block number in x- and y- directions from 1 to 12 and 64 to 75, which means grid
blocks with x and y numbers from 13 to 63 are inside the injection region. The bulk volumes
39

80×125 �
2×3000 �
30×100 �
2×3000 �
2×3000 �
30×100 �
4×1000 �
20×2.5 �
30×100 �
log10 permeability (mD)
30×100 �
4×1000 �
80×125 �
20×2.5 �
2×3000 �
x
Figure 2.4: Dual-grid model for 3D coupled flow and geomechanics simulations. Left-hand
side is the geomechanics grid which encapsulates the reservoir domain. Right-hand side is
the reservoir domain grid. The two homogeneous permeability regions on left and right sides
of the heterogeneous injection region are the aquifer regions for dissipating injection-induced
pressure.
of the aquifer elements are increased using keyword ‘VOLMOD’ in CMG-GEM to model a
constant pressure boundary effect on the injection region. The boundary conditions for the
whole reservoir domain (including the aquifer and the injection region) are set to be zeroflux. For geomechanics domain, fixed displacement boundary conditions are imposed on the
four sides of the geomechanics domain (in the x and y directions) and a fixed displacement
boundary condition in z-direction is imposed on the bottom boundary. The injection region is
further subdivided into homogeneous and heterogeneous regions. We assign heterogeneous
flow and geomechanical properties to the heterogeneous CO2 injection region (grid block
numbers 26-50 in x- and y-direction) using a sequential indicator simulation tool (SISIM).
Other parts of the reservoir region are set to be homogeneous with a permeability of 20
mD. Wells are only allowed to be placed in the heterogeneous region to avoid CO2 plume’s
interaction with the aquifer. Additional details about the model are summarized in Table 2.2.
The average permeability in the injection region is 22 mD.
We define four lithofacies within the heterogenous region. Properties for each facies are
listed in Table 2.3 and plotted in Figure 2.6. In this example, the optimum areal locations
to drill four CO2 injection wells are perforated at the depth of 18th layer of the reservoir.
40

Table 2.1: Simulation model settings for numerical model 1.
Wells are separated from each other by a minimum horizontal distance of 500 meters. We
assume that the wells are injecting CO2 into the reservoir for 20 years at a constant BHP
of 24000 kPa. For this model, the average simulation run time is about an hour. We use
parallel computing with 40 nodes to implement the optimization.
2.4.3 Model 3: Modeling with Caprock and Aquifer Layers
In this example, we build a synthetic 3D coupled flow-geomechanics-fracturing simulation
model in CMG-GEM. A dual-grid system (reservoir and geomechanics domain grids) is
incorporated in the model as shown in Figure 2.6. The reservoir domain is confined by
under-burden and side-burdens of the geomechanical domain to account for the pressureinduced stress changes around the injection reservoir. The reservoir domain is divided into
three vertical regions: an upper aquifer at a depth of 750 meter with thickness of 125 meter,
41

Permeability (mD)
Reservoir Domain
12×200 �
12×200 �
50×100 �
50×100 �
Heterogeneous Region Homogeneous Region
Aquifer Region 13×200 �
13×200 �
13×200 �
12×200 �
13×200 �
12×200 �
�
�
�
Figure 2.5: Dual-grid model for 3D coupled flow and geomechanics simulations; (Left): the
geomechanics grid with the encapsulated reservoir grid. (Right): the reservoir grid with the
homogeneous aquifer surrounding the heterogeneous injection region at the center.
Figure 2.6: Heterogeneous injection region reservoir properties for numerical model 2.
a caprock layer mainly consisting of shale below the upper aquifer with depth of 875 meters
and thickness of 125 meter, and an injection layer at a depth of 1000 meter with a thickness
of 50 meter. The upper aquifer and the caprock layers have a coarser vertical grid refinement
with 10 grid blocks × 25 m compared to the injection layer which has a finer vertical grid
refinement of 20 grid blocks × 2.5 m. The storage reservoir is further divided into two subregions horizontally along the x-axis: the middle injection subregion, which is the region where
injectors are placed, and the two side-aquifers which are used to allow pressure dissipation
and to minimize the effects of no-flow boundaries. The stress-dependent permeability model
42

Table 2.2: Simulation model settings for numerical model 2.
Table 2.3: Geomechanical and flow properties for the multi-facies model.
Facies-1 Facies-2 Facies-3 Facies-4
Young’s modulus (GPa) 17 18 19 20
Poisson’s Ratio 0.25 0.25 0.25 0.25
Cohesion (kPa) 2000 3000 4000 5000
Friction angle (degree) 20 19 16 14
Biot’s coefficient 0.8 0.77 0.75 0.73
Permeability (mD) 50 20 1 1e-4
Porosity 0.35 0.3 0.2 0.1
Corresponding geo-materials Berea Sandstone Less Shaly Berea Sandstone Shaly Berea Sandstone Marcellus Shale
(modified Barton-Bandis model in CMG-GEM with keyword ‘GPERMBB’) is implemented
in the caprock layers. We set a constant permeability of 20 mD in the upper aquifer and the
side-aquifers. The injection subregion has heterogeneous permeability values in the range
from 0.1 mD to 2000 mD. The initial matrix and fracture permeability for the caprock layer
is constant at 1E − 07 mD. Zero-flux boundary conditions are imposed on the reservoir
domain boundaries. For the geomechanics domain, fixed displacement boundary conditions
in x- and y-directions are imposed on the four lateral boundaries and a fixed displacement
43

boundary condition is applied on the bottom boundary. The model settings for both domains are summarized in Table 2.4. The modified Barton-Bandis model settings for the
caprock layers are shown in Figure 2.6 and Table 2.5. Since tensile failure is assumed to be
the primary fracture initiation criterion, the fracture opening stress (σ
′
frs) can be related to
as the tensile strength of the rock. The geo-material and the rock properties of the caprock
layer are referred to shale and the tensile strength of shale varies with burial depth [61] [133]
and bedding plane orientations [65]. In our model, we assume the reference fracture opening
stress σ
′
frs = 7580 kPa at the depth of 887.5 meters (grid-centered depth of the first layer of
the caprock) and the tensile strength gradient to be 1.131 kPa/m [133].
Reservoir Domain
Geomechanics Domain
x y
z
Pressure (kPa)
1×500 �
10×25 �
20×2.5 �
1×4000 �
30×100 �
1×4000 �
1×4000 �
30×100 �
1×4000 �
4×1000 �
20×2.5 �
30×100 �
30×100 �
4×1000 �
5×25 �
5×25 �
Upper Aquifer
Caprock
Injection Reservoir
Injection Reservoir
Injection
Subregion Side Aquifer
Side Aquifer
4×1000 � 30×100 � 4×1000 �
30×100 �
Injector
Figure 2.7: Numerical Simulation model for model 3.
2.5 CO2 Surrogate Modeling: U-Net Architecture
In this section, we present the architecture of the U-Net surrogate model that has been
developed for the purpose of mapping 3-D input data fields encompassing well locations and
permeability fields to 3-D output data, which includes pressure, vertical displacement, and
CO2 plume fields, at the last time step of simulation. Additionally, we outline the schematic
44

Table 2.4: Simulation model settings for numerical model 3.
workflow employed to map from well locations as input data to the ultimate generation of
a 3-D plastic strain field as the desired output data, involving an intermediate step in the
process. Figure 2.8 illustrates the U-Net architecture we have created. The model’s input
data dimension is configured as 56×56×20 in the x-, y-, and z-directions, and the output
data dimension mirrors this setup. In the input dataset, we combine the well location data
with the deterministic permeability field, represented in the fourth dimension. To encode
the well location dataset, we utilize one-hot encoding, where cells are marked “1” at the
perforation locations within a 3-D grid, while the remaining cells are designated as “0”. We
apply min-max normalization to standardize the output dataset.
The network is designed with a U-shaped architecture, comprising a contracting path
and an expansive path. The contracting path is a standard convolutional network, featuring
45

Table 2.5: Modified Barton-Bandis Model Parameters.
Depth
(m)
Initial
fracture
aperture
(m)
Initial
normal
fracture
stiffness
(kPa/m)
Fracture
opening
stress
Maximum
fracture
permeability
(md)
Fracture
closure
permeability
(md)
Residual
value of
fracture
closure
permeability
(md)
Layer 1 887.5 7580
Layer 2 912.5 7608
Layer 3 937.5 1.98E-05 6.79E+07 7637 100 100 33
Layer 4 965.5 7665
Layer 5 987.5 7693
a series of convolution operations followed by a Leaky Rectified Linear Unit (Leaky-ReLU).
Notably, a single down-sampling operation occurs at the third convolution layer within the
contracting path to reduce dimensionality. During the contraction phase, spatial information
is reduced while feature information is amplified. Conversely, the expansive pathway integrates both feature and spatial information by employing deconvolutions and concatenating
high-resolution features from the contracting path. In this case, a single deconvolutional
operation is carried out, with Leaky-ReLU serving as the activation function. In our design,
we opt for a single down-sampling and up-sampling operations on the input data. This
choice is deliberate, aimed at retaining a significant portion of the spatial information while
maintaining a manageable number of model parameters. This balance between preserving
spatial information and managing model parameters ensures a trade-off between model accuracy and training efficiency. To enhance network depth while maintaining efficiency, we
incorporate four residual blocks into the U-Net structure. These residual blocks, with skip
connections, contribute to creating a deeper convolutional neural network. This approach
has been shown to yield improved performance [25] [62]. Figure 2.9 provides an in-depth
view of the model structure. The U-Net model we have crafted comprises three primary
layer types: convolutional layers, residual layers, and deconvolutional layers. Figure 2.9
46

encapsulates the layer types, their corresponding output shapes, and the parameter count
for each layer. The complete model comprises approximately 1.5 million parameters. We
apply Leaky Rectified Linear Unit (Leaky ReLU) as the activation function after each convolutional layers to avoid the dying ReLU problem [88]. The Gaussian Error Linear Units
(GELUs) activation function is applied in the last layer to prevent negative output values
[54].
Pressure
U-Net Architecture
Input: 20x56x56x2
1 16
32
32 64
16 1
conv3d
down sampling
deconv3d
concatenate
Input fields Output: 20x56x56x1(multi-model)
64 128
20x56x56
10x28x28
CO2 Plume
res net …
64 x4 128
64
20x56x56
20x56x56
10x28x28
Dz field
Well locations
+
Permeability field
Output fields
Plastic strain
or 20x56x56x4(multi-output)
Figure 2.8: U-Net architecture for mapping from 3-D input data including well locations and
permeability fields to 3-D output data including pressure, vertical displacement, CO2 plume
and plastic strain fields.
2.6 General Workflow of Surrogate Modeling
In this study, we have trained a surrogate model to predict four simulation outputs: pressure,
vertical displacement, CO2 saturation, and plastic strain, using two distinct approaches.
The first approach is a multi-model learning method, where we train four separate U-Net
47

models, each dedicated to predicting one of the four outputs. The second approach employs
multi-output learning, where a single U-Net model is trained to predict all four outputs
simultaneously.
2.6.1 Multi-model Learning
We train three different surrogate models using the same proposed U-Net architecture, each
with one type of output data: pressure, vertical displacement, or CO2 plume/saturation.
Each model maps from well location information combined with permeability field directly
to the 3-D state variables of interest. The mapping process of the U-Net model can be
formulated as
Yˆ = Funet(x, m), Yˆ ∈ R
H×W×D, x ∈ R
H×W×D, m ∈ R
H×W×D (2.29)
where Yˆ represents the model output which has three-dimension: H, W, D that represent
the first, second, and third coordinate direction of the output field, respectively. Funet
represents the mapping process of the U-Net model. The input, x, represents the well
location field in 3-D with dimensions of H, W, D that represent the first, second, and third
coordinate direction, respectively. Input m represents the flow property or permeability field
in 3-D with dimensions of H, W, D that represent the first, second, and third coordinate
direction, respectively. The loss function of the training is defined using MSE loss as
L(Y, Yˆ ) = 1
n
Xn
i=1
(Yi − Yˆ
i)
2
(2.30)
where n is the total number of sample size, Yˆ is the predicted U-Net model output field,
and Y is the observed values. The training hyperparameter settings are listed in Table 2.6.
48

Table 2.6: Training hyperparameter settings.
Training-testing split ratio 0.8
Learning rate 1E-4
Step size 500
Batch size 8
Epochs 1000
Optimizer Adam
For predicting the 3-D plastic strain field, a direct mapping from well locations with
permeability values to the plastic strain field proves to be ineffective due to the highly localized and sparse nature of the output. We will elaborate on this challenge in the subsequent
results section. To overcome this limitation, we have devised a workflow, illustrated in Figure 2.10, that involves a two-step process. In the first step, we employ a U-Net model to
learn the mapping from well locations to pressure, while the second step involves training
another U-Net model to map the derived pressure field to the plastic strain field. It’s worth
noting that both U-Net models share the same structural architecture, differing only in their
respective inputs and outputs.
2.6.2 Multi-output Learning
In this approach, the surrogate model is trained using the U-Net architecture presented
in Figure 2.8 to map from well location data and permeability field into pressure, CO2
saturation, vertical displacement, and plastic strain. The mapping process of the multioutput learning can be formulated as
Yˆ = Funet(x, m), Yˆ ∈ R
N×H×W×D, x ∈ R
H×W×D, m ∈ R
H×W×D (2.31)
where Yˆ represents the model output which has four-dimension: N, H, W, D that represent the total number of output types, the first, second, and third coordinate direction of
49

the output field, respectively. In our case, N = 4. Funet represents the mapping process of
the U-Net model. The input, x, represents the well location field in 3-D with dimensions of
H, W, D that represent the first, second, and third coordinate direction, respectively. Input
m represents the flow property or permeability field in 3-D with dimensions of H, W, D that
represent the first, second, and third coordinate direction, respectively. The loss function of
the training is defined in Equation 2.30. The training hyperparameters are summarized in
Table 2.6.
2.6.3 General Workflow of Surrogate Modeling
A schematic workflow as depicted in Figure 2.11 is described for the construction of U-Net
models aimed at predicting various simulation results, including pressure, CO2 saturation,
vertical displacement, and effective plastic strain, from well location information and permeability field data using two distinct approaches. We start with data preparation step where
we establish a coupled-physics reservoir simulation model to facilitate data generation. Subsequently, we generate a comprehensive training dataset by executing reservoir simulations
under a variety of well configurations, resulting in 3-D outputs of pressure, CO2 saturation,
vertical displacement, and effective plastic strain. The detailed data generation process is
elaborated in Appendix E. Following data generation, we perform min-max normalization
on each of the diverse outputs to prepare them for training.
In the training step of the multi-model learning, we first construct and train three U-Net
models as illustrated on Figure 2.8 to establish the mapping from well location information
with permeability field to the final time step simulation results of pressure, CO2 saturation,
and vertical displacement. We then perform fine-tuning of the model structure and adjusting
training hyperparameters to optimize the model’s performance. Subsequently, we undertake
a sensitivity analysis which is presented in Appendix D to determine the minimum required
sample size for effective training. Moving on to the training step 2, in this step, we build and
train the same U-Net structure with input as the predicted pressure from the first U-Net
50

and output as the effective plastic strain field. The fine tuning and sensitivity analysis on
training sample size is also performed for this model.
In the training phase of the multi-output learning approach, we develop and train a single
U-Net model to establish a mapping between well location information and the permeability
field, and the resulting simulation outcomes at the final time step. These outcomes encompass pressure, CO2 saturation, vertical displacement, and plastic strain. Subsequently,
we conduct fine-tuning of model hyperparameters and perform sensitivity analysis on the
training dataset size in order to determine the optimal model configuration and identify an
efficient training dataset size.
The training performance of the two approaches are then compared and summarized in
Chapter 5.
51

# of parameters Output shape Layer type Block type
880 20) 56, 56, (16, Conv3d
13,856 20) 56, 56, (32, Conv3d Convolution
block 64 20) 56, 56, (32, GroupNorm
0 20) 56, 56, (32, LeakyReLU
27,680 20) 56, 56, (32, Conv3d Convolution
block 64 20) 56, 56, (32, GroupNorm
0 20) 56, 56, (32, LeakyReLU
55,360 10) 28, 28, (64, Conv3d Convolution
block 128 10) 28, 28, (64, GroupNorm
0 10) 28, 28, (64, LeakyReLU
110,656 10) 28, 28, (64, Conv3d Convolution
block 128 10) 28, 28, (64, GroupNorm
0 10) 28, 28, (64, LeakyReLU
# of parameters Output shape Layer type Block type
110,656 10) 28, 28, (64, Conv3d
Residual
block
128 10) 28, 28, (64, GroupNorm
0 10) 28, 28, (64, LeakyReLU
110,656 10) 28, 28, (64, Conv3d
128 10) 28, 28, (64, GroupNorm
110,656 10) 28, 28, (64, Conv3d
Residual
block
128 10) 28, 28, (64, GroupNorm
0 10) 28, 28, (64, LeakyReLU
110,656 10) 28, 28, (64, Conv3d
128 10) 28, 28, (64, GroupNorm
110,656 10) 28, 28, (64, Conv3d
Residual
block
128 10) 28, 28, (64, GroupNorm
0 10) 28, 28, (64, LeakyReLU
110,656 10) 28, 28, (64, Conv3d
128 10) 28, 28, (64, GroupNorm
110,656 10) 28, 28, (64, Conv3d
Residual
block
128 10) 28, 28, (64, GroupNorm
0 10) 28, 28, (64, LeakyReLU
110,656 10) 28, 28, (64, Conv3d
128 10) 28, 10, (64, GroupNorm
# of parameters Output shape Layer type Block type
0 10) 28, 28 (128, Upsample
Deconvolution
block
221,248 10) 28, 28 (64, Conv3d
128 10) 28, 28 (64, GroupNorm
0 10) 28, 28 (64, LeakyReLU
0 20) 56, 56 (128, Upsample
Deconvolution
block
110,624 20) 56, 56 (32, Conv3d
64 20) 56, 56 (32, GroupNorm
0 20) 56, 56 (32, LeakyReLU
0 20) 56, 56 (64, Upsample
Deconvolution
block
55,328 20) 56, 56 (32, Conv3d
64 20) 56, 56 (32, GroupNorm
0 20) 56, 56 (32, LeakyReLU
0 20) 56, 56 (32, Upsample
Deconvolution
block
27,664 20) 56, 56 (16, Conv3d
32 20) 56, 56 (16, GroupNorm
0 20) 56, 56 (16, LeakyReLU
433 20) 56, 56 (1, Conv3d
0 20) 56, 56 (1, GELU
0 20) 56, 56 (1, Decoder
1,510,673
Convolutional Layers
Residual Layers Deconvolutional Layers
Figure 2.9: A breakdown of the U-Net model’s layer specifications including output shapes
and the parameter counts for each layer. The output shape is represented in four dimensions,
capturing the number of filter layers and input dimensions along the x-, y-, and z-directions,
respectively. The kernel size is 3×3 and the paddings are not shown in the table.
52

Well locations
+
Permeability field U-Net 1
Pressure field
U-Net 2
Plastic strain field
Figure 2.10: Schematic workflow to map from well location into plastic strain field through
two U-Net models.
53

Data Preparation
Data
Build a coupled-physics
reservoir simulation for CO2
storage
Data
Generate training data by running
simulations with distinct scenarios of
well configurations.
Data
Data cleaning,
normalization.
3D Pressure Plum
e
Dz Plastic
strain
Pressure field
U-Net Architecture
Input: 56x56x20x2
1 16
32
32 64
16 1
conv3d
down sampling
deconv3d
concatenate
Output: 56x56x20x1
Input fields
64 128
56x56x20
28x28x10
CO2 Plume
res net
…
x4 64 128
64
56x56x20
56x56x20
28x28x10
Dz field
Well locations
+ Permeability field
Output fields
Multi-model Training
Train
Build & train three U-Nets for
mapping from well location to
pressure, plume, and vertical
displacement.
Well locations
+
Permeability field U-Net 1
Pressure field
U-Net 2
Plastic strain field
Train
Build & train U-Net for
mapping from predicted
pressure to plastic strain.
Multi-output Training
Train
Build & train a U-Net for mapping from
well location to pressure, plume, vertical
displacement, plastic strain.
Pressure field
U-Net Architecture
Input: 56x56x20x2
1 16
32
32 64
16 1
conv3d
down sampling
deconv3d
concatenate
Output: 56x56x20x4
Input fields
64 128
56x56x20
28x28x10
CO2 Plume
res net
…
x4 64 128
64
56x56x20
56x56x20
28x28x10
Dz field
Well locations
+ Permeability field
Output fields
Plastic strain field
Train
Fine tuning, sensitivity analysis
on training data size.
1.6E-05
1.4E-05
1.2E-05
1.0E-05
8.0E-06
6.0E-06
4.0E-06
2.0E-06
0.0E+00
3.0E-03
2.5E-03
2.0E-03
1.5E-03
1.0E-03
5.0E-04
0.0E+00
0 200 400 600 800 1000
l2 error
Training Size
Model Sensitivity Analysis
Fine Tuning,
Sensitivity Analysis
Figure 2.11: General workflow of multi-model learning. 54

Chapter 3
Well-placement Optimization of CO2 Storage under
Geomechanical Risks
3.1 Case 3.1: Gaussian-type model
We perform the CO2 storage optimization based on the general optimization framework
proposed in section 2.3 using numerical simulation model 1 defined in section 2.4.1, and
investigate a series of injection scenarios by considering two to eight injection wells with a
constant BHP of 24000 kPa for each well. Figure 3.1 shows the 3D (a) and 2D (b & c)
views of the Pareto-front for the scenario with three injection wells. Over 4000 simulations
were performed before the optimization was stopped due to lack of noticeable improvements.
The black dots show all the simulations performed during the optimization and the red dots
depict the simulations on the Pareto-front. We use the scenario with three wells as an
example to show the Pareto-front and the Pareto-fronts of other multi-well scenarios show a
similar behavior in terms of the shape of the solution space.
From the Pareto-front, we observe that the scenario with two wells has the widest span
of CO2 storage (from 1.5 to ∼9 mega-ton) while the scenario with eight wells exhibits the
narrowest span (from 6.5 to ∼9.5 mega-ton). From these results, we select three cases from
the scenario with three wells and analyze them in more detail to understand the physical
55

Figure 3.1: The 3-D and 2-D views of Pareto-front (PF) for the three-well scenario. Figure
(a), (b) and (c) show the solutions with non-zero effective plastic strain. Figure (d), (e) and
(f) show the solutions with zero effective plastic strain, and the vertical axis is measured by
negative of the safety factor defined in section 2.2.1.
mechanisms that play an important role during optimization under geomechanical risks. Figure 3.2 shows the plots of well configuration on the injection layer’s log10 permeability map,
as well as plots of reservoir gas plume, pressure, effective plastic strain, and z-displacement
(positive numbers indicate uplift). All plots are at year 15 (the end of injection year). In
Figure 3.2, Case 1 (a) and Case 2 (b) are having the same CO2 storage amount of 5.97 megaton. The effective plastic strain of Case 1 is above zero, which indicates that rock failure
happens during injection. For Case 2 and Case 3, the effective plastic strain is below zero
during injection. Case 2 (b) and Case 3 (c) are having the same level of safety factor and
the effective plastic strain of both cases is zero. The CO2 storage amount for Case 2 and 3
are 5.97 mega-ton and 5.4 mega-ton, respectively. From Figure 3.2, we observe that a higher
injection block permeability allows for a greater pressure build-up and a larger gas plume as
well as higher uplift. This is because when the BHP of the well is fixed, the injection rate
is proportional to the difference between BHP and the grid block pressure. We also observe
56

that for different well configurations, when the total CO2 storage amount is the same, the
effective plastic strain level can be different. For simulations of different well configurations
that are having the same level of effective plastic strain, the total CO2 storage amount can
be different.
Figure 3.2: For the three-well scenario, plots of well configuration on the injection layer’s
log10 permeability map, 3-D plots of gas plume, pressure (kPa), z-displacement and effective
plastic strain for the selected three cases. All plots are at year 15 (the end of injection year).
Case 1 (a) and Case 2 (b) are having the same CO2 storage amount of 5.97 mega-ton. Case
2 (b) and Case 3 (c) are having the same level of safety factor and the effective plastic strain
of both cases is zero. The CO2 storage amount for Case 2 and 3 are 5.97 mega-ton and 5.4
mega-ton, respectively.
57

In case of a deterministic permeability map, as considered in this example, the optimal
number of wells will depend on the availability of high permeability regions. In our numerical
experiments, the high permeability region could only accommodate two to three wells for
reasonable values of inter-well spacing, which means two to three CO2 injection wells in the
high permeability region are the most effective wells. From the cumulative gas mass for four-,
five-, six-, seven-, and eight-well cases, we observe that one of the wells assumes significantly
higher injection levels relative to other wells, and the corresponding well is always located
in the high permeability region. Considering the cost of drilling, the number of wells must
be determined based on economics and injection capacity of each well. The injection well
locations can be identified based on the trade-off between minimizing geomechanical risks
and maximizing CO2 storage.
Figure 3.3 shows a plot of the maximum storage vs. number of wells with and without
the no-failure constraint. The blue curve shows the maximum storage for different number of
wells without any constraints and the black curve shows the maximum storage for different
number of wells with no-failure constraint. For a scenario with Nw wells, an exhaustive
search is performed by fixing the locations of (Nw − 1) well based on optimization with
(Nw − 1) wells (with or without constraint) and considering a new infill well. With a fixed
BHP control and 15 years of injection, the maximum CO2 storage increases with the number
of wells. However, the incremental storage for three and more wells in the field is limited.
Moreover, the reservoir’s storage potential is smaller when a no-failure constraint is added,
which highlights the importance of considering geomechanical risks in evaluating the storage
potential of the field.
3.2 Case 3.2: Multi-facies model
In this section, we present the simulation and optimization results of the multi-facies model
defined in section 2.4.2. Figure 3.4 shows the 3-D and 2-D views of the Pareto-front for the
58

Figure 3.3: Plot of the maximum storage vs. the number of wells with and without the
zero-failure constraint. The difference in the two curves highlights the role of geomechanics
in finding optimal well locations for maximum storage.
scenario with four wells. From the optimization results, we select and analyze three cases
to understand the physical mechanisms that play important roles during optimization under
geomechanical risks. Figure 3.4 shows the well configuration over the log10 permeability
map in the injection layer, as well as the CO2 plume, pressure, effective plastic strain, and zdisplacement (positive numbers indicate uplift). All plots are at Year 20 (the end of injection
year).
In Figure 3.5, Case 1 (a) and Case 2 (b) have the same CO2 storage amount of 127 megaton, while the storage amount for Case 3 is 77 mega-ton. Case 2 (b) and Case 3 (c) have the
same level of safety factor of 0.15 and the effective plastic strain of both cases is zero. From
Figure 13, we observe that CO2 storage is high when the well is placed in high permeability
regions. This is consistent with the observations made in section 3.1. When the well is placed
59

in a high permeability grid block, the high injection rate leads to a greater pressure build-up,
a larger gas plume, greater uplift, lower safety factor, and higher effective plastic strain. As in
section 3.1, we observe that for simulations with different well configurations, when the total
CO2 storage amount is the same the effective plastic strain level can be different. Similarly,
when different well configurations are used with the same level of effective plastic strain,
the total CO2 storage amount can be different. Optimization can search within different
well configurations with the same level of geomechanical risks and identify the solution with
highest storage capacity. Alternatively, it can search within well configurations with the
same level of storage to identify solutions that minimize the geomechanical risk. If coupled
fluid flow and geomechanics simulation is not used in formulating the optimization problem,
the resulting well configuration solutions may lead to unacceptable risks of rock failure.
Figure 3.4: The 3-D and 2-D views of Pareto-front (PF) for the four-well scenario using
model 2 in section 2.4.2. Plots (a), (b) and (c) show the solutions with non-zero effective
plastic strain. Plots (d), (e) and (f) show the solutions with zero effective plastic strain, and
the vertical axis is measured by negative of the safety factor defined in section 2.2.1.
60

Figure 3.5: Three simulations selected from the optimization results. Case 1 (a) and Case
2 (b) are having the same CO2 storage amount of 127 mega-ton. Case 2 (b) and Case 3
(c) are having the same level of safety factor of 0.15 and the effective plastic strain of both
cases is zero. The CO2 storage amount for Case 2 and 3 are 127 mega-ton and 77 mega-ton,
respectively.
3.3 Case 3.3: Robust Optimization for Multi-facies
model
We also investigate the effect of uncertainty in rock’s hydraulic and geomechanical properties on the geomechanical response to CO2 injection. In total, 100 realizations of the aquifer
61

model are generated using the SISIM algorithm. Figure 3.6 shows five samples of the permeability models that exhibit distinct spatial variability. The ensemble mean and standard
deviation are also plotted over the 100 realizations. We perform 100 forward simulations using the well configurations of Case 2 in Figure 3.5. The realizations reflect the uncertainty in
flow properties (rock permeability and porosity) as well as geomechanical properties (Young’s
modulus, cohesion, friction angle, and Biot’s coefficient). The mean map of gas plume, pressure field, effective plastic strain, and vertical displacement over the entire ensemble are
plotted in Figure 3.6 (a). The histograms of CO2 storage, maximum effective plastic strain,
and maximum vertical displacement for the realizations are plotted in Figure 8. In Figure 3.6, we observe that the mean effective plastic strain will no longer remain zero during
injection (compared to Case 2 as the base case in 3.5). This is an expected behavior when
multiple realizations are used to represent the uncertainty in formation properties because
in some realizations the injection wells may intersect grid cells with low permeability, which
can lead to higher pressure build-up and geomechanical risks.
mean(log10(perm)) std(log10(perm))
Figure 3.6: Five realizations generated as examples using SISIM representing the uncertainty
in permeability. The ensemble mean and standard deviation are also plotted over the 100
realizations.
The injection limit and security of a CO2 storage reservoir depends on the properties of
both reservoir and overburden such as permeability, porosity, Young’s Modulus, Poisson’s
ratio, Biot coefficient, cohesion, and friction angle. In general, the geomechanical risks and
constraints can limit the storage capacity of aquifers and influence engineering design of
CO2 injection. Coupled-physics models can be used to optimize the number, locations and
62

Figure 3.7: Simulation results for the ensemble mean (a) of gas plume, pressure, effective
plastic strain, and vertical displacement (top to bottom); simulation results for Realization
91 (b) and Realization 73 (c);
bottomhole pressure of CO2 injection wells based on the flow and mechanical properties
of the storage formation. However, the flow and geomechanical properties of deep saline
aquifers are highly uncertain. To reflect this important issue, we evaluated the optimal
well configuration in our examples using 100 realizations to represent the uncertainty in
flow and geomechanical properties. We observed significant variability in the geomechanical
response of the formation and the occurrence of failure in multiple realizations. Therefore, to
63

Figure 3.8: Histograms of (a) CO2 storage, (b) maximum effective plastic strain, and (c)
maximum Z-displacement for the 100 individual simulations.
properly account for the geomechanical risks of GCS, stochastic description of rock properties
is necessary.
3.4 Summary
In this chapter, We present the results of a new optimization framework for CO2 storage
under geomechanical risks, including ground surface uplift and rock failure that can lead
to induced seismicity. Specifically, optimization strategies are sought for maximizing the
amount of CO2 that is stored in a geologic formation, while minimizing geomechanical risk
measures. Rigorous coupled flow and geomechanical simulation is adopted to quantify the
amount of CO2 storage and the risks associated with ground surface uplift and rock plastic
strain. The CO2 storage and the two geomechanical risk measures constitute the three
objective functions that are included in our multi-objective optimization framework. Several
numerical experiments are conducted to demonstrate the application of the developed multiobjective optimization framework to CO2 storage problem with for different number of wells.
The results indicate that different well configurations are selected with and without including geomechanical risks in the optimization problem. When geomechanical risks are
considered, the well locations are determined to avoid unaccepted reservoir uplift and rock
failure while maximizing the amount of CO2 stored. In this case, the optimization seeks a
trade-off solution (well locations) to minimize the risks while maximizing the CO2 storage.
64

The results further show that the use of flow-only simulation models is not suitable for predicting the induced geomechanical risks even if the fracturing pressure is known because the
feedback from geomechanical coupling varies with time and with distance from the well. In
addition to flow properties of the storage formation, such as permeability and porosity, the
injection limit and security of a CO2 storage reservoir depends also on the geomechanical
properties of the reservoir and the over-burden layers such as, Young’s Modulus, Poisson’s
ratio, Biot’s coefficient, cohesion, and friction angle. The complex coupling between the flow
and geomechanical effects calls for coupled-physics models to be integrated into state-ofthe-art optimization frameworks to determine site development plans to maximize storage
capacity while maintaining the risks at an acceptable level.
The main contribution of the study is the development of a new optimization framework
to determine well configurations that balance the trade-off between storage capacity and geomechanical risk. This was accomplished by integrating two-way coupled physics modeling
(multiphase flow and rock deformation and failure processes) for geologic CO2 storage with
multi-objective optimization. We considered ground uplift and plastic failure as geomechanical risks; however, the optimization framework is general and can be used when alternative
risk measures are adopted. The choice of multi-objective optimization is to determine a range
of solutions (well configurations) that offer different options (depending on the acceptable
risk levels) to resolve the trade-off between maximizing storage performance and minimizing
the geomechanical risks to ensure the geomechanical integrity of the geologic formation.
The significance of the proposed framework is in the ability to consider the non-trivial
impact that multiphysics coupling has on the optimization results. The Pareto-front that
identified in multi-objective optimization allows the user to adjust the risk level depending
on existing local or national regulatory geomechanical safety constraints. Incorporation
of geomechanical risks can also result in optimal well configurations covering both high
and low permeability regions of the aquifer in a way that cannot be predicted by a flowonly approach excluding geomechanics. It is also important to note that the optimal well
65

configuration depends on the distribution of flow and mechanical properties of the formation,
which control the dynamic response of the storage formation and the surrounding rock.
66

Chapter 4
Well-controlled Optimization and CO2 Leakage
In this section, we present the results from our numerical experiments based on the proposed optimization framework and the coupled flow-geomechanics-fracture simulation model
described above. The simulation model consists of three injection wells and the locations
of the wells are taken from the optimal results of well placement in section 3.1, which are
shown in Figure 3.2 (b). The total injection duration is 15 years. The optimization results
based on the flow-only model is discussed in section 4.1, while the optimization results for
the coupled flow-geomechanics-fracturing model are presented in section 4.2. In both cases,
the optimization is performed assuming constant injection rates for each well (3 control
variables, u = (uw1
, uw2
, uw3
)). This low-dimensional problem allows us to perform linear
interpolations under the total injection volume constraint (Eqn. 2.24) to construct a response
surface (solution space) for illustration purpose.
4.1 Case 4.1: Optimization with Flow-only Simulation
In this section, we use flow-only simulation and minimize the free CO2 as described in
section 2.2.2. The simulation model has the same settings as the model described in section 2.4.3 except that the geomechanical simulation section in CMG-GEM is disabled so
that the geomechanical simulator and the fracturing model is deactivated. To construct
the solution space, we obtain about 1800 sets of injection rates u = (uw1
, uw2
, uw3
) so that
67

uw1 + uw2 + uw3 = 583914 m3/day and 0 ≤ u ≤ 583914 m3/day. The injection rate vector
u is discretized with increments of du = 11000 for each well as shown in Figure 4.1 (a).
Figure 4.1 (b) shows the solution space (total free CO2 amount in moles) for the flow-only
simulation model using the set of control variables in Figure 4.1 (a). With the solution space,
the maximum and minimum values of free CO2 moles are 9.11E+10 for u = (0, 583914, 0)
and 7.67E+10 for u = (31857, 10000, 542057). The objective function has a general decreasing trend from u = (0, 583914, 0) and u = (583914, 0, 0) towards u = (31857, 10000, 542057),
and then the objective value increases about 0.6% towards u = (0, 0, 583914).
(a) Control variables discretization (b) Solution space
Figure 4.1: (a) Distribution of well controls u discretized using du = 11000 representing the
solution space. (b) Total free CO2 (moles) for flow-only simulation using the controls in (a).
Figure 4.2 shows the control optimization results with flow-only simulation using Initialization 1 as described in Table 4.1. Figure 4.2 (a) shows the objective function value vs.
number of optimization iterations. The optimization converges in five iterations to the value
of 7.88E+10 moles of free CO2 amount. Figure 4.2 (b) shows the 3-D plot of gas plume
distribution for the optimized well injection rates u = (93748, 59052, 431113). The optimization solution allocates 70% of the total injection to Well 3 to minimize the total free gas in
the aquifer. Since the permeability around Well 3 is the highest among the three injectors,
the injected CO2 has a wider spread to create more contact area with the resident brine.
More gas can then be dissolved and trapped by injecting into Well 3, thereby minimizing the
68

free CO2. The optimal well controls are subsequently employed in a coupled-physics simulation involving a caprock fracturing model. In Figure 4.2 (c), we could observe the resulting
fracture permeability within the caprock layers. The red cells indicate fracture openings,
with fracture permeability increased from 0 to 100 mD. This indicates when conducting
flow-only simulations, the optimization process converges to a solution that is unfavorable
from a geomechanical perspective. In the next section, we perform well control optimization
that incorporates coupled flow-geomechanics-fracturing and CO2 leakage to quantify and
minimize the geomechanical risk.
Table 4.1: Well flow rates for different initializations and the corresponding local solutions.
Initial well controls Initial Optimal well controls Local Solution
Init. u =

u
0
w1
, u0
w2
, u0
w3

f(u
0
,m) u =

u
⋆
w1
, u⋆
w2
, u⋆
w3

f(u
⋆
,m)
(m3/day) Coupled Flow (m3/day) Coupled Flow
Coupled (moles) Coupled (moles)
1 (194638, 194638, 194638) 2.43E+13 8.41E+10 (84747, 208606, 290560) -19446.54 7.88E+10
2 (58391, 58391, 467131) 6.79E+14 - (187135, 32944, 363834) -19952.56 -
(a) (b) (c)
Figure 4.2: (a) Optimization objective function, (b) 3-D gas plume, and (c) the injection
rates corresponding to optimization solution with the flow-only simulation and Initialization
2.

4.2 Case 4.2: Deterministic Optimization
In this section, we present the results from our numerical experiments based on the proposed
optimization framework with the advanced numerical optimization algorithm (Eqn. 2.24)
and the coupled flow-geomechanics-fracture simulation model described in section 2.4.3.
We first construct the solution space by obtaining about 1800 sets of injection rates u =
(uw1
, uw2
, uw3
) so that uw1 + uw2 + uw3 = 583914 m3/day and 0 ≤ u ≤ 583914 m3/day. The
injection rate vector u is again discretized with increments of du = 11000 for each well as
shown in Figure 4.1 (a).
Figure 4.3 shows the solution space using injection controls in Figure 4.1 (a). The xand y-axes are the injection controls for uw1 and uw2, respectively. The corresponding value
of uw3 can be computed using the linear constraint. The z-axis is plotting the log10 of the
objective function values. The smooth surface plot is generated using linear interpolation. In
Figure 4.3, the color scheme of the ‘plateau’ region appears differently from the rest of the region (as indicated by the use of a different colorbar on the right). From the surface plot of the
total CO2 leakage amount, we observe that the solution space is non-convex, and there exist
many local minima. We also observe that the simulation case with the highest objective value
happens when Well 1 has injection volume equals the injection constraint, u=(583914,0,0).
In this case, the objective function value is 2.57E+20. When Well 2 has injection volume
equals to the injection constraint, u=(0,583914,0), the objective function value is 1.61E+20.
For Well 3 as the main injection well, u=(0,0,583914), and the objective function value is
2.11E+17. For evenly distributed injection controls, u=(194638,194638,194638), and the
objective function value is 2.43E+13. The minimum point on the response surface plot has
an objective function value of -20156.16 with controls u=(143973,30000,409940). For this
minimum point, no caprock fracturing or CO2 leakage occurs, and the stress state is far away
from the fracture opening stress.
Table 4.1 summarizes the initial well controls, initial objective value, the optimal well
controls, and the local solutions for two initializations of the well control optimization. For
70

the first initialization, we use an average injection rate of 194638 m3/day, distributed evenly
across all wells, while ensuring that the total injection does not exceed the constraint of
583914 m3/day for the sum of the injection rates of all three wells. We use an evenly
distributed injection rate for initialization because this approach facilitates a more uniform
pressure dissipation across the reservoir. For initialization 2, we allocate most of the injection
rate to well 3, because it is located in a higher permeability region which can facilitate more
effective pressure dissipation.
Figure 4.4 shows the optimization results of the objective function vs. iteration number for initialization 1 (blue) and initialization 2 (black) in log10 scale. The optimization
algorithm successfully reduces the magnitude of the objective function from 1013 − 1014 to
10−4
, a reduction of approximately 17 to 18 orders of magnitude. Different initializations
lead to distinct local solutions, summarized in Table 4.1. Figure 4.3 depicts the optimization
trajectories of the two initializations, revealing that both of the resulting solutions converge
towards the ’plateau’ region of the solution space. In the case of initialization 1, the resulting solution is trapped in a local minimum, as demonstrated in the subplot in the upper
right corner of Figure 4.3. To prevent such occurrences, it is advisable to conduct multiple
optimizations using diverse initializations.
Figure 4.5 compares the initial and final solutions for the gas plume and injection schedules for the two initializations. For Initialization 1, the injection rate for Well 1 is brought
down from 194638 m3/day to 84747.16 m3/day, the injection rate for Well 2 is almost unchanged, and the injection rate for Well 3 has increased from 194638 m3/day to 290560.56
m3/day. The CO2 leakage has been minimized from 4.93E+06 moles to 0 moles. From
Figure 4.5, we can observe that the caprock fracturing in the caprock layer in plot (a) is
reduced to zero in plot (b). This is because the high injection rate in Well 1 result in high
pressure build-up at the initialization, which can trigger fracturing and potential leakage
in the caprock layer above Well 1. In Figure 4.5 (b), since Well 3 gets a higher injection
volume, pressure build-up is no longer concentrated around Well 1. Hence, pressure can be
71

Local minimum
x 105
Figure 4.3: Solution space and optimization trajectory in log10 scale. The blue square, blue
curve and cross, and blue star are the initial solution, optimization trajectory, and optimal
solution for initialization 1, respectively. The black square, black curve and cross, and black
star are the initial solution, optimization trajectory, and optimal solution for initialization 2,
respectively. The subplot in the upper right corner illustrates an instance of a local minimum
located within the red circle on the response surface.
dissipated more evenly in the reservoir to allow for a lower average pressure build-up, resulting in a smaller stress change in the reservoir during injection to prevent rock fracturing.
The same situation can be observed in the optimization results with the other initialization. In optimization results with Initialization 2 (Figure 4.5 (d)), the caprock layers above
Well 3 are fractured. In the corresponding local solution, caprock fracturing is prevented by
distributing a portion of the injection volume from Well 3 to Well 1. This observation indicates that the injection allocation distributes the pressure in the reservoir according to the
72

Figure 4.4: Optimization results of objective function vs. iteration number for initialization
1 (blue) and initialization 2 (black) in log10 scale. Solid lines and filled circles are plotted
against the left y-axis, while dashed lines and unfilled circles are plotted against the right
y-axis.
heterogenous flow properties to manage and minimize caprock fracturing and the total CO2
leakage amount. Figure 4.6 is showing plots of the initial (solid curves) and optimal (dashed
curves) pressure (black curves) and minimum effective stress (blue curves) over time for two
gridblock locations: (a) a location at well 1 with index={22, 8, 27}, and (b) a location at
well 3 with index={9, 18, 27}. Two observations can be found from the plots. First, we find
that the initial pressure and minimum stress are consistently higher and lower, respectively,
compared to their optimal values in both cases examined. Second, we observe non-monotonic
trends in pressure and minimum stress for the initial cases, with caprock fracturing occurring
at approximately year 14 for both. No caprock fracturing is observed in the optimal cases,
which exhibit monotonic trends in pressure and minimum stress. These results suggest that
optimization effectively reduces pressure buildup, leading to improved pressure and stress
dissipation, ultimately preventing caprock fracturing.
73

w1 w2
w3
w1 w2
w3
w1 w2
w3
w1 w2
w3
Figure 4.5: Optimization results for the two initializations. Initial gas plume (plot a, d) and
injection rate (plot c, f) are compared with the optimal (local) gas plume (plot b, e) and
injection rate (plot c, f) for the two cases in year 15. The blue gridblocks are caprock regions
with fractures, visualized as regions with high fracture permeability in the k-direction. The
red blocks underneath represent the gas plume in the injection reservoir.
0 5 10 15
Year
1
1.5
2
2.5
Pressure (kPa)
104
5000
6000
7000
8000
9000
10000
min effective stress
(a) At well 1
0 5 10 15
Year
1
1.5
2
2.5
Pressure (kPa)
104
5000
6000
7000
8000
9000
10000
min effective stress
(b) At well 3
initial pressure
optimal pressure
initial effective min stress
optimal effective min stress Figure 4.6: Plots of the initial (solid curves) and optimal (dashed curves) pressure (black
curves) and minimum effective stress (blue curves) over time for two gridblock locations: (a) a
location at well 1 with index={22, 8, 27}, and (b) a location at well 3 with index={9, 18, 27}.
74

4.3 Case 4.3: Robust Optimization under Geological
Uncertainties
The robust optimization, which is based on the proposed workflow (section 2.2.3) and the
robust optimization formulation (Eqn. 2.27), is performed using the same model described
in section 2.4.3. We generated 100 realizations of 3-D permeability fields conditioned on the
three wells presented in Figure 4.7 (a) as the hard data, the porosity map in Figure 4.7 (b)
as the soft data, and the permeability distribution in Figure 4.7 (c). Figure 4.8 shows
the mean and standard deviation permeability fields over the ensemble of 100 realizations.
Figure 4.8 (c) shows the permeability of the first five realizations of the ensemble. The total
injection duration is 15 years, and the optimization is performed assuming constant injection
rates for each well. This gives us three control variables, u = (uw1, uw2, uw3). In this case,
100 simulation runs are conducted with 300 gradient evaluations per optimization iteration.
The set of initial controls of initialization 1 in Table 4.1 is used as the initialization of the
robust optimization.
(a) Well locations on
log10(permeability) map
at injection layer
w1 w2
w3
(b) Well locations on porosity map
w1 w2
w3
(c) Histogram of permeability
Figure 4.7: (a) Well locations on log10(permeability) map at the injection layer. (b) well
locations on the porosity map. (c) histogram of the log10 (permeability) map. (a) Well
locations on log10(permeability) map at the injection layer. (b) well locations on the porosity
map. (c) histogram of the log10 permeability map.
We first perform an exploratory search using about 50 control schedules to construct the
response surface illustrated in Figure 12 4.9. Figure 4.9 (a) displays the robust optimization
solution space in log10 scale. The smooth surface shown in the figure was generated through
75

(c) log10(permeability) maps for the first 5 realizations
Figure 4.8: Permeability maps of ensemble mean (a) and standard deviation (b) over 100
realizations. (c) log10(permeability) maps for the first 5 realizations.
linear interpolation of the filled ∼50 circles. Figure 4.9 (b) presents a zoomed-in view of
the plateau region in (a), where no caprock fracturing occurred over the entire ensemble
with the specific control schedules. We observed that the plateau region occurred within the
approximate ranges of uw1 = [2.6E + 04, 9E + 04], uw2 = [1.4E + 05, 2E + 05], and uw3 =
[3.35E + 05, 3.7E + 05], with the injection rate allocation for well 3 being higher than that of
wells 1 and 2. This observation aligns with the findings from the deterministic optimization
results, which indicate that the non-caprock-fracturing ‘plateau’ region emerges when well
3 exhibits a higher injection rate compared to wells 1 and 2. Moreover, this observation
supports the conclusion that the injection allocation strategy effectively distributes pressure
throughout the reservoir according to its heterogeneous flow properties, thereby controlling
and reducing caprock fracturing and the overall CO2 leakage amount.
The robust optimization results and trajectory are shown in Figure 4.10 and Figure 4.9 (ab), respectively. The objective function value for robust optimization decreases from 2.97E+19
76

with the initial injection rates u
0 = (194638, 194638, 194638) to -33933.92 with the optimal
injection rates u
⋆ = (31385.86, 183236.88, 369291.26). This represents a minimization of approximately 23 orders of magnitude. Figure 4.11 presents three distinct result plots derived
from the robust optimization. Plot (a) is a box plot illustrating the ensemble of objective
function values at each iteration, with black boxes plotted against the left y-axis and blue
boxes plotted against the right y-axis. This plot reveals that the ensemble results at iteration 0 are clustered at high objective values, with orders of magnitude ranging from 13 to
21, indicating that caprock fracturing occurs for all realizations. However, from iteration 2
to iteration 10, the ensemble results are consistently clustered at negative values, signifying
that no caprock fracturing takes place in any realization. As the optimization advances, the
maximum value, third quartile, median, first quartile, and minimum value all exhibit a decreasing trend. Plot (b) features a cumulative density function (CDF) plot for the ensemble
of objective function values (in log10 scale) for all iterations, while plot (c) offers a magnified
view of plot (b), concentrating on the x-axis range between -4.56 and -4.515. Table 4.2 shows
the 95th percentile on CDF in Figure 4.11 for each optimization iteration. Through the CDF
plots and the table, we observe that the ensemble of objective values exhibits a decreasing
trend as the optimization progresses.
Table 4.2: True value of the 95th percentile on CDF in Figure 4.11 for each optimization
iteration.
Iteration # 0 1 2 3 4 5 6 7 8
F
−1(0.95) 1.32E+20 -31702.97 -32277.51 -32299.81 -32353.41 -32357.13 -32383.22 -32386.95 -32389.18
Figure 4.12 showcases three plots, each displaying different variables against time at the
grid cell of well 3 with index={9, 18, 27}: (a) pressure, (b) minimum effective stress, and (c)
cumulative CO2 leakage amount. In all three plots, the grey and cyan lines represent the
ensemble values for initial and optimal results, respectively, while the black and blue lines
illustrate the ensemble average values for initial and optimal results, respectively. Plot (a)
shows that the pressure vs. time ensemble for the optimal control schedules has a significantly
better pressure dissipation than the initial case, as indicated by the much higher values
77

(a) Robust optimization solution space (b) Plateau region
Figure 4.9: (a) Robust optimization solution space and optimization trajectory in log10 scale.
Linear interpolation is used to generate the smooth surface from the ∼50 filled circles. The
black square, black curve and cross, and black star are the initial solution, optimization
trajectory, and optimal solution for initialization 1, respectively. (b) A zoomed-in view for
the ‘plateau’ region.
(c) Optimization Results
Figure 4.10: Robust optimization results shown for initialization 1. The objective function
in log10 scale vs. the optimization iteration number. The black solid lines and filled circles
are plotted against the left y-axis, while the blue dashed lines and unfilled circles are plotted
against the right y-axis.
for the latter. This suggests that the optimal control schedules are effective in reducing
average pressure buildup. Conversely, plot (b) displays lower values of minimum effective
78

(a) Box plot of ensemble vs. iteration # (b) CDF plot of ensemble (c) Zoomed-in CDF plot of ensemble
Figure 4.11: (a) a box plot depicting the ensemble of objective function values at each
iteration, with black boxes plotted against the left y-axis and blue boxes plotted against
the right y-axis; (b) a CDF plot illustrating the ensemble objective function values for each
iteration. The light color box has a modified x-axis for the ease of plotting; and (c) a
magnified view of plot (b) focusing on the x-axis range between -4.56 and -4.515.
stresses vs. time for the initial ensemble, while the ensemble of optimal results shows higher
values. This trend aligns with the pressure behavior according to Eqn. 2.3 and suggests
that the minimum effective stress for the optimal ensemble is far from the fracture opening
stress and the caprock has not undergone fracturing. The non-monotonic trends observed
in the initial ensembles of pressure and minimum effective stress are consistent with our
previous observations of fracturing. Finally, plot (c) reveals a correspondence between total
CO2 leakage, pressure, and minimum effective stress. The optimal ensemble shows zero
leakage, whereas the initial ensemble shows leakage accumulation, further emphasizing the
effectiveness of the optimal control schedules.
Figure 4.13 presents the results of the robust optimization for initialization 1. The average
initial gas plume (a) and injection rate (c) are compared with the average optimal (local) gas
plume (b) and injection rate (c) over the ensemble at year 15. The blue grid blocks represent
the caprock fractures, and the lower red gridblocks represent the gas plume in the injection
reservoir. The plots demonstrate that caprock fracturing is suppressed in the optimal case
because some of the injection volume from well 1 is redirected into well 3, while the injection
volume in well 2 remains relatively unchanged.
79

0 5 10 15
Year
1
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
Pressure (kPa)
104 (a) Pressure (kPa)
initial ensemble
initial average over ensemble
optimal ensemble
optimal average over ensemble
0 5 10 15
Year
4000
5000
6000
7000
8000
9000
10000
Minimum Stress (kPa)
(b) Minimum Stress (kPa)
0 5 10 15
Year
0
0.5
1
1.5
2
2.5
Total leakage (mole)
1010 (c) Total leakage (mole)
Figure 4.12: This figure showcases three plots, each displaying different variables against
time at the grid of well 3 with index={9, 18, 27}: (a) pressure, (b) minimum effective stress,
and (c) cumulative CO2 leakage amount. In all three plots, the grey and cyan lines represent
the ensemble values for initial and optimal cases, respectively, while the black and blue lines
illustrate the ensemble average values for initial and optimal cases, respectively.
w1 w2
w3
w1 w2
w3
Figure 4.13: Robust optimization results for initialization 1. Averaged initial gas plume
(a) and injection rate (c) is compared with the averaged optimal (local) gas plume (b) and
injection rate (c) over ensemble at year 15. The blue blocks represent caprock fractures in
terms of their fracture permeability in the k-direction, and the lower red blocks represent
the gas plume in the injection reservoir.
4.4 Summary
In this chapter, we introduces a general optimization framework to minimize CO2 leakage
risk due to caprock fracturing by using a coupled flow-geomechanics-fracturing modeling approach under geological uncertainty. A new CO2 leakage metric is formulated to incorporate
an optimization penalty function together with a term used to quantify failure potential as
the difference between the minimum effective stress and the fracture opening stress. Our
findings highlight the limitations of flow-only simulation models for optimizing storage under
80

leakage risk, emphasizing the importance of coupled flow and geomechanics simulation for
accurately modeling CO2 leakage potential and addressing caprock integrity and CO2 storage optimization challenges. The proposed optimization algorithm determines the optimal
well injection schedule (flow rate vs. time profile) for a fixed total injection amount. The
optimal schedule distributes pressure within the storage formation with heterogeneous flow
properties to minimize the caprock fracturing potential and the total CO2 leakage. The optimization process continually seeks improved solutions, maximizing the difference between
the minimum effective stress and the fracture opening stress to reduce fracturing potential
further and to maximize CO2 storage capacity. The paper bridges the gap between the
CO2 caprock integrity problem and the CO2 storage optimization problem by integrating a
new CO2 leakage metric into a leakage minimization problem to enhance the CO2 storage
efficiency.
In this work, the objective function was formulated under a total injection constraint.
As a result, specifying an excessive total CO2 injection can lead to no solution, that is,
with an overly excessive injection amount, all strategies will result in caprock failure and
leakage. In such cases, one can repeat the optimization by reducing the total CO2 injection
to find solution with minimum caprock failure risk. Alternatively, one can employ other
optimization techniques, including multi-objective, sequential, two-step, or weighted sum
formulations, to maximize total injection while minimizing CO2 leakage amount. Moreover,
the risk metric is a used defined function and can be formulated in different ways. For
example, if fracturing is viewed as a discrete or extreme event, statistical measures other
than expectation might be more appropriate for characterization and quantification of risk.
As an alternative, in robust optimization, the confidence level can be used as the statistical
measure [43], for instance, by ensuring only less than 5% of the total number of realizations
experience fracturing. This implies that optimization results could have leakage occur within
very few scenarios using a pre-defined confidence level. It is important to recognize that such
risks can originate from various sources. In this study, we only considered the uncertainty
81

in describing the heterogeneity of the storage formation as one of the common sources of
uncertainty.
An important consideration in implementing the proposed framework is the high computational cost associated with robust optimization, involving a large set of geological realizations and the use of computationally demanding coupled flow and geomechanics simulation. To solve a large-scale optimization problem where space and time complexities can
become the limiting factor to prevent converging within a realistic time frame for industrial application, stochastic gradient descent [11] [12], mini-batch gradient descent [34], and
semi-stochastic gradient descent [87] algorithms are typically applied in the field of Machine
Learning. To further save the computational cost, the reduced sampling method [139] can
be introduced where a much smaller subset of the model realizations is randomly sampled to
represent the full ensemble and thus reducing the total number of simulation runs. Another
alternative is the use of proxy models in place of full-scale coupled flow and geomechanics
simulation, as will be shown in the next Chapter. Given the emerging trend in the use
of coupled physics simulation, the computational considerations associated with modeling
workflows that require many simulation runs presents an important challenge. Advances
in proxy model development and machine learning techniques can be leveraged to address
the computational demand of intensive coupled physics simulators, by introducing different
levels of approximation error.
82

Chapter 5
Surrogate Modeling of CO2 Coupled-physics
Simulation
5.1 Multi-model Learning
5.1.1 U-Net with permeability
In this section, we present the results of Multi-model learning (section 2.6.1) and show
the performance of the proposed U-Net model with permeability concatenated with well
location information as the input data. A total of 200 unseen testing data is used to test
the performance of the four U-Net models. The testing MSEs are summarized in Table D.1.
Figure 5.1 shows the plots of the results for pressure prediction using the U-Net model
for 11 selected cases. In the first row, the maps display the locations of wells denoted by
green and black dots on the permeability map for each of the 11 cases. In the second
row, we present the pressure predictions at the injection layer. The first plot within the
dashed box identifies the region of heterogeneous permeability. The third row illustrates the
actual pressure observations at the same layer. Lastly, the fourth row depicts the mismatch
between the predicted pressure and the true pressure observations. In the fourth row, it
is evident that the proposed U-Net model demonstrates superior predictions, exhibiting
minimal spatial mismatch values. A similar trend emerges when applying the proposed
83

U-Net model to vertical displacement, as depicted in Figure 5.2 for vertical displacement
prediction. It’s worth noting that the mismatch in the case of vertical displacement is notably
smaller than that observed for pressure, primarily due to the smoother nature of the vertical
displacement field. Figure 5.3 provides a visual representation of the prediction results for
CO2 saturation. The comparison between predicted and observed values generally reveals a
good match to capture both horizontal and vertical spatial correlations. Although a relatively
higher mismatch is noticeable, particularly in the boundary cells. This discrepancy can be
attributed to the inherently less smooth [98] [20] and sparser nature of the CO2 saturation
field.
Figure 5.1: U-Net model prediction for pressure (with permeability as one of the input
data). The black vertical lines indicate wells on the maps. The first row plots the wells on
permeability map for each of the 11 cases. The second row shows the 3-D pressure prediction.
The third row shows the true observation of pressure at the same layer, and the fourth row
shows the difference between predicted pressure and the true observation of pressure.
We perform multiple scenarios for the training of effective plastic strain with different
input structure and loss functions. The testing MSE errors are summarized in Table 5.1.
Figure 5.4 provides a visual representation of the prediction results for effective plastic strain
across multiple training scenarios in Table 5.1. Both the table and the plots clearly indicate
that the testing error is minimized when using predicted pressure (case 4) as the training
input for the same training hyperparameter settings. We observe that when Mean Squared
84

Figure 5.2: U-Net model prediction for vertical displacement, Dz (with permeability as one
of the input data). The black vertical lines indicate wells on the maps. The first row plots
the wells on permeability map for each of the 11 cases. The second row shows the 3-D Dz
prediction. The third row shows the true observation of Dz at the same layer, and the fourth
row shows the difference between predicted Dz and the true observation of Dz.
Figure 5.3: U-Net model prediction for CO2 Saturation (with permeability as one of the input
data). The black vertical lines indicate wells on the maps. The first row plots the wells on
permeability map for each of the 11 cases. The second row shows the 3-D CO2 Saturation
prediction. The third row shows the true observation of CO2 Saturation at the same layer,
and the fourth row shows the difference between predicted and the true observation of CO2
Saturation.
Error (MSE) is used as the loss function, the sparsity of the output data can lead to a very
85

small MSE value during the initial training iterations, causing the training to converge too
rapidly without sufficiently capturing the input-output relationship.
Additionally, it is worth noting that incorporating the permeability map alongside well
location data as input results in deteriorated training performance. This phenomenon can
be attributed to the fact that the output variable, effective plastic strain, exhibits strong
localization around the well locations. When permeability information is introduced as part
of the input, it tends to overshadow the significance of well locations. This observation
becomes evident when analyzing the feature maps of the convolutional layers, as depicted
in Figure 5.5. In Figure 5.5, we focus on the feature maps within convolutional block no.
3, using a model input comprising well location information concatenated with permeability
maps. The first row displays the 16 feature map layers for well configuration 1, with the well
locations depicted in the final column. The second row illustrates the difference between
feature maps at the same locations between well configuration 1 and 0 (not shown here).
The third row showcases the feature maps at the same locations for well configuration 2, and
the fourth row illustrates the difference between feature maps at the same locations between
well configuration 1 and 2.
From the plots, it becomes evident that the inclusion of permeability data in the input
leads to a near disappearance of the well location information in the feature maps, which
is only discernible when subtracting the permeability information. In contrast, the feature
maps that have undergone permeability subtraction, as shown in the second and fourth rows,
successfully preserve the original well locations and exhibit noticeable distinctions among different well configurations. However, when permeability information remains uncompensated,
as depicted in the first and third rows, the localized well information becomes indistinguishable, and the feature maps display minimal differences among various well configurations.
Furthermore, an analysis of the colorbar reveals that the well location information represents
only 0.2% of the permeability information. This underscores the ease with which well location information can be overshadowed when combined with permeability data. Considering
86

the localized and sparse nature of the effective plastic strain output variable, it is imperative to emphasize local information in the input to enable the model to effectively capture
the input-output relationship during training. With these observations and conclusion, we
conducted another study in which we assign greater importance to the well location input.
Table 5.1 case 2 displays the testing MSE, and it clearly demonstrates that increasing the
weight on well location information leads to improved training performance than that of case
1 and case 3. The weight factor must be determined through a sensitivity analysis.
Table 5.1: Testing errors across various training scenarios for effective plastic strain with the
same hyperparameter settings.
Case Input 1 Weight for Input 1 Input 2 Weight for Input 2 Testing MSE Error
1 well location 1 permeability map 1 1.51E-05
2 well location 1E4 permeability map 1 1.709E-06
3 well location 1 - - 9.384E-06
4 Predicted Pressure 1 - - 8.810E-07
87

Well locations on
permeability map
True
Observation
Case 1
Prediction
Case 1
Difference
Case 2
Difference
Case 3
Difference
Case 4
Difference
Case 2
Prediction
Case 3
Prediction
Case 4
Prediction
Figure 5.4: Different cases on Table 5.1 of U-Net model prediction for effective plastic strain.
The green and black dots are well locations shown on the map. The first row shows the well
locations on permeability maps. The second row plots the true observation for each of the
11 cases. The rest of the rows show the prediction and mismatch at the first layer for each
case.
88

Feature maps at convolutional block no.3
Well config 1 Channel 1
layer 1
Channel 1
layer 16
Difference between Well config 1 and well config 0
Well config 2
Difference between Well config 2 and well config 1
Well
locations
Figure 5.5: Feature maps at convolutional block no.3 with model input of well location
information concatenated with permeability maps. The first row plots the 16 layers on
channel 1 for well configuration 1 with well locations plotted on the last column. The
second row shows the difference between the feature maps at the same location between well
configuration 1 and 0 (which is not shown here). The third row shows the feature maps at
the same location for well configuration 2 and the fourth row shows the difference between
the feature maps at the same location between well configuration 1 and 2.
89

5.1.2 U-Net without permeability
When training the model for prediction of non-sparse and smoother parameters, such as
pressure, vertical displacement, and CO2 saturation, the absence of permeability results
in degraded training performance compared to when permeability is included, under the
same model training settings. As demonstrated by the discrepancies in vertical displacement (Figure 5.6) and CO2 saturation predictions (Figure 5.7), some testing cases exhibit
significant underestimation, indicating a substantial gap between predicted values and true
observations. This underscores the importance of incorporating the permeability field into
the model input, as it enables the model to learn more quickly about the input-output
relationship with a better representation of global spatial information.
x
y
Figure 5.6: U-Net model prediction for vertical displacement, Dz (WITHOUT permeability
as one of the input data). The green and black dots are well locations shown on the map.
The first row plots the wells on permeability map for each of the 11 cases. The second row
shows the Dz prediction at the injection layer. The heterogeneous permeability region is
labelled at the first plot inside the dashed box. The third row shows the true observation of
Dz at the same layer, and the fourth row shows the difference between predicted Dz and the
true observation of Dz.
90

x
y
Figure 5.7: U-Net model prediction for CO2 Saturation (WITHOUT permeability as one of
the input data). The green and black dots are well locations shown on the map. The first
row plots the wells on permeability map for each of the 11 cases. The second row shows the
CO2 Saturation prediction at the injection layer. The heterogeneous permeability region is
labelled at the first plot inside the dashed box. The third row shows the true observation of
CO2 Saturation at the same layer, and the fourth row shows the difference between predicted
and the true observation of CO2 Saturation.
91

5.2 Multi-output Learning
In this section, we present a comprehensive analysis of the results obtained from the proposed
multi-output learning approach (Section 2.6.2). We systematically explore four distinct
cases characterized by variations in input data weights, with the primary goal of finding an
optimal balance between the global effects represented by the permeability input and the
localized effects stemming from the well location input. Across these four cases, we maintain
a consistent weight of 1 for the permeability input, while we progressively vary the weights
assigned to the well location input, setting them at 1, 10, 100, and 10,000, respectively. The
training results of the four cases and the MSE losses are summarized in Table 5.2. To enhance
the visual understanding of our findings, we also present 3-D representations of the training
outputs including pressure, vertical displacement, CO2 saturation, and plastic strain. These
visual representations are depicted in Figures 5.8, 5.9, 5.10, and 5.11, respectively.
In Figure 5.8, we present a visual representation of the 3-D pressure field mismatch for
each of the four cases. A notable observation from the figure is that the degree of mismatch
becomes more pronounced as the weight factor associated with well location input increases.
This observation aligns with the increasing trend of the MSE loss as the weight factor of
well location increases for pressure, as depicted in Table 5.2. Figure 5.9 shows the 3-D
vertical displacement field mismatch for the four cases. A similar behavior can be observed
that the degree of mismatch becomes more evident as the weight factor of well location
input increases. This observation also aligns with the increasing trend of the MSE loss for
vertical displacement in Table 5.2. A similar observation can be made when examining
the 3-D CO2 saturation field mismatch in Figure 5.10. An increase in mismatch can be
observed as the weight factor associated with well location input rises. This trend mirrors
the increasing MSE loss of CO2 saturation, as also shown in Table 5.2. Figure 5.11 presents
the prediction of the 3-D effective plastic strain field across the four cases. Unlike the other
three outputs, we observe a contrasting trend where the degree of mismatch decreases as
the weight factor assigned to the well location input increases. When the weight factor for
92

well location input is set to 1, the model fails to generate meaningful predictions for plastic
strain, as evidenced by the empty plots in Figure 5.11, row 3. This pattern corresponds to
the decreasing trend in MSE loss as the weight factor associated with well location input
increases. Based on these observations, we can draw the conclusion that outputs such as
pressure, vertical displacement, and CO2 saturation exhibit a significant global influence
and rely more on the permeability input. In contrast, the localized and sparse output, such
as effective plastic strain, relies more on well location information as input. To strike an
optimal balance between these influences, a well location weight factor of 10 appears to be
a promising choice. It’s worth noting that this weight factor can be further fine-tuned and
optimized through optimization techniques.
Table 5.2: MSE losses of four training cases with different weight factor of well location
input.
Output MSE Loss
Pressure �!
CO2
Saturation
Plastic
strain
Permeability
weight
Well location
weight Case
1.87E-05 3.44E-04 3.04E-05 1.09E-05
1 1 1
bad good good good
4.50E-06 3.63E-04 3.03E-05 1.04E-05
2 1.00E+01 1
good good good good
6.58E-06 3.96E-04 3.78E-05 1.41E-05
3 1.00E+02 1 underestimate
underestimate good good
4.84E-06 6.39E-04 1.30E-04 6.11E-05
4 1.00E+04 1 underestimate
underestimate
underestimate good 93

Wells on permeability map
True Observation (Pressure)
Mismatch
Wells weight factor = 1
Wells weight factor = 1E+1
Wells weight factor = 1E+2
Wells weight factor = 1E+4
Figure 5.8: Seven examples of the 3-D pressure field mismatch for each of the four cases in
Table 5.2. The first row shows the well locations on permeability map, the second row shows
the true observations, and the third to the sixth rows shows the pressure mismatch of the
four cases.
94

Wells on permeability map
True Observation (�!)
Mismatch
Wells weight factor = 1
Wells weight factor = 1E+1
Wells weight factor = 1E+2
Wells weight factor = 1E+4
Figure 5.9: Seven examples of the 3-D vertical displacement field mismatch for each of the
four cases in Table 5.2. The first row shows the well locations on permeability map, the
second row shows the true observations, and the third to the sixth rows shows the vertical
displacement mismatch of the four cases.
95

Wells on permeability map
True Observation (CO2 saturation)
Mismatch
Wells weight factor = 1
Wells weight factor = 1E+1
Wells weight factor = 1E+2
Wells weight factor = 1E+4
Figure 5.10: Seven examples of the 3-D CO2 saturation field mismatch for each of the four
cases in Table 5.2. The first row shows the well locations on permeability map, the second
row shows the true observations, and the third to the sixth rows shows the CO2 saturation
mismatch of the four cases.
96

Wells on permeability map
True Observation (�!"#
$ )
Prediction (�!"#
$ )
Wells weight factor = 1
Wells weight factor = 1E+1
Wells weight factor = 1E+2
Wells weight factor = 1E+4
Figure 5.11: Seven examples of the 3-D effective plastic strain field mismatch for each of
the four cases in Table 5.2. The first row shows the well locations on permeability map, the
second row shows the true observations, and the third to the sixth rows shows the effective
plastic strain mismatch of the four cases. Row 3 is empty because the model fails to predict
the plastic strain output.
97

5.3 Summary
In this chapter, we presented the results of our study on the performance of the proposed
U-Net model, which incorporates permeability data concatenated with well location information as input, with two distinct workflows of multi-model learning and multi-output learning.
We tested this model using 200 unseen testing data points and evaluated its performance
across multiple parameters, including pressure, vertical displacement, CO2 saturation, and
effective plastic strain. For pressure, vertical displacement, and CO2 saturation predictions,
the inclusion of permeability data in the input significantly improved the model’s training
performance. The training results demonstrated that the proposed U-Net model produced
predictions that closely matched the true observations for both training approaches, especially in regions with heterogeneous permeability. This suggested that the inclusion of
permeability data allowed the model to capture global spatial information more effectively.
In the case of CO2 saturation prediction, while there was generally a good match between
predicted and observed values, some discrepancies were observed, particularly in boundary
cells. This discrepancy could be attributed to the inherent sparsity and less smooth nature of the CO2 saturation field. Nonetheless, the proposed U-Net model demonstrated its
ability to capture both horizontal and vertical spatial correlations in this parameter. We
also explored different training scenarios for effective plastic strain, considering various input structures and loss functions. The results revealed that using predicted pressure as the
training input yielded the lowest testing error. This indicated the importance of selecting
appropriate training inputs to effectively capture the input-output relationship. Notably,
using the Mean Squared Error as the loss function led to rapid convergence during training,
potentially sacrificing the model’s ability to capture complex relationships. Additionally, we
observed that including the permeability map alongside well location data as input could
lead to deteriorated training performance for effective plastic strain prediction. This was
primarily due to the sparse and localized nature of the effective plastic strain field and that
permeability data overshadowing the significance of well locations in the feature maps of
98

the convolutional layers. The importance of preserving local well information in the input
was underscored. In light of these observations and conclusions, we conducted a follow-up
study in which we assigned greater importance to the well location input by adjusting the
weight factor. We observe that increasing the weight on well location information resulted
in improved training performance compared to other cases, highlighting the significance of
carefully determining this weight factor through sensitivity analysis. In summary, our study
demonstrated the effectiveness of the proposed U-Net model in predicting various subsurface
parameters, with the inclusion of permeability data enhancing performance for non-sparse
parameters. The choice of training inputs and the balance between well location and permeability information were found to be crucial in achieving accurate predictions, emphasizing
the need for thoughtful model configuration and sensitivity analysis in subsurface modeling
applications.
This study focuses on the development of a deep-learning-based surrogate model aimed
at rapidly predicting geomechanical risks in various well configuration scenarios. The model
takes as input well location information and a deterministic permeability map, and it predicts
the final time step values of key variables of interest, including pressure, CO2 saturation, vertical displacement, and effective plastic strain within the field. Our research demonstrates the
superior performance of the proposed surrogate models compared to directly mapping well
locations to the state variables of interest. Notably, the predictive accuracy for geomechanical outputs, such as vertical displacement and plastic deformation, is particularly strong.
This is significant as vertical displacement, which predominantly remains in the elastic deformation regime during simulation, exhibits superior predictive results. Furthermore, the
model’s excellent performance in predicting plastic deformation, characterized by its higher
degree of non-linearity according to governing equations, underscores the effectiveness of our
approach.
99

However, it’s essential to clarify the limitations and future directions of this study. Firstly,
our scope does not encompass the prediction of time-series data or addressing geologic uncertainties with robust flow properties, which can pose challenges in problem setup, model
structure, and training methods. These aspects are earmarked for future work. Additionally,
we recognize that parameters associated with CO2 storage, featuring an even higher degree
of non-linearity, such as fault reactivation, caprock fracturing, and CO2 leakage, could lead
to a more intricate input-output relationship. Addressing these complexities in building
surrogate models for CO2 storage will be a future research endeavor. We anticipate that
incorporating physics-informed deep learning techniques will aid in constraining the inputoutput relationship during training, and this will also be explored in our future work. In
conclusion, this study advances the development of deep learning-based surrogate models
for geomechanical risk prediction, demonstrating their effectiveness in various well configuration scenarios. We acknowledge the study’s limitations and lay out a clear roadmap for
future research, including addressing time-series data prediction, geologic uncertainties, and
more complex CO2 storage scenarios, thus contributing to the evolving field of subsurface
modeling.
100

Chapter 6
Conclusion
In conclusion, this research covers three distinct but interconnected topics in the field of
CO2 storage optimization, geomechanics, and deep-learning based surrogate modeling for
CO2 storage.
In the first topic, we introduced a novel optimization framework for CO2 storage under
geomechanical risks. This framework aimed to balance the maximization of CO2 storage capacity with the minimization of geomechanical risks such as ground surface uplift and rock
failure leading to induced seismicity. We integrated coupled flow and geomechanical simulations into a multi-objective optimization framework, considering three objective functions:
CO2 storage, ground uplift risk, and rock plastic strain risk. Our research demonstrated
that when geomechanical risks are taken into account, the optimization process yields well
configurations that maximize CO2 storage while avoiding unacceptable risks. We highlighted
the critical importance of coupling flow and geomechanics for accurate predictions and emphasized that the choice of geological properties and well locations can significantly impact
optimal solutions. This chapter’s novelty lies in the integration of geomechanical risk considerations into CO2 storage optimization, providing a valuable framework for addressing the
trade-off between storage capacity and geomechanical integrity.
The second topic introduced an innovative optimization framework focused on minimizing CO2 leakage risk due to caprock fracturing. This novel approach incorporated coupled flow-geomechanics-fracturing modeling under geological uncertainty. We formulated a
101

unique CO2 leakage metric that combined optimization penalty functions and failure potential, addressing the challenge of caprock integrity in CO2 storage optimization. Our research
highlighted the limitations of relying solely on flow-only simulation models, emphasizing the
necessity of coupled flow and geomechanics simulations for accurate CO2 leakage predictions.
The chapter introduced an optimization algorithm to determine optimal well injection schedules, promoting efficient pressure distribution within the storage formation while minimizing
leakage potential. The novelty of this chapter lies in its bridging of the caprock integrity
problem with CO2 storage optimization, providing a comprehensive framework to enhance
storage efficiency while addressing leakage risks.
In the third topic, we presented a cutting-edge deep-learning-based surrogate model, the
U-Net, for predicting subsurface parameters relevant to CO2 storage. This model incorporated permeability data and well location information as inputs and demonstrated superior
performance, particularly for non-sparse parameters such as pressure, vertical displacement,
and CO2 saturation. The chapter emphasized the importance of well-balanced training
inputs and sensitivity analysis for achieving accurate predictions. Notably, the model’s effectiveness in predicting geomechanical outputs, such as vertical displacement and plastic
deformation, highlighted its potential in geomechanical risk assessment. The novel contribution of this chapter is the development of a deep-learning surrogate model tailored to model
subsurface CO2 storage, offering valuable insights into the role of input configurations in
prediction accuracy and the computational efficiency it brings to the table.
Overall, these three chapters represent significant contributions to the field of subsurface
modeling and CO2 storage optimization. They introduce novel frameworks, optimization
techniques, and deep learning approaches that address the complex challenges associated
with geomechanical risks and caprock integrity. The chapters collectively advance our understanding of CO2 storage optimization and geomechanical risk assessment while setting
the stage for future research directions in these critical areas.
102

Chapter 7
Ongoing and Future Works
In terms of future work, several avenues are identified across the topic of deep learning-based
surrogate modeling. For the surrogate modeling of CO2 storage, there is a need to expand the
networks to address time-series data and geologic uncertainties with robust flow properties.
Furthermore, parameters like fault reactivation, caprock fracturing, and CO2 leakage, which
introduce additional non-linearity and complexities, will be explored.
Figure 7.1 provides a visual representation of the underlying physics governing the deeplearning U-Net neural network CO2 storage surrogate model. As depicted in the plot, several
key relationships exist between various subsurface parameters:
1. Well location information can be related to pressure through the source and sink terms
in the fluid mass-conservation equation.
2. Permeability can be related to pressure primarily through Darcy’s equation.
3. Vertical displacement is intricately related to pressure through the poro-elasticity equation couples flow and geomechanics governing equations by incorporating the effective
stress equation. Additionally, the application of Hooke’s law in the constitutive equation links the stress and strain solutions.
4. CO2 saturation is explicitly extracted from pressure solution, as illustrated by the
poro-elasticity equation, which establishes connections between these parameters.
103

5. Effective plastic strain is linked to displacement and pressure through plasticity theory
and the poro-elasticity equation, emphasizing the interplay between geomechanics and
fluid flow in subsurface modeling.
This comprehensive visual representation highlights the intricate relationships between
these key subsurface parameters, providing valuable insights into the complex physics that
underlie the CO2 storage surrogate model developed in this study.
Future endeavors will be dedicated to fine-tuning the configurations of the deep-learning
model and integrating physics-informed approaches that leverage the governing equations
outlined in Figure 7.1. These equations serve to couple the input-output relationship of
subsurface parameters, and by incorporating them, we aim to enhance the model’s accuracy
and reliability. Addressing the challenges posed by time-dependent geomechanical and flow
properties in more intricate subsurface scenarios will be a key research direction.
In summary, this research contributes valuable insights and methodologies for optimizing
CO2 storage while mitigating geomechanical risks and advancing the field of subsurface
modeling. Future work will build upon these foundations to tackle increasingly complex
challenges in this critical domain.
104

�!"#
$
CO2 Saturation
�!
Constitutive equations
(Hooke’s law):
�̇ = ��%
̇
Effective stress equation
�$
= � − ��,
Effective stress equation
�$
= � − ��,
�̇ = � �̇ − �!̇
Constitutive equation of poroelasticity
Pressure
well locations
+
permeability
Darcy’s equation
�& = − �(!
"
)!
∇�& − �&�
��&
�� + ∇ 4 �& = �&�&,
Physics Behind
Source and sink term in fluid
mass-conservation equation:
��!
�!
=
�
�"#
��$ +
�%
�"#
�� + ��& �! + �''��' + �'()!��()!
Figure 7.1: Physics behind the deep-learning U-Net neural network CO2 storage surrogate
model.
105

References
[1] Ajith Abraham and Lakhmi Jain. Evolutionary multiobjective optimization. Springer,
2005.
[2] Simeon Agada, Sebastian Geiger, Ahmed Elsheikh, and Sergey Oladyshkin. Datadriven surrogates for rapid simulation and optimization of wag injection in fractured
carbonate reservoirs. Petroleum Geoscience, 23(2):270–283, 2017.
[3] Watheq J Al-Mudhafar and Dandina N Rao. Proxy-based metamodeling optimization
of the gas-assisted gravity drainage gagd process in heterogeneous sandstone reservoirs.
In SPE Western Regional Meeting. OnePetro, 2017.
[4] Watheq J Al-Mudhafar, Dandina N Rao, Sanjay Srinivasan, Hung Vo Thanh, and
Erfan M Al Lawe. Rapid evaluation and optimization of carbon dioxide-enhanced
oil recovery using reduced-physics proxy models. Energy Science & Engineering,
10(10):4112–4135, 2022.
[5] Faisal Alenezi and Shahab Mohaghegh. A data-driven smart proxy model for a comprehensive reservoir simulation. In 2016 4th Saudi International Conference on Information Technology (Big Data Analysis)(KACSTIT), pages 1–6. IEEE, 2016.
[6] Steven T Anderson. Cost implications of uncertainty in co 2 storage resource estimates:
A review. Natural Resources Research, 26:137–159, 2017.
[7] Stefan Bachu. Review of co2 storage efficiency in deep saline aquifers. International
Journal of Greenhouse Gas Control, 40:188–202, 2015.
[8] Peyman Bahrami, Farzan Sahari Moghaddam, and Lesley A James. A review of proxy
modeling highlighting applications for reservoir engineering. Energies, 15(14):5247,
2022.
[9] Nick Barton, Strength Bandis, and K Bakhtar. Strength, deformation and conductivity
coupling of rock joints. In International journal of rock mechanics and mining sciences
& geomechanics abstracts, volume 22, pages 121–140. Elsevier, 1985.
[10] Martin Blunt. Carbon dioxide storage. 2010.
[11] L´eon Bottou. Stochastic learning. In Summer School on Machine Learning, pages
146–168. Springer, 2003.
106

[12] L´eon Bottou. Large-scale machine learning with stochastic gradient descent. In
Proceedings of COMPSTAT’2010: 19th International Conference on Computational
StatisticsParis France, August 22-27, 2010 Keynote, Invited and Contributed Papers,
pages 177–186. Springer, 2010.
[13] Royal Harvard Brooks. Hydraulic properties of porous media. Colorado State University, 1965.
[14] Michael S Bruno, Kang Lao, Julia Diessl, Bill Childers, Jing Xiang, Nicky White,
and Ellen van der Veer. Development of improved caprock integrity analysis and risk
assessment techniques. Energy Procedia, 63:4708–4744, 2014.
[15] Abdulrahman Bubshait and Birendra Jha. Revisiting 2013–2014 azle seismicity to
understand the role of barnett production on stress propagation and fault stability.
Geophysics, 87(4):M127–M149, 2022.
[16] McMillan Burton, Navanit Kumar, and Steven L Bryant. Time-dependent injectivity
during co2 storage in aquifers. In SPE Improved Oil Recovery Conference?, pages
SPE–113937. SPE, 2008.
[17] Andreas Busch, Alexandra Amann, Pieter Bertier, M Waschbusch, and Bernhard M
Krooss. The significance of caprock sealing integrity for co2 storage. In SPE International Conference on CO2 Capture, Storage, and Utilization, pages SPE–139588. SPE,
2010.
[18] David A Cameron and Louis J Durlofsky. Optimization of well placement, co2 injection
rates, and brine cycling for geological carbon sequestration. International Journal of
Greenhouse Gas Control, 10:100–112, 2012.
[19] Wenzhuo Cao, Ji-Quan Shi, Sevket Durucan, and Anna Korre. Evaluation of shear
slip stress transfer mechanism for induced microseismicity at in salah co2 storage site.
International Journal of Greenhouse Gas Control, 107:103302, 2021.
[20] Yuan Cao, Zhiying Fang, Yue Wu, Ding-Xuan Zhou, and Quanquan Gu. Towards
understanding the spectral bias of deep learning. arXiv preprint arXiv:1912.01198,
2019.
[21] Fr´ed´eric Cappa and Jonny Rutqvist. Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep underground injection of co2.
International Journal of Greenhouse Gas Control, 5(2):336–346, 2011.
[22] Marco A Cardoso, Louis J Durlofsky, and Pallav Sarma. Development and application
of reduced-order modeling procedures for subsurface flow simulation. International
journal for numerical methods in engineering, 77(9):1322–1350, 2009.
[23] Michael Anthony Celia, Stefan Bachu, Jan Martin Nordbotten, and Karl W Bandilla.
Status of co2 storage in deep saline aquifers with emphasis on modeling approaches
and practical simulations. Water Resources Research, 51(9):6846–6892, 2015.
107

[24] Bailian Chen, Dylan R Harp, Youzuo Lin, Elizabeth H Keating, and Rajesh J Pawar.
Geologic co2 sequestration monitoring design: A machine learning and uncertainty
quantification based approach. Applied energy, 225:332–345, 2018.
[25] Sirui Chen, Shengjie Zhao, and Quan Lan. Residual block based nested u-type architecture for multi-modal brain tumor image segmentation. Frontiers in Neuroscience,
16:832824, 2022.
[26] Judith C Chow, John G Watson, Antonia Herzog, Sally M Benson, George M Hidy,
William D Gunter, Stanley J Penkala, and Curt M White. Separation and capture of
co2 from large stationary sources and sequestration in geological formations. Journal
of the Air & Waste Management Association, 53(10):1172–1182, 2003.
[27] Chieh Chu. Prediction of steamflood performance in heavy oil reservoirs using correlations developed by factorial design method. In SPE California Regional Meeting.
OnePetro, 1990.
[28] Abdullah Cihan, Jens T Birkholzer, and Marco Bianchi. Optimal well placement
and brine extraction for pressure management during co2 sequestration. International
Journal of Greenhouse Gas Control, 42:175–187, 2015.
[29] CMG. Gem user manual, 2019.
[30] Kalyanmoy Deb. Multi-objective optimization. search methodologies. Search
Methodol, 2014:403–449, 2014.
[31] Kalyanmoy Deb and Tushar Goel. Controlled elitist non-dominated sorting genetic
algorithms for better convergence. In International conference on evolutionary multicriterion optimization, pages 67–81. Springer, 2001.
[32] Kalyanmoy Deb, Amrit Pratap, Sameer Agarwal, and TAMT Meyarivan. A fast and
elitist multiobjective genetic algorithm: Nsga-ii. IEEE transactions on evolutionary
computation, 6(2):182–197, 2002.
[33] JP Dejean and G Blanc. Managing uncertainties on production predictions using
integrated statistical methods. In SPE Annual Technical Conference and Exhibition?,
pages SPE–56696. SPE, 1999.
[34] Ofer Dekel, Ran Gilad-Bachrach, Ohad Shamir, and Lin Xiao. Optimal distributed
online prediction using mini-batches. Journal of Machine Learning Research, 13(1),
2012.
[35] Alfhild L Eide, Lars Holden, Edel Reiso, and Sigurd I Aanonsen. Automatic history
matching by use of response surfaces and experimental design. In ECMOR IV-4th
European conference on the mathematics of oil recovery, pages cp–233. European Association of Geoscientists & Engineers, 1994.
[36] Damsleth Elvind, Hage Asmund, and Volden Rolf. Maximum information at minimum
cost: a north sea field development study with an experimental design. Journal of
Petroleum Technology, 44(12):1350–1356, 1992.
108

[37] Bruno Figueiredo, Chin-Fu Tsang, Jonny Rutqvist, Jac Bensabat, and Auli Niemi.
Coupled hydro-mechanical processes and fault reactivation induced by co2 injection
in a three-layer storage formation. International Journal of Greenhouse Gas Control,
39:432–448, 2015.
[38] Jeffrey P Fitts and Catherine A Peters. Caprock fracture dissolution and co2 leakage.
Reviews in Mineralogy and Geochemistry, 77(1):459–479, 2013.
[39] Erling Fjaer, Rune Martin Holt, Per Horsrud, and Arne Marius Raaen. Petroleum
related rock mechanics. Elsevier, 2008.
[40] F Friedmann, A Chawathe, and DK Larue. Assessing uncertainty in channelized reservoirs using experimental designs. SPE Reservoir Evaluation & Engineering, 6(04):264–
274, 2003.
[41] Priya Ravi Ganesh and Srikanta Mishra. Reduced physics modeling of co2 injectivity.
Energy Procedia, 63:3116–3125, 2014.
[42] Sarah E Gasda, Jan M Nordbotten, and Michael A Celia. Application of simplified
models to co2 migration and immobilization in large-scale geological systems. International Journal of Greenhouse Gas Control, 9:72–84, 2012.
[43] Roger Ghanem, Christian Soize, and Charanraj Thimmisetty. Optimal well-placement
using probabilistic learning. Data-Enabled Discovery and Applications, 2:1–16, 2018.
[44] Mohammadreza Ghasemi, Yanfang Yang, Eduardo Gildin, Yalchin Efendiev, and Victor Calo. Fast multiscale reservoir simulations using pod-deim model reduction. In
SPE reservoir simulation symposium. OnePetro, 2015.
[45] Kieran A Gilmore, Chunendra K Sahu, Graham P Benham, Jerome A Neufeld, and
Mike J Bickle. Leakage dynamics of fault zones: experimental and analytical study
with application to co2 storage. Journal of Fluid Mechanics, 931:A31, 2022.
[46] Takashi Goda and Kozo Sato. Global optimization of injection well placement toward
higher safety of co2 geological storage. Energy Procedia, 37:4583–4590, 2013.
[47] David E Goldberg. Genetic algorithms in search, optimization and machine learning
addison welssey publishing company. Reading, MA, 1989.
[48] S Goodarzi, A Settari, MD Zoback, and DW Keith. Optimization of a co2 storage
project based on thermal, geomechanical and induced fracturing effects. Journal of
Petroleum Science and Engineering, 134:49–59, 2015.
[49] Somayeh Goodarzi, Antonin Settari, M Zoback, and DW Keith. Thermal aspects of
geomechanics and induced fracturing in co2 injection with application to co2 sequestration in ohio river valley. In SPE International Conference on CO2 Capture, Storage,
and Utilization. OnePetro, 2010.
109

[50] Dylan R Harp, Rajesh Pawar, J William Carey, and Carl W Gable. Reduced order
models of transient co2 and brine leakage along abandoned wellbores from geologic
carbon sequestration reservoirs. International Journal of Greenhouse Gas Control,
45:150–162, 2016.
[51] CD Hawkes, PJ McLellan, and S Bachu. Geomechanical factors affecting geological
storage of co2 in depleted oil and gas reservoirs. Journal of Canadian Petroleum
Technology, 44(10), 2005.
[52] Jincong He and Louis J Durlofsky. Reduced-order modeling for compositional simulation by use of trajectory piecewise linearization. SPE Journal, 19(05):858–872, 2014.
[53] Jason E Heath, Sean A McKenna, Thomas A Dewers, Jesse D Roach, and Peter H
Kobos. Multiwell co2 injectivity: impact of boundary conditions and brine extraction
on geologic co2 storage efficiency and pressure buildup. Environmental science &
technology, 48(2):1067–1074, 2014.
[54] Dan Hendrycks and Kevin Gimpel. Gaussian error linear units (gelus). arXiv preprint
arXiv:1606.08415, 2016.
[55] Ahmed Khalil Jaber, Sameer Noori Al-Jawad, and Ali K Alhuraishawy. A review of
proxy modeling applications in numerical reservoir simulation. Arabian Journal of
Geosciences, 12(22):701, 2019.
[56] Jayanth Jagalur-Mohan, Birendra Jha, Zheng Wang, Ruben Juanes, and Youssef
Marzouk. Inferring fault frictional and reservoir hydraulic properties from injectioninduced seismicity. Geophysical Research Letters, 45(3):1313–1320, 2018.
[57] Birendra Jha, Francesca Bottazzi, Rafal Wojcik, Martina Coccia, Noah Bechor, Dennis
McLaughlin, Thomas Herring, Bradford H Hager, Stefano Mantica, and Ruben Juanes.
Reservoir characterization in an underground gas storage field using joint inversion of
flow and geodetic data. International Journal for Numerical and Analytical Methods
in Geomechanics, 39(14):1619–1638, 2015.
[58] Birendra Jha and Ruben Juanes. Coupled multiphase flow and poromechanics: A
computational model of pore pressure effects on fault slip and earthquake triggering.
Water Resources Research, 50(5):3776–3808, 2014.
[59] Wei Jia, Brian McPherson, Feng Pan, Zhenxue Dai, Nathan Moodie, and Ting Xiao.
Impact of three-phase relative permeability and hysteresis models on forecasts of storage associated with co2-eor. Water Resources Research, 54(2):1109–1126, 2018.
[60] Wei Jia, Brian J McPherson, Feng Pan, Ting Xiao, and Grant Bromhal. Probabilistic
analysis of co2 storage mechanisms in a co2-eor field using polynomial chaos expansion.
International Journal of Greenhouse Gas Control, 51:218–229, 2016.
[61] Guanghui Jiang, Jianping Zuo, Yulin Li, and Xu Wei. Experimental investigation
on mechanical and acoustic parameters of different depth shale under the effect of
confining pressure. Rock Mechanics and Rock Engineering, 52:4273–4286, 2019.
110

[62] Zhihao Jiang, Pejman Tahmasebi, and Zhiqiang Mao. Deep residual u-net convolution
neural networks with autoregressive strategy for fluid flow predictions in large-scale
geosystems. Advances in Water Resources, 150:103878, 2021.
[63] Zhongyi Jiang, Min Zhu, Dongzhuo Li, Qiuzi Li, Yanhua O Yuan, and Lu Lu. Fouriermionet: Fourier-enhanced multiple-input neural operators for multiphase modeling of
geological carbon sequestration. arXiv preprint arXiv:2303.04778, 2023.
[64] Zhaoyang Larry Jin and Louis J Durlofsky. Reduced-order modeling of co2 storage
operations. International Journal of Greenhouse Gas Control, 68:49–67, 2018.
[65] Zhefei Jin, Weixin Li, Congrui Jin, James Hambleton, and Gianluca Cusatis.
Anisotropic elastic, strength, and fracture properties of marcellus shale. International
Journal of Rock Mechanics and Mining Sciences, 109:124–137, 2018.
[66] Ruben Juanes, EJ Spiteri, FM Orr Jr, and MJ Blunt. Impact of relative permeability
hysteresis on geological co2 storage. Water resources research, 42(12), 2006.
[67] CS Kabir, A Chawathe, SD Jenkins, AJ Olayomi, C Aigbe, and DB Faparusi. Developing new fields using probabilistic reservoir forecasting. SPE Reservoir Evaluation &
Engineering, 7(01):15–23, 2004.
[68] Subhash Kalla and Christopher D White. Efficient design of reservoir simulation
studies for development and optimization. SPE Reservoir Evaluation & Engineering,
10(06):629–637, 2007.
[69] Elizabeth Keating, Diana Bacon, Susan Carroll, Kayyum Mansoor, Yunwei Sun,
Liange Zheng, Dylan Harp, and Zhenxue Dai. Applicability of aquifer impact models
to support decisions at co2 sequestration sites. International Journal of Greenhouse
Gas Control, 52:319–330, 2016.
[70] Joseph Kestin, H Ezzat Khalifa, and Robert J Correia. Tables of the dynamic and
kinematic viscosity of aqueous nacl solutions in the temperature range 20–150 c and the
pressure range 0.1–35 mpa. Journal of physical and chemical reference data, 10(1):71–
88, 1981.
[71] Aaditya Khanal and Md Fahim Shahriar. Physics-based proxy modeling of co2 sequestration in deep saline aquifers. Energies, 15(12):4350, 2022.
[72] Dieter Kraft. A software package for sequential quadratic programming.
Forschungsbericht- Deutsche Forschungs- und Versuchsanstalt fur Luft- und Raumfahrt, 1988.
[73] Ajitabh Kumar, R Ozah, M Noh, Gary A Pope, Steven Bryant, Kamy Sepehrnoori,
and Larry W Lake. Reservoir simulation of co2 storage in deep saline aquifers. Spe
Journal, 10(03):336–348, 2005.
[74] Jorge L Landa and Baris Guyaguler. A methodology for history matching and the
assessment of uncertainties associated with flow prediction. In SPE Annual Technical
Conference and Exhibition?, pages SPE–84465. SPE, 2003.
111

[75] Abdeldjalil Latrach, Mohamed Lamine Malki, Misael Morales, Mohamed Mehana, and
Minou Rabiei. A critical review of physics-informed machine learning applications in
subsurface energy systems. arXiv preprint arXiv:2308.04457, 2023.
[76] B Li and F Friedmann. Novel multiple resolutions design of experiment/response
surface methodology for uncertainty analysis of reservoir simulation forecasts. In SPE
Reservoir Simulation Conference?, pages SPE–92853. SPE, 2005.
[77] Mingliang Liu, Divakar Vashisth, Dario Grana, and Tapan Mukerji. Joint inversion
of geophysical data for geologic carbon sequestration monitoring: A differentiable
physics-informed neural network model. Journal of Geophysical Research: Solid Earth,
128(3):e2022JB025372, 2023.
[78] Pascal Longuemare, M Mainguy, P Lemonnier, A Onaisi, Ch G´erard, and N Koutsabeloulis. Geomechanics in reservoir simulation: overview of coupling methods and field
case study. Oil & gas science and technology, 57(5):471–483, 2002.
[79] Y Zee Ma. Uncertainty analysis in reservoir characterization and management: How
much should we know about what we don’t know? 2011.
[80] Emmanuel Manceau, M Mezghani, I Zabalza-Mezghani, and F Roggero. Combination
of experimental design and joint modeling methods for quantifying the risk associated with deterministic and stochastic uncertainties-an integrated test study. In SPE
Annual Technical Conference and Exhibition?, pages SPE–71620. SPE, 2001.
[81] Rafael March, Florian Doster, and Sebastian Geiger. Assessment of co2 storage potential in naturally fractured reservoirs with dual-porosity models. Water Resources
Research, 54(3):1650–1668, 2018.
[82] Marcia McMillan, Robert Will, William Ampomah, Robert Balch, and Paige Czoski.
Coupled geomechanical modeling to assess cap rock integrity and mechanical fault stability: Application to farnsworth field unit project. In SPE Western Regional Meeting.
OnePetro, 2019.
[83] Saro Meguerdijian and Birendra Jha. 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, 2021.
[84] Saro Meguerdijian, Rajesh J Pawar, Bailian Chen, Carl W Gable, Terry A Miller,
and Birendra Jha. Physics-informed machine learning for fault-leakage reduced-order
modeling. International Journal of Greenhouse Gas Control, 125:103873, 2023.
[85] Saro Meguerdijian, Rajesh J Pawar, Dylan R Harp, and Birendra Jha. Thermal and
solubility effects on fault leakage during geologic carbon storage. International Journal
of Greenhouse Gas Control, 116:103633, 2022.
[86] Richard S Middleton, Gordon N Keating, Hari S Viswanathan, Philip H Stauffer, and
Rajesh J Pawar. Effects of geologic reservoir uncertainty on co2 transport and storage
infrastructure. International Journal of Greenhouse Gas Control, 8:132–142, 2012.
112

[87] Yadong Mu, Wei Liu, Xiaobai Liu, and Wei Fan. Stochastic gradient made stable:
A manifold propagation approach for large-scale optimization. IEEE Transactions on
Knowledge and Data Engineering, 29(2):458–471, 2016.
[88] Vinod Nair and Geoffrey E Hinton. Rectified linear units improve restricted boltzmann
machines. In Proceedings of the 27th international conference on machine learning
(ICML-10), pages 807–814, 2010.
[89] Long Nghiem, Peter Sammon, Jim Grabenstetter, and Hiroshi Ohkuma. Modeling co2
storage in aquifers with a fully-coupled geochemical eos compositional simulator. In
SPE Improved Oil Recovery Conference?, pages SPE–89474. SPE, 2004.
[90] Myeong Noh, Larry Wayne Lake, Steven Lawrence Bryant, and A Araque-Martinez.
Implications of coupling fractional flow and geochemistry for co2 injection in aquifers.
SPE Reservoir Evaluation & Engineering, 10(04):406–414, 2007.
[91] Yagna Deepika Oruganti and Srikanta Mishra. An improved simplified analytical
model for co2 plume movement and pressure buildup in deep saline formations. International Journal of Greenhouse Gas Control, 14:49–59, 2013.
[92] Indranil Pan, Masoud Babaei, Anna Korre, and Sevket Durucan. Artificial neural
network based surrogate modelling for multi-objective optimisation of geological co2
storage operations. Energy Procedia, 63:3483–3491, 2014.
[93] Peng-Zhi Pan, Jonny Rutqvist, Xia-Ting Feng, and Fei Yan. Modeling of caprock
discontinuous fracturing during co2 injection into a deep brine aquifer. International
Journal of Greenhouse Gas Control, 19:559–575, 2013.
[94] Ding-Yu Peng and Donald B Robinson. A new two-constant equation of state. Industrial & Engineering Chemistry Fundamentals, 15(1):59–64, 1976.
[95] W-J Plug and J Bruining. Capillary pressure for the sand–co2–water system under
various pressure conditions. application to co2 sequestration. Advances in Water Resources, 30(11):2339–2353, 2007.
[96] Bruce E Poling, John M Prausnitz, and John P O’connell. Properties of gases and
liquids. McGraw-Hill Education, 2001.
[97] Karsten Pruess and Julio Garcia. Multiphase flow dynamics during co 2 disposal into
saline aquifers. Environmental Geology, 42:282–295, 2002.
[98] Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville. On the spectral bias of neural networks.
In International Conference on Machine Learning, pages 5301–5310. PMLR, 2019.
[99] Nicolas Remy, Alexandre Boucher, and Jianbing Wu. Applied geostatistics with
SGeMS: A user’s guide. Cambridge University Press, 2009.
113

[100] Antonio P Rinaldi, Pierre Jeanne, Jonny Rutqvist, Fr´ed´eric Cappa, and Yves
Guglielmi. Effects of fault-zone architecture on earthquake magnitude and gas leakage
related to co2 injection in a multi-layered sedimentary system. Greenhouse Gases:
Science and Technology, 4(1):99–120, 2014.
[101] Antonio P Rinaldi, Jonny Rutqvist, and Fr´ed´eric Cappa. Geomechanical effects on co2
leakage through fault zones during large-scale underground injection. International
Journal of Greenhouse Gas Control, 20:117–131, 2014.
[102] RJ Rosenbauer and B Thomas. Carbon dioxide (co2) sequestration in deep saline
aquifers and formations. In Developments and innovation in carbon dioxide (CO2)
capture and storage technology, pages 57–103. Elsevier, 2010.
[103] Allen M Rowe Jr and James CS Chou. Pressure-volume-temperature-concentration
relation of aqueous sodium chloride solutions. Journal of Chemical and Engineering
Data, 15(1):61–66, 1970.
[104] J Rutqvist, JT Birkholzer, and Chin-Fu Tsang. Coupled reservoir–geomechanical analysis of the potential for tensile and shear failure associated with co2 injection in multilayered reservoir–caprock systems. International Journal of Rock Mechanics and
Mining Sciences, 45(2):132–143, 2008.
[105] Jonny Rutqvist, Donald W Vasco, and Larry Myer. Coupled reservoir-geomechanical
analysis of co2 injection and ground deformations at in salah, algeria. International
Journal of Greenhouse Gas Control, 4(2):225–230, 2010.
[106] Prasad Saripalli and Peter McGrail. Semi-analytical approaches to modeling deep well
injection of co2 for geological sequestration. Energy Conversion and Management,
43(2):185–198, 2002.
[107] Mohammad Sayyafzadeh. Reducing the computation time of well placement optimisation problems using self-adaptive metamodelling. Journal of Petroleum Science and
Engineering, 151:143–158, 2017.
[108] Klaus Schittkowski. The nonlinear programming method of wilson, han, and powell
with an augmented lagrangian type line search function: Part 1: Convergence analysis.
Numerische Mathematik, 38:83–114, 1982.
[109] Klaus Schittkowski. On the convergence of a sequential quadratic programming
method with an augmented lagrangian line search function. Mathematische Operationsforschung und Statistik. Series Optimization, 14(2):197–216, 1983.
[110] Jared Schuetter, Priya Ravi Ganesh, Doug Mooney, et al. Building statistical proxy
models for co2 geologic sequestration. Energy Procedia, 63:3702–3714, 2014.
[111] Drew R Schwab, Tandis S Bidgoli, and Michael H Taylor. Characterizing the potential
for injection-induced fault reactivation through subsurface structural mapping and
stress field analysis, wellington field, sumner county, kansas. Journal of Geophysical
Research: Solid Earth, 122(12):10–132, 2017.
114

[112] Hossein Shamshiri and Behnam Jafarpour. Controlled co2 injection into heterogeneous
geologic formations for improved solubility and residual trapping. Water Resources
Research, 48(2), 2012.
[113] Parisa Shokouhi, Vikas Kumar, Sumedha Prathipati, Seyyed A Hosseini, Clyde Lee
Giles, and Daniel Kifer. Physics-informed deep learning for prediction of co2 storage
site response. Journal of Contaminant Hydrology, 241:103835, 2021.
[114] Per Arne Slotte and E Smørgrav. Response surface methodology approach for history
matching and uncertainty assessment of reservoir simulation models. In SPE Europec
featured at EAGE Conference and Exhibition?, pages SPE–113390. SPE, 2008.
[115] Yen A Sokama-Neuyam, Wilberforce N Aggrey, Patrick Boakye, Kwame Sarkodie,
Sampson Oduro-Kwarteng, and Jann R Ursin. The effect of temperature on co2 injectivity in sandstone reservoirs. Scientific African, 15:e01066, 2022.
[116] Jerzy Stopa, Damian Janiga, Pawe l Wojnarowski, and Robert Czarnota. Optimization
of well placement and control to maximize co2 trapping during geologic sequestration.
AGH Drilling, Oil, Gas, 33(1):93–104, 2016.
[117] Anna L Stork, James P Verdon, and J-Michael Kendall. The microseismic response at
the in salah carbon capture and storage (ccs) site. International Journal of Greenhouse
Gas Control, 32:159–171, 2015.
[118] Fahad I Syed, Abdulla AlShamsi, Amirmasoud K Dahaghi, and Shahin Neghabhan.
Application of ml & ai to model petrophysical and geomechanical properties of shale
reservoirs–a systematic literature review. Petroleum, 8(2):158–166, 2022.
[119] Meng Tang, Xin Ju, and Louis J Durlofsky. Deep-learning-based coupled flowgeomechanics surrogate model for co2 sequestration. International Journal of Greenhouse Gas Control, 118:103692, 2022.
[120] Pietro Teatini, Giuseppe Gambolati, Massimiliano Ferronato, A Tony Settari, and
Dale Walters. Land uplift due to subsurface fluid injection. Journal of Geodynamics,
51(1):1–16, 2011.
[121] David Tran, Long Nghiem, and Antonin Settari. New iterative coupling between a
reservoir simulator and a geomechanics module. Spe Journal, 9(03):362–369, 2004.
[122] David Tran, Long Nghiem, Vijay Shrivastava, and Bruce Kohse. Study of geomechanical effects in deep aquifer co2 storage. In ARMA US Rock Mechanics/Geomechanics
Symposium, pages ARMA–10. ARMA, 2010.
[123] Minh Tran and Birendra Jha. Effect of poroelastic coupling and fracture dynamics on solute transport and geomechanical stability. Water Resources Research,
57(10):e2021WR029584, 2021.
[124] Jorn FM Van Doren, Renato Markovinovi´c, and Jan-Dirk Jansen. Reduced-order optimal control of water flooding using proper orthogonal decomposition. Computational
Geosciences, 10:137–158, 2006.
115

[125] James P Verdon, J-Michael Kendall, Anna L Stork, R Andy Chadwick, Don J White,
and Rob C Bissell. Comparison of geomechanical deformation induced by megatonnescale co2 storage at sleipner, weyburn, and in salah. Proceedings of the National
Academy of Sciences, 110(30):E2762–E2771, 2013.
[126] JP Verdon, J-M Kendall, DJ White, and DA Angus. Linking microseismic event observations with geomechanical models to minimise the risks of storing co2 in geological
formations. Earth and Planetary Science Letters, 305(1-2):143–152, 2011.
[127] PA Vermeer. Non-associated plasticity for soils, concrete and rock. In Physics of dry
granular media, pages 163–196. Springer, 1998.
[128] St´ephanie Vialle, Jennifer L Druhan, and Katharine Maher. Multi-phase flow simulation of co2 leakage through a fractured caprock in response to mitigation strategies.
International Journal of Greenhouse Gas Control, 44:11–25, 2016.
[129] Sandrine Vidal-Gilbert, Jean-Francois Nauroy, and Etienne Brosse. 3d geomechanical modelling for co2 geologic storage in the dogger carbonates of the paris basin.
International Journal of Greenhouse Gas Control, 3(3):288–299, 2009.
[130] Victor Vilarrasa, Sebastia Olivella, and Jesus Carrera. Geomechanical stability of the
caprock during co2 sequestration in deep saline aquifers. Energy Procedia, 4:5306–5313,
2011.
[131] Rafael Villamor Lora, Ehsan Ghazanfari, and Enrique Asanza Izquierdo. Geomechanical characterization of marcellus shale. Rock mechanics and rock engineering,
49:3403–3424, 2016.
[132] Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy,
David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan
Bright, et al. Scipy 1.0: fundamental algorithms for scientific computing in python.
Nature methods, 17(3):261–272, 2020.
[133] G Waters, J Heinze, Randy Jackson, A Ketter, J Daniels, and Doug Bentley. Use of
horizontal well image tools to optimize barnett shale reservoir exploitation. In SPE
Annual Technical Conference and Exhibition?, pages SPE–103202. SPE, 2006.
[134] Gege Wen, Zongyi Li, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M
Benson. U-fno—an enhanced fourier neural operator-based deep-learning model for
multiphase flow. Advances in Water Resources, 163:104180, 2022.
[135] Christopher D White and Steve A Royer. Experimental design as a framework for
reservoir studies. In SPE Reservoir Simulation Conference?, pages SPE–79676. SPE,
2003.
[136] Elizabeth J Wilson, Timothy L Johnson, and David W Keith. Regulating the ultimate
sink: managing the risks of geologic co2 storage. Environmental Science & Technology,
37(16):3476–3483, 2003.
116

[137] Bicheng Yan, Dylan Robert Harp, Bailian Chen, and Rajesh Pawar. A physicsconstrained deep learning model for simulating multiphase flow in 3d heterogeneous
porous media. Fuel, 313:122693, 2022.
[138] Burak Yeten, Alexandre Castellini, Baris Guyaguler, and WH Chen. A comparison
study on experimental design and response surface methodologies. In SPE Reservoir
Simulation Conference?, pages SPE–93347. SPE, 2005.
[139] Junjie Yu, Atefeh Jahandideh, and Behnam Jafarpour. Efficient robust production
optimization with reduced sampling. SPE Journal, 27(04):1973–1988, 2022.
[140] Zheming Zhang and Ramesh Agarwal. Numerical simulation and optimization of co2
sequestration in saline aquifers. Computers & Fluids, 80:79–87, 2013.
[141] Xiaoxi Zhao and Birendra Jha. Role of well operations and multiphase geomechanics
in controlling fault stability during co2 storage and enhanced oil recovery. Journal of
Geophysical Research: Solid Earth, 124(7):6359–6375, 2019.
[142] Xiaoxi Zhao and Birendra Jha. A new coupled multiphase flow–finite strain
deformation–fault slip framework for induced seismicity. Journal of Computational
Physics, 433:110178, 2021.
[143] F Zheng, B Jha, and B Jafarpour. Controlled co2 injection into storage reservoirs to
minimize geomechanical risks under geologic uncertainty. In SPE Annual Technical
Conference and Exhibition?, page D021S020R005. SPE, 2023.
[144] Fangning Zheng, Atefeh Jahandideh, Birendra Jha, and Behnam Jafarpour. Geologic
co2 storage optimization under geomechanical risk using coupled-physics models. International Journal of Greenhouse Gas Control, 110:103385, 2021.
[145] Zhi Zhong, Alexander Y Sun, and Hoonyoung Jeong. Predicting co2 plume migration
in heterogeneous formations using conditional deep convolutional generative adversarial network. Water Resources Research, 55(7):5830–5851, 2019.
[146] Mark D Zoback. Reservoir geomechanics. Cambridge university press, 2010.
[147] Mark D Zoback and Steven M Gorelick. Earthquake triggering and large-scale geologic storage of carbon dioxide. Proceedings of the National Academy of Sciences,
109(26):10164–10168, 2012.
117

Appendices
A Evolutionary Multi-objective Optimization (EMO)
Method
Multi-objective evolutionary (MOE) algorithms have largely sought to characterize nondominated solutions in an evolutionary population. A general flowchart of an evolutionary
algorithm is shown on Figure A.1 [1]. The algorithm starts with a feasible initial population and proceed to the evaluation step to assign a fitness value to each individual in the
initial population. If the convergence criteria are not met, the algorithm proceeds to the
reproduction step which includes selection, crossover and mutation to create the offspring
generation. The iteration continues until the convergence criteria are satisfied. In this section, we describe one of the state-of-the-art EMO solutions approaches using the genetic
algorithm (GA), known as the Non-dominated Sorting GA or NSGA-II procedure [31] [32].
A.1 Elitist Non-dominated Sorting GA (NSGA-II)
The NSGA-II has three main features that enable fast convergence and characterization
of the Pareto-front: elitist principle application, diversity preservation, and non-dominated
118

Figure A.1: Flowchart of multi-objective genetic algorithm.
solutions utilization [32]. We begin with a summary of the important definitions and terminologies used in MOE algorithms. A point x dominates a point y for objective functions f
if
fi(x) ≤ fi(y) for all i, (A.1a)
fj (x) < fj (y) for some j. (A.1b)
The individuals on the feasible solution space are iteratively ranked based on domination to
identify the Pareto-front. Rank 1 individuals are non-dominated by any other individuals
119

and Rank 2 individuals are only dominated by Rank 1 individuals and so on. The crowding
distance is a measure of how far an individual is from its surrounding neighbors for dimension
m and can be calculated as
distancei =
X
m
x(m, i + 1) − x(m, i − 1), (A.2)
where distancei
is the crowding distance of an individual xi
from its sorted neighbors xi+1
and xi−1 in each dimension. The spread is a measure of the movement of the Pareto set and
is calculated as
spread = (µ + σ)/(µ + Qd) (A.3)
where µ is the sum over the k objective function indices of the norm of the difference between
the current minimum-value Pareto point for that index and the minimum point for that index
in the previous iteration, σ is the standard deviation of the crowding distance measure of
points that are on the Pareto-front with finite distance, Q is the number of these points, and
d is the average distance measure among these points.
Figure A.2 (a) illustrates the NSGA-II procedure while Table A.1 Presents a pseudo code
for it [32]. The algorithm starts with a parent population (Pt), which is used to create an
offspring population (Qt) by selection, crossover and mutation operators [47]. The parent
population (Pt) are then combined with the offspring population (Qt) to form an entire
population (Rt) that has a size of 2N. A non-dominated sorting is performed to rank the
individuals in Rt and fill in a new population of size 2N, starting with individuals from the
best non-dominated front (F1) and continuing with the second non-dominated front (F2)
until the worst non-dominated front (Fk). Next, by creating a new empty population (Pt+1)
of size N, individuals in Rt are accommodated into Pt+1 according to front ranking starting
from F1. All the fronts which cannot be counted into Pt+1 due to the limiting size N will
be discarded. If there exists a last accepted front where only part of the individuals can be
filled into Pt+1, a crowding distance sorting will be performed to sort the individuals in this
120

Table A.1: pseudocode of NSGA-II
Rt = Pt ∪ Qt
F = Fast-non-dominated-sort(Rt),
Pt+1 = ϕ and i = 1
until |Pt+1| + |Fi ≤ N
Crowding-distance-assignment(Fi)
Pt+1 = Pt+1 ∪ F1
i = i + 1
Sort(Fi
, ≺n)
Pt+1 = Pt+1 ∪ Fi
[1 : (N − |Pt+1|)]
Qt+1 = Make-new-pop(Pt+1)
t = t + 1
Combine parent and offspring population
F = (F1, F2, . . .) all non-dominated fronts of Rt
sorted by non-dominated sorting until the parent
population is filled
calculate crowding-distance in Fi
include i-th non-dominated front in the parent pop
sort in descending order using crowded-comparison
operator ≺n
choose the first (N − |Pt+1) elements of Fi use selection,
crossover and mutation to create a new population
Qt+1
increment the generation counter
front according to the spreading distance on the front. The most widely spread individuals
are counted into Pt+1 to preserve population diversity until total number of N individuals in
Pt+1 is reached. The new parent population (Pt+1) is formed and a new offspring population
(Qt+1) is created from Pt+1 by selection, crossover, and mutation.
A.2 Controlled NSGA-II
In NSGA-II, all the elite solutions are preserved and accepted as solutions in the population.
This could create a situation where too many elite solutions may prevent new solutions from
being accepted as the parent population. When this scenario happens, the global search
algorithm may converge to a local optimal solution set. To deal with this issue, Deb and Goel
(2001) [31] proposed the controlled NSGA-II to restrict the number of elite solutions in the
parent population. They proposed to use a geometric distribution to maintain the number
of individuals in each non-dominated front. In the controlled NSGA-II algorithm, after the
non-dominated sorting step, the number of individuals (ni) allowed to be accommodated
in the new parent population (Pt+1) in the i-th non-dominated front (i = 1, 2, . . . , k) is
calculated as
ni = N
1 − r
1 − r
k
r
i−r
(A.4)
121

where N is the size of Pt+1, r(< 1) is the reduction rate which is a user defined parameter. In
this way, the new parent population in controlled NSGA-II has an improved diversity than
that in NSGA-II (Figure A.2 (b)).
Rejected
Non-dominated
Sorting Crowding Distance
Sorting
Pt
Qt
Rt
ℱ1
Pt+1
…
ℱ!
ℱ"
ℱ#
Pt
Rt Pt+1
Controlled NSGA-II
Pt+1
NSGA-II
…
ℱ1
ℱ!
ℱ"
ℱ#
(a) Procedure of NSGA-II.
(b) Comparison between Controlled NSGA-II and NSGA-II.
Figure A.2: NSGA-II algorithm.
122

B Preliminary Tests on Numerical Simulations and Objective
Functions
B.1 Model 1: Gaussian-type
A preliminary test on the coupled fluid flow and geomechanics simulation model is performed
using CMG-GEM to test the solution accuracy for Model 1 (2.4.1). The sensitivity analysis
of this model on the aquifer settings as well as the model grid resolution for both the
geomechanical domain and the reservoir domain are included in Appendix C. Figure B.1
shows the simulation results of the preliminary test for a three-well, 15-year CO2 injection
example of well locations on the injection layer permeability map, gas plume, pressure,
effective mean stresses, Von Mises stresses (VMS), effective plastic strain, and displacement
in x, y, and z directions. The total CO2 injection is 6.23 mega-ton for this specific well
configuration and the highest reservoir pressure reaches to around 23800 kPa. The maximum
vertical displacement at the end of injection is around 0.046 meter. Figure B.2 shows the
cumulative gas mass (mega-ton) and injection block pressure for each well. During injection,
the CO2 plume of a well evolves and potentially interferes with other plumes. The effective
plastic strain map in Figure B.1 (f) shows that the maximum plastic deformation occurs
between Well 1 and Well 2. We illustrate this behavior in Figure B.3. The threshold,
effective plastic strain > 0, is used to show only failed grid blocks. The invisible grid blocks
in Figure B.3 are grid blocks with zero plastic strain. In Figure B.3, we select two groups
of grid blocks to plot the stress path. The first group has three grid blocks (Well 1 top grid
block, Well 2 top grid block and top ‘a’ grid block). The stress paths for these grid blocks
are shown in Figure B.4 (a). In Figure B.4, we plot the Von Mises stress (VMS) vs. effective
mean stress (I1) which can be calculated, respectively, as:
1
3
V MS2 = J2 =
1
6
h
(σ
′
11 − σ
′
22)
2 + (σ
′
22 − σ
′
33)
2 + (σ
′
33 − σ
′
11)
2
i
+ σ
2
12 + σ
2
23 + σ
2
31 (B.1)
123

and,
I1 =
1
3
(σ
′
11 + σ
′
22 + σ
′
33) (B.2)
The Von Mises stress is a measure of the deviatoric stress (responsible for plastic failure)
and the effective mean stress is a measure of the volumetric stress (responsible for volumetric
expansion/contraction). From Figure B.4 (a), the three stress paths start from the same
point (labeled year 0) because the initial stresses are the same for grid blocks at the same
layer in the reservoir. When injection starts, the stress paths evolve until they reach the
failure envelope. The rate of change of VMS for ‘a’ is faster than that of Well 1 and Well
2, meaning that the deviatoric stress change is the fastest at location ‘a’. The increase in
deviatoric stress is mainly due to the contribution from both Wells 1 and Well 2 and leads
to an earlier contact with the failure envelope (at around Year 11 for ‘a’). As injection
continues, plastic strain continues to accumulate, and the stress path remains on the failure
envelope until injection stops.
For Well 1 and Well 2, the stress path touches the failure envelope at around Year 12
and remains on the failure envelope until injection stops. Figure B.4 (b) shows the stress
paths for Well 1 at the reservoir top layer and at the injection layer. For the injection grid
block, the VMS evolution is faster than that of the top grid block. However, the top grid
block touches the failure line faster than the injection grid block because the initial stresses
at the top block are closer to failure.
B.2 Model 3: CO2 Storage with Caprock and Aquifer layers
In this section, we perform a preliminary test on the coupled flow-geomechanics-fracturing
model presented in section 2.4.3 using CMG-GEM for the solution accuracy check. The
purpose of the preliminary test is to qualitatively validate the results of the simulation
model for the CO2 caprock leakage system. Figure B.5 shows the simulation results of the
preliminary test for a one-well, 15-year CO2 injection scenario using the model described in
124

Figure B.1: Simulation results for a three-well CO2 injection example. Well locations on the
injection layer’s permeability map, gas plume shape after injecting a certain mass, pressure,
effective mean stress, Von Mises Stress, effective plastic strain, and displacements in the x-,
y-, and z-directions.
section 2.4.3. The injector is placed in the middle grid block (block index {19, 15, 27}) of
the injection reservoir at the 27th layer (z=27 in the reservoir domain) at a depth of 1042.5
m. The grid shown in Figure B.5 is the injection reservoir subregion. The well location on
a log10(permeability) map at the injection layer, gas plume, pressure field, z-displacement,
effective mean stress, and the Von Mises stress are plotted. The total CO2 mass injected
is about 6 mega-ton. The maximum pressure is 17817 kPa, and the maximum vertical
displacement on the ground surface is 0.056 m.
Figure B.6 (a), (b), (c), and (d) shows the CO2 gas plume in the upper aquifer (upper
red blocks), middle caprock fractures (middle blue blocks), and the lower injection reservoir
(lower red blocks) at year 1, 5, 10, and 15, respectively. Figure B.6 (e), (f), (h), and (h)
125

Figure B.2: Plots of (a) cumulative gas mass vs. injection time and (b) injection block
pressure vs. injection time for each of the three CO2 injection wells. The three wells are
injecting under BHP control of 24000 kPa for 15 years.
Figure B.3: Plots of (a) 3D plot of effective plastic strain with a threshold > 0 (only showing
grid blocks that failed). Wells are labeled by the black stick and numbers from ‘1’ to ‘3’.
The blue stick labeled ‘a’ is used for comparison later. (b) sliced 3D plot on x-direction. It
shows that the failure region is shaped like a bowl with larger and more widespread failure
observed at shallower depths.
represents the corresponding years of fracture permeability in the z-direction (kfz), respectively. We observe caprock fracturing at year 1 and gas leaks from the injection reservoir
to the caprock’s fractures. As injection continues, higher pressure build-up in the reservoir
occurs to allow more fractures to open in the caprock layers. When the fractures reach the
126

Year 0
Year 15 Failure Line Failure Line
Year 15
Year 0
(a) (b)
Figure B.4: Von Mises stress vs. mean effective normal stress plots showing the stress
trajectories of selected grid blocks. Marker locations correspond to time steps in years,
starting from Year 0 (labelled on the plot) to Year 15. (a) black circle is the stress path for
the reservoir top grid block of well 1. Red star identifies the stress path of the reservoir top
grid block of ‘a’ in Figure 5. Blue square identifies the stress path of the reservoir top grid
block of well 2. (b) black square identifies the stress path of the reservoir top grid block of
well 1. Black circle identifies the stress path of the reservoir injection grid block of well 1.
top layer of the caprock, a leakage pathway is formed to allow CO2 leaking from the injection
reservoir to the upper aquifers (Year 10).
To better understand the changes of pressure and stress states of fracturing, we define an
observation point called ‘CapPoint1’ located at the bottom of the caprock layer above the
injector at gridblock having index of {19, 15, 10} (the injector block index is {19, 15, 27}).
Figure B.7 shows a line plot of the fracture permeability, minimum effective stress, and pressure vs. time for ‘CapPoint1’. As pressure increases due to injection, fracturing is initiated
in the grid block at the 2nd month after injection starts. Fracture pressure immediately
increases from about 9700 kPa to 15500 kPa. Mean total stress suddenly increases from
about 18500 kPa, to over 20500 kPa and the effective minimum stress decreases from about
8880 kPa to below 8,300 kPa, which is the fracture opening stress. In the year 2004, month
10 (four years and ten months after injection started), another fracture was initiated nearby
the ‘CapPoint1’. This resulted in a sudden pressure increase in the newly fractured grid
block and allowed some pressure release in the ‘CapPoint1’. The effective minimum stress
127

Figure B.5: Preliminary test results on the injection layers for a one-well, 15-year injection
coupled flow-geomechanics-fracturing simulation.
rises to above 8300 kPa, and fractures start to close. Since it cannot return completely to
the original state, the fracture permeability reaches the fracture closure permeability of 33
mD. As injection continues after the year 2004 month 10, pressure continues to increase until
the secondary fracturing event happens at the ‘CapPoint1’. Fracture permeability suddenly
rises to the fracture opening permeability of 100 mD again and remains open afterward. The
preliminary test of the simulation model successfully demonstrates the fracturing mechanism
and the application of the modified Barton-Bandis model described in section 2.1.4.
B.3 Model 3: Objective Function Behavior
We examine the objective function behavior before performing optimization to show the
behavior of the new equation (Eqn. 2.24) we developed. We place a single injector in the
middle of the reservoir (Figure 2.7), injecting CO2 for 15 years. Figure B.8 (a) and (b)
128

Figure B.6: CO2 Gas plumes (a,b,c,d) and kfz (e,f,g,h) at year 1, 5, 10, and 15, respectively.
The 3D plots are in x − z view. In CO2 Gas plume plots (a,b,c,d), the upper red blocks
represent the gas plume in the upper aquifer, the middle blue blocks represent the gas plume
in the caprock fractures, and the lower red blocks represent the gas plume in the injection
reservoir. For plots of kfz (e,f,g,h), the z-axis range is from 875 to 1000 meters in depth,
which is the caprock layers. The colored grid blocks represent the locations of caprock
fractures where the kfz has reached 100 mD.
shows a plot of the objective function f(u, m) vs. injection rate u of the injector located
at index {19, 15, 27}. Figure B.8 (a) shows the data points inside the blue dashed-line
box of Figure B.8 (b) on a different y-scale. From the figure, we observe that for a single
well scenario, as the injection rate increases, the objective function value increases. Caprock
fracturing occurs when the injection rate increases to 30000 m3/day. For the cases of injection
rates lower than 30000 m3/day, the objective function f(u, m) has a linearly increasing trend.
While for the cases of injection rates that are equal to or higher than 30000 m3/day, the
objective function f shows an exponentially increasing trend, which is in accordance with
the effect of the barrier term in the objective function. From Figure B.8 (a) and (b), we
conclude that the greater the total injection volume, the larger the likelihood and extent of
the caprock fracturing for a single injector scenario.
129

Figure B.7: Time evolution of kfz, minimum Effective Stress, and pressure for ‘CapPoint1’.
To further investigate the non-linearity for simulations with caprock fracturing, we plot
the corresponding pressure and minimum effective stress for each data point of Figure B.8 (ab) in Figure B.8 (c-f). In Figure B.8 (c-f), each color of data points represents the pressure/minimum effective stress for a given simulation at injection location (grid index={19, 15, 27},
blue dots) or caprock location (grid index={19, 15, 10}, orange dots) at year 15. From the
plots, we observe that the pressure or the minimum effective stress of the first 11 data points
with small injection rates show a linear trend corresponding to the linear trend shown in
Figure B.8 (a-b) for simulations without caprock fracturing. Pressure shows a monotonically
increasing trend, while the minimum effective stress shows a monotonically decreasing trend,
and the monotonicity is because the rock continues to be in an elastic state. We observe a
non-linear trend induced by caprock fracturing shown by the last five data points for both
pressure and minimum effective stress due to the inelastic behavior of the rock undergoing
fracturing.
130

Fracture Opens:
exponential trend
No fracturing: linear
trend
No fracturing:
linear trend
linear
linear
non-linear
linear
non-linear
linear
non-linear
u=1e4 m3/day
u=7e5 m3/day
u=7e5 m3/day
u=1e4 m3/day
linear
non-linear
19
19
(c) (d) (e) (f)
Figure B.8: Objective function f vs. injection rate u. (a) Before fracture, a near-linear trend
is observed. (b) After fracture opens at caprock, the trend of objective function becomes
exponential. (a) is a plot of the data points inside the blue dashed-line box of (b) on a
different y-scale. (c) Pressure at year 15 vs. injection rate u at the 27th layer (injection
reservoir) for the blue dots and at the 10th layer (the bottom layer of the caprocks) for the
orange dots, respectively. (d) Effective minimum stress at year 15 vs. injection rate u at
the 27th layer (blue dots) and the 10th layer (orange dots), respectively. (e) Pressure vs.
effective minimum stress at year 15 at the 27th layer (blue dots). (f) Pressure vs. effective
minimum stress at year 15 at the 10th layer (orange dots). The y-axes are correlated with
the same colors of points.
Figure B.9 plots the pressure and the minimum effective stress vs. time for each simulation, respectively, with different injection rates indicated in different colors at grid {19, 15, 27},
which is the injector’s perforation location. From the figure, we observe that for the first 11
injection rates (u=10000 to 200000 m3/day), the plots of pressure vs. time and minimum
effective stress vs. time present a linear and monotonically increasing trend. This is because
the reservoir rock remained within its elastic deformation range and did not experience any
fracturing. This explains the linear trends in Figure B.8. For the last five injection rates
(u=300000 to 700000 m3/day), the plots of pressure and minimum effective stress vs. time
are showing a non-linear and non-monotonic trend. This is because caprock fracturing has
occurred during the injection. The fracturing of the caprock will open a highly permeable
131

pathway within caprock layers, the pressure of the reservoir will be released, and this will
cause the minimum effective stress to increase (become more compressive) as per Eqn. 2.3 to
result in the non-linear and non-monotonic behavior for both pressure and effective stresses.
For the case of u=700000 m3/day, the abrupt release of high pressures and stresses is an
indication of caprock fracturing, which has opened a pathway for leakage to the shallow
aquifer due to a large amount of CO2 injection.
Linear and
monotonically
increasing
Non-linear and
non-monotonic
Non-linear and
non-monotonic
Linear and
monotonically
decreasing
Figure B.9: Plots of pressure (a) and minimum effective stress (b) at grid {19, 15, 27} vs.
time. Different colors of lines are representing different simulation results with the corresponding injection rate u specified with the legend on the right.
A complicating factor in coupled-physics models is the complexity and non-linearity of
the simulation model. For a given grid block, if the minimum stress reaches the fracture
opening stress, a fracture in the given grid block opens, and the pressure of the grid block
will increase suddenly. The pressure of the surrounding rock will decrease because of the
pressure release through the fracture. As injection continues, fracture initiation and growth
occur in other places. In this case, the opening of a fracture at a point will cause the
pressure of the surrounding region to decrease and might also cause the minimum stress
of some gridblocks to fall below the fracture closure stress. The fracture closure will then
increase the pressure of the surrounding gridblocks to enable additional fracture openings in
the surrounding gridblocks.
132

C Sensitivity Analysis of Numerical Simulation Model
In this section, we present the sensitivity analysis with respect to model grid size for the
aquifer in the two numerical simulation examples in section 2.4.1 and section 2.4.2. The main
purpose is to select a reasonable aquifer grid size that balances accuracy and computationally
efficiency. Table C.1 summarizes different cases that are considered for sensitivity analysis of
aquifer grid resolution for numerical models 1 and 2. The x and y dimensions of the aquifer
for numerical model 1 is 4000-meter by 1500-meter on the two sides of the reservoir along
the x-direction. We perform the sensitivity analysis with different numbers of grid block and
sizes. We use CMG GEM’s grid refinement capability to refine the aquifer grid based on the
largest grid size of Case 1 in Table C.1. The simulation model settings are the same as in
Table 2.1 and Table 2.2, except for the aquifer grid resolutions. For comparison, we perform
single-well CO2 injection simulation under a constant injection rate of 105 m3/day (surface
condition) for 15 years. The injection well is located at the center of the reservoir for the
cases listed in Table C.1. We calculate and plot the relative errors in pressure and effective
mean stress (I1) in Figure C.1.
The relative error for the entire simulation grid is defined as
Error% = ∥Pi − Pf inest∥2
∥Pf inest∥2
× 100, (C.1)
where Pi
is the reservoir pressure or mean stress I1 of Case i and Pf inest is the reservoir
pressure or mean stress I1 of the case with the finest grid resolution. The relative error for
the injection grid block is defined as
Injection Region Error% = ∥Pinj,i − Pinj,f inest∥2
∥Pinj,f inest∥2
× 100, (C.2)
where Pinj,i is the Pressure or stress I1 of Case i at the injection block and Pinj,f inest is the
Pressure or mean stress I1 of the case that has the finest grid resolution at the injection
block.
133

Figure C.1 shows the bar chart of pressure and effective mean stress errors for different
cases of model 1. The relative errors in pressure for aquifer grid block numbers greater
than 8 (4+4 on two side of reservoir in x-direction) are less than 0.5% and the relative
errors in effective mean stress for aquifer grid block numbers greater than 8 are less than
0.2%. Figure C.2 shows the cross-sectional plots of pressure and effective mean stress at the
injection layer across the injection well located at the center of the reservoir for model 1.
We observe a converging trend in both pressure and effective mean stress as the aquifer grid
resolution increases towards the finest grid. This indicates that the solution accuracy does
not improve significantly for aquifer grid block numbers greater than 20. The aquifer grid
block number 20 is then used as the reference case to calculate the relative errors. The choice
of aquifer grid resolution is a trade-off between solution accuracy and simulation runtime.
To save the computational cost of optimization, we use an aquifer grid block number of 8 to
have reasonable runtime and acceptable solution error.
The same sensitivity analysis procedure is performed for numerical model 2 in section 2.4.2. Different cases of grid resolution are summarized in Table C.1. Figure C.3 show
the pressure and I1 errors for difference cases, respectively. From these two plots, we observe
that relative errors are less than 0.2% for the case with grid block numbers greater than 12
(6+6). Figure C.4 shows the cross-sectional plots of pressure and effective mean stress at
the injection layer across the injection well, which is located at the center of the reservoir. A
similar converging trend in pressure and effective mean stress results is observed when the
aquifer grid resolution increases towards the finest grid. This suggests that no significant
improvement in accuracy will be made for grid block numbers greater than 30. We choose
an aquifer grid block number of 24 (12+12) for the simulation model in numerical model 2
as a reasonable trade-off between simulation runtime and accuracy.
134

Table C.1: Aquifer grid resolution for model 1 in section 2.4.1 and model 2 in section 2.4.2.
Model 1 Model 2
Aquifer Grid Resolutions in x-direction
number of grid blocks × grid size (meter)
on left- and right-sides of injection region
Case 1 1×4000+1×4000 2×1200+2×1200
Case 2 2×2000+2×2000 4×600+4×600
Case 3 4×1000+4×1000 6×400+6×400
Case 4 6×666.67+6×666.67 8×300+8×300
Case 5 8×500+8×500 10×240+10×240
Case 6 10×400+10×400 12×200+12×200
Case 7 15×266.67+15×266.67 15×160+15×160
Case 8 20×200+20×200 -
135

(a)
(b)
Figure C.1: (a) Plot of the pressure relative error and simulation run time for the sensitivity
analysis of aquifer grid resolution for numerical model 1 (section 2.4.1). The x-axis represents
different cases of aquifer grid resolution listed in Table C.1. The left y-axis represents the
percentage error relative to the finest case (20×200). The right y-axis represents simulation
CPU run time. (b) Plot of the mean effective stress (I1) relative error from the sensitivity
analysis of aquifer grid resolution for numerical model 1. The x-axis represents different
cases of aquifer grid resolution listed in Table C.1. The y-axis represents the percentage
error relative to the finest case (20×200).
136

A
B
Figure C.2: Numerical Model 1 (section 2.4.1). Cross-sectional plots of (a) pressure and (b)
mean effective stress I1 at the injection layer of the reservoir for different grid resolutions.
Position of the cross-section A-B in 3D grid is shown in (c).
137

(a)
(b)
Figure C.3: (a) Plot of pressure relative error and simulation run time for the sensitivity
analysis of aquifer grid resolution for numerical model 2 (section 2.4.2). The x-axis represents
different cases of aquifer grid resolution listed in Table C.1. The left y-axis represents the
percentage error relative to the finest case (15×160). The right y-axis represents simulation
CPU run time the cases. (b) Plot of I1 relative errors for the sensitivity analysis of aquifer
grid resolution for numerical model 2. The x-axis represents different cases of aquifer grid
resolution listed in Table C.1. The y-axis represents the percentage error relative to the
finest case (15×160).
138

A
B
Figure C.4: Numerical model 2 (section 2.4.2). Cross-sectional plots of (a) pressure and
(b) I1 at the injection layer of the reservoir for different grid resolutions. Position of the
cross-section A-B in the 3D grid is shown in (c) which is extended along the x-axis from x=0
to x=15 km. Grid sizes are not to scale.
139

D Sensitivity Analysis of Surrogate Models
D.1 Training Data Size of Multi-model Learning
In this section, we present a sensitivity analysis of the proposed surrogate model concerning
the determination of training data size. The training is conducted using the U-Net models outlined in Section 2.5 and the training approach described in Section 2.6.1 to predict
the final-time step values for pressure, CO2 plume, vertical displacement (Dz), and effective
plastic strain (ε
′
p
). We explore the impact of varying training data sizes, ranging from 100
to 1000 data points, while maintaining a consistent 4:1 ratio between training and testing
datasets. For testing purposes, we employ a fixed set of 200 additional, previously unseen
data points for each training scenario. Table D.1 provides a comprehensive summary of
the Mean Squared Error (MSE) for different prediction parameters, encompassing pressure,
CO2 plume saturation, vertical displacement, and effective plastic strain. We also investigate
the influence of including or excluding permeability information. Concurrently, Figure D.1
depicts a graphical representation of the MSE error in relation to the training data size for
all four parameters. Upon examining the results presented in Table D.1 and Figure D.1,
it becomes evident that there is a consistent trend of decreasing MSE error as the training
data size increases for all four parameters. Specifically, a convergence in the MSE error is
observed for pressure and vertical displacement when the training size reaches 500, suggesting that little additional improvement in training performance is achieved beyond this point.
Therefore, a training size of 500 data points is deemed sufficient to yield highly accurate
testing results for predicting pressure and vertical displacement. Similarly, the convergence
for CO2 saturation and effective plastic strain is observed at a training size of 700, indicating
that there is minimal room for improvement in training performance beyond this threshold.
Consequently, we establish a minimum data size requirement of 700 data points for training models aimed at predicting CO2 saturation and effective plastic strain. Additionally,
the availability of permeability information played a crucial role in enhancing prediction
140

Table D.1: A summary of the Mean Squared Error (MSE) for four U-Net models including mapping to pressure, CO2 plume/saturation, vertical displacement, and effective plastic
strain, considering different training data sizes and the presence (y) or absence (n) of permeability information.
accuracy. When permeability information was not included, the predictive performance for
pressure, vertical displacement, and CO2 saturation experienced a decline.
D.2 Training Data Size of Multi-output Learning
In this section, a sensitivity analysis of the multi-output learning concerning the determination of training data size is presented. The training is conducted using the U-Net models
outlined in Section 2.5 and the training approach described in Section 2.6.2 to predict the
final-time step values for pressure, CO2 plume, vertical displacement (Dz), and effective
plastic strain (ε
′
p
). The weights for the training input of well location and permeability are
10 and 1, respectively (Case 2 in Table 5.2). We explore the impact of varying training data
sizes, ranging from 100 to 1000 data points, while maintaining a consistent 4:1 ratio between
141

0.0E+00
2.0E-06
4.0E-06
6.0E-06
8.0E-06
1.0E-05
1.2E-05
1.4E-05
1.6E-05
0.0E+00
5.0E-04
1.0E-03
1.5E-03
2.0E-03
2.5E-03
3.0E-03
100 300 500 700 900
MSE Loss
Training Size
Multi-model Learning Model Sensitivity Analysis
Series1
Series2
Series3
Series4 Effective Plastic Strain (�!"#
$ )
Vertical Displacement (�%)
Pressure
CO2 Saturation
Figure D.1: Multi-model learning model sensitivity analysis on training sample size. The
y-left axis is the MSE loss plotted for CO2 plume, vertical displacement, and pressure while
the y-right axis is the MSE error plotted for plastic strain.
training and testing datasets. For testing purposes, we employ a fixed set of 200 additional,
previously unseen data points for each training scenario.
The results of the sensitivity analysis are summarized in Table D.2 and plotted in Figure D.2. The plot reveals that the MSE losses for vertical displacement, pressure, and effective plastic strain exhibit convergent behavior when the training dataset reaches a size of 500.
In contrast, the MSE losses for CO2 saturation begin to converge at a training dataset size of
900. This suggests that a minimum dataset size of 900 is required for achieving satisfactory
performance with the proposed U-Net model and training methodology. Figures D.3, D.4,
D.5, and D.6 present 3-D plots illustrating the outcomes of the sensitivity analysis conducted
in the context of multi-output learning. It is evident from these visual representations that
there is a noticeable disparity in the results when the dataset size is less than 500, signifying
that the model tends to underestimate the training data. As the training dataset size surpasses 500, this visual disparity diminishes, indicating improved performance of the results
with a dataset size exceeding 500.
142

Table D.2: A summary of the Mean Squared Error (MSE) for pressure, CO2 saturation,
vertical displacement, and effective plastic strain of the multi-output U-Net model with well
location input weight of 10 (Case 2 in Table 5.2) under different training data sizes.
MSE Loss
Vertical
Displacement CO2 Saturation Pressure Effective plastic
strain Case Training Size
2.1 100 2.0698E-05 2.7600E-03 4.8505E-04 2.3137E-04
2.2 200 1.8912E-05 1.8316E-03 2.5286E-04 1.1494E-04
2.3 300 7.5814E-06 1.5942E-03 1.9350E-04 8.4369E-05
2.4 400 8.0313E-06 9.8534E-04 6.2673E-05 2.5058E-05
2.5 500 6.1849E-06 8.1404E-04 4.0452E-05 1.6108E-05
2.6 600 5.2279E-06 6.3809E-04 3.0952E-05 1.0066E-05
2.7 700 4.9967E-06 5.4373E-04 2.5499E-05 8.0929E-06
2.8 800 9.5344E-06 4.5889E-04 2.5224E-05 7.2994E-06
2.9 900 8.3585E-06 3.7076E-04 2.3680E-05 7.0834E-06
2.10 1000 4.5016E-06 3.6267E-04 3.0268E-05 1.0385E-05
0.00E+00
5.00E-04
1.00E-03
1.50E-03
2.00E-03
2.50E-03
3.00E-03
0.0E+00
1.0E-04
2.0E-04
3.0E-04
4.0E-04
5.0E-04
6.0E-04
100 300 500 700 900
MSE Loss
Training size
Multi-output Learning Model Sensitivity Analysis
Series2
Series3
Series4
Series1
Effective Plastic Strain (�!"#
$ )
Vertical Displacement (�%)
Pressure
CO2 Saturation
Figure D.2: Multi-output learning model sensitivity analysis on training sample size. The
y-left axis is the MSE loss plotted for pressure, vertical displacement, and effective plastic
strain while the y-right axis is the MSE error plotted for CO2 Saturation.
143

Well location on permeability map
True Observation (Pressure)
Mismatch
datasize=100
datasize=200
datasize=300
datasize=400
datasize=500
datasize=600
datasize=700
datasize=800
datasize=900
datasize=1000
Figure D.3: Sensitivity analysis training results of pressure. Seven examples are chosen for
3-D plotting. The first row plots the well locations on permeability map. The second row
plots the true observation of pressure. The third to the twelfth rows plot the mismatch
between true observation and prediction.
144

Well location on permeability map
True Observation (�!)
Mismatch
datasize=100
datasize=200
datasize=300
datasize=400
datasize=500
datasize=600
datasize=700
datasize=800
datasize=900
datasize=1000
Figure D.4: Sensitivity analysis training results of vertical displacement. Seven examples
are chosen for 3-D plotting. The first row plots the well locations on permeability map.
The second row plots the true observation. The third to the twelfth rows plot the mismatch
between true observation and prediction.
145

Well location on permeability map
True Observation (CO2 saturation)
Mismatch
datasize=100
datasize=200
datasize=300
datasize=400
datasize=500
datasize=600
datasize=700
datasize=800
datasize=900
datasize=1000
Figure D.5: Sensitivity analysis training results of CO2 saturation. Seven examples are
chosen for 3-D plotting. The first row plots the well locations on permeability map. The
second row plots the true observation. The third to the twelfth rows plot the mismatch
between true observation and prediction.
146

Well location on permeability map
True Observation (�!"#
$ )
Mismatch
datasize=100
datasize=200
datasize=300
datasize=400
datasize=500
datasize=600
datasize=700
datasize=800
datasize=900
datasize=1000
Figure D.6: Sensitivity analysis training results of plastic strain. Seven examples are chosen
for 3-D plotting. The first row plots the well locations on permeability map. The second
row plots the true observation. The third to the twelfth rows plot the mismatch between
true observation and prediction.
147

E Surrogate Modeling Data Generation
In this section, we describe the process to generate training dataset for the surrogate model.
The numerical simulation model described in Section 2.4.2 is used for the data generation. A
total of 1000 reservoir simulations are executed, each featuring distinct well locations under
a consistent Bottom Hole Pressure (BHP) control setting. Each simulation run is configured
with four wells. To encompass a broad spectrum of well configurations, we partition the
heterogeneous region into four quadrants, as illustrated in Figure E.1. Each well is then
randomly assigned to one of these quadrants. Furthermore, we maintain a minimum interwell distance of 500 meters, which spans over three grid blocks to ensure an appropriate
separation between the wells.
-4
-3
-2
-1
0
1
Well 1 Well 2
Well 4
Well 3
log(permeability)
Figure E.1: Four quadrants on the heterogeneous permeability region.
148

Abstract (if available)

Abstract Carbon Capture and Storage (CCS) has been recognized as one of the major mitigation strategies for industrial CO2 emission and atmospheric CO2 accumulation. Numerical simulation of fluid flow, transport, and trapping mechanisms have been used to optimize CO2 storage in geologic formations by improving trapping efficiency through injection strategies. Flow models, however, do not capture the geomechanical deformation that can occur during CO2 injection, including reservoir expansion, ground surface uplift, caprock fracturing, and induced seismicity. The coupled flow and geomechanical simulation models are increasingly used to study the geomechanical effects during CO2 injection. An existing gap in the literature is related to the impact of geomechanical effects mentioned above on optimizing the storage performance of geologic formations based on designing the well configuration and control schedules. I present a novel optimization framework for geologic CO2 storage under geomechanical risks, where coupled flow-geomechanics simulations are used to quantify the risks of injection-induced ground surface deformation and rock failure in reservoir and caprock layers. A multi-objective optimization problem is formulated and solved to maximize CO2 storage while minimizing the two forms of geomechanical risks. Caprock failure could result in creating a leakage pathway for CO2 to leak from the storage reservoir to shallower aquifers and even to the atmosphere. To address this concern, I have developed a novel optimization approach aimed at maintaining the caprock integrity during the storage of CO2 in geologic formations under geological uncertainty. The developed workflow integrates advanced numerical optimization algorithms with coupled multiphase flow-geomechanics-fracturing models for simulating the response of the storage reservoir to CO2 injection. Using the geomechanical response of the simulation, I define and quantify the potential caprock failure and CO2 leakage risks. An optimization formulation is used to minimize the risk of caprock fracturing and CO2 leakage by finding the optimal distribution of dynamically changing CO2 injection rates across several wells throughout the injection period. Multiple numerical experiments with increasing complexity are presented to demonstrate the performance of the proposed framework. The flow and geomechanical properties of deep saline aquifers are highly uncertain and can lead to significant variability in the geomechanical responses. To properly account for the geomechanical risks of geologic CO2 storage, a stochastic description of rock properties is incorporated into the optimization framework. Ensemble optimization and reduced sampling method are examined and compared for the stochastic optimization under geologic uncertainties. Simulating CO2 storage under geomechanical risks involves substantial computational costs due to the coupling between multiphase flow and geomechanics. When identification of optimal well locations is considered, the computational burden can escalate hundred-fold, rendering industrial-scale applications of the technology impractical. In this work, we introduce a novel deep-learning framework designed to significantly reduce the computational overhead associated with simulating and quantifying the geomechanical risks of CO2 storage. Our approach leverages deep learning-based surrogate modeling to significantly enhance the efficiency of coupled flow-geomechanics simulations for identifying suitable injection well locations for CO2 storage. We design a novel U-Net convolutional neural network to streamline the process of mapping well location indices combined with permeability to the simulation outputs. Subsequently, we train the U-Net to establish mappings from different well location scenarios to corresponding pressure fields, CO2 saturation and geomechanical outputs, including vertical displacement and plastic strain. A coupled flow-geomechanics CO2 storage simulation model is used to generate the training datasets. The U-Net model is subsequently adopted as an efficient tool to replace the coupled flow-geomechanics simulation needed for identification of injection well locations to minimize geomechanical risks. The trained U-Net model drastically reduces the computational cost by replacing coupled physics simulation with a fast proxy model that can be used to approximately identify the geomechanical risk associated with different well locations.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Dynamics of CO₂ leakage through faults

PDF

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

$A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow$

PDF

A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow

PDF

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

$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

Deep learning for subsurface characterization and forecasting

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

Stochastic oilfield optimization under uncertainty in future development plans

PDF

Optimization of coupled CO₂ sequestration and enhanced oil recovery

PDF

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

PDF

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

PDF

Modeling and simulation of complex recovery processes

$Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells$

PDF

Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells

$History matching for fractured reservoirs with Hough-transform-based parameterization$

PDF

History matching for fractured reservoirs with Hough-transform-based parameterization

PDF

Modeling and simulation of multicomponent mass transfer in tight dual-porosity systems (unconventional)

PDF

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

$Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure$

PDF

Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure

PDF

Stress and deformation analysis on fluid-exposed reservoir rocks

PDF

Inverse modeling and uncertainty quantification of nonlinear flow in porous media models

$An extended finite element method based modeling of hydraulic fracturing$

PDF

An extended finite element method based modeling of hydraulic fracturing

Asset Metadata

Creator Zheng, Fangning (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Degree Conferral Date 2023-12

Publication Date 12/19/2023

Defense Date 12/18/2023

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

Tag carbon capture and storage,carbon storage,convolutional neural network,deep learning,genetic algorithm,geologic CO2 storage,geomechanics,machine learning,OAI-PMH Harvest,optimization,proxy model,reservoir engineering

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jafarpour, Behnam (committee chair), Ghanem, Roger (committee member), Jha, Birendra (committee member)

Creator Email fangning93@gmail.com,fangninz@usc.edu

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

Unique identifier UC113797584

Identifier etd-ZhengFangn-12573.pdf (filename)

Legacy Identifier etd-ZhengFangn-12573

Document Type Dissertation

Format theses (aat)

Rights Zheng, Fangning

Internet Media Type application/pdf

Type texts

Source 20231221-usctheses-batch-1117 (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.

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