suppReferences

# Optimizing Carbon Storage Operations for Long-Term Safety

Yizheng Wang
Earth and Planetary Sciences
Stanford University
Stanford, CA 94305
<EMAIL_ADDRESS>
Markus Zechner
Earth and Planetary Sciences
Stanford University
Stanford, CA 94305
<EMAIL_ADDRESS>
Gege Wen
Energy Science and Engineering
Stanford University
Stanford, CA 94305
<EMAIL_ADDRESS>
Anthony Louis Corso
Aeronautics and Astronautics
Stanford University
Stanford, CA 94305
<EMAIL_ADDRESS>
John Michael Mern
Aeronautics and Astronautics
Stanford University
Stanford, CA 94305
<EMAIL_ADDRESS>
Mykel J. Kochenderfer
Aeronautics and Astronautics
Stanford University
Stanford, CA 94305
<EMAIL_ADDRESS>
Jef Karel Caers
Earth and Planetary Sciences
Stanford University
Stanford, CA 94305
<EMAIL_ADDRESS>

###### Abstract

To combat global warming and mitigate the risks associated with climate
change, carbon capture and storage (CCS) has emerged as a crucial technology.
However, safely sequestering CO2 in geological formations for long-term
storage presents several challenges. In this study, we address these issues by
modeling the decision-making process for carbon storage operations as a
partially observable Markov decision process (POMDP). We solve the POMDP using
belief state planning to optimize injector and monitoring well locations, with
the goal of maximizing stored CO2 while maintaining safety. Empirical results
in simulation demonstrate that our approach is effective in ensuring safe
long-term carbon storage operations. We showcase the flexibility of our
approach by introducing three different monitoring strategies and examining
their impact on decision quality. Additionally, we introduce a neural network
surrogate model for the POMDP decision-making process to handle the complex
dynamics of the multi-phase flow. We also investigate the effects of different
fidelity levels of the surrogate model on decision qualities.

## Main

On our current trajectory, climate models project the global temperature will
be $2.8^{\circ}$C hotter at the end of this century [1], which almost doubles
the $1.5^{\circ}$C target set by the Paris agreement. Failing to curb
greenhouse emissions by transitioning to a low-carbon economy is catastrophic.
A report published by the Intergovernmental Panel on Climate Change (IPCC) in
2018 suggests even a half a degree increase in global temperature from
$1.5^{\circ}$C to $2^{\circ}$C can have significant negative impacts on human
and natural systems [2, 3]. As temperatures continue to rise, the implications
and risks associated with climate change become more severe, affecting a wide
range of areas from flood to drought, with disastrous consequences on food
supply and ecosystems [4, 5, 6].

Many technologies can help to reduce global warming, such as renewable energy
sources, technologies to improve efficiency, nuclear, and carbon capture and
storage (CCS) [7, 8, 9, 10]. Among these technologies, CCS is the only one
that can reduce $\mathrm{CO_{2}}$ emissions from large centralized sources
[11]. Many mitigation portfolios and emission reduction approaches, such as
those assessed by IPCC and International Energy Agency (IEA), include CCS [12,
13, 14]. It is essential for addressing climate change and decarbonizing of
the global energy system.

CCS is designed to capture $\mathrm{CO_{2}}$ emissions from large industrial
sources, such as power plants and industrial facilities, and then transport
and store them underground. The captured $\mathrm{CO_{2}}$ is typically
injected into underground geological formations, such as depleted oil and gas
reservoirs, saline aquifers, and coal seams [15, 16]. Saline aquifers have the
highest $\mathrm{CO_{2}}$ storage potential compared to other geological
storage options, making them important prospects to scale CCS [17].

By 2022, 20 countries are already developing and operating more than 130
$\mathrm{CO_{2}}$ storage sites [18]. However, sequestering $\mathrm{CO_{2}}$
into geological formations for long-term storage comes with many subsurface
technological challenges and safety concerns. Although $\mathrm{CO_{2}}$
sequestration projects may seem similar to other subsurface applications that
involve working with the flow in porous media, such as groundwater or oil and
gas [19, 20], they pose unique challenges. First, less is known about saline
aquifers than oil and gas reservoirs, resulting in significantly higher
subsurface uncertainty. Second, one must account for the complex physical and
chemical processes such as trapping mechanisms, mineral precipitations, and
geomechanics that result from interactions between supercritical
$\mathrm{CO_{2}}$, brine, and rock [21]. Injection of $\mathrm{CO_{2}}$ into
geological formations comes with significant risks, such as induced
earthquakes, fractured cap rocks, and reactivation of faults that can lead to
$\mathrm{CO_{2}}$ leakage, harming human life and the environment [22]. Thus,
any decisions made in the course of a carbon storage project, such as injector
locations, injection rates, and monitoring and mitigation strategies, are
safety critical. Overall, the inherent subsurface uncertainty combined with
complex physio-chemical processes imposes challenging decision problems for
safe and effective carbon storage operations that require tight coupling of
information gathering and sequential planning.

The Partial Observable Markov Decision Process (POMDP) is a framework for
making sequential decisions under uncertainties (see Methods) [23]. The
framework has been successfully applied to many domains, such as autonomous
driving, gaming, and healthcare [24, 25, 26, 27]. Optimization of decision
strategies within this framework has resulted in systems that perform better
than human decision-makers and other state-of-the-art techniques [28].

Figure 1: Optimizing for carbon storage operations.

As shown in Fig 1, optimizing carbon storage operations requires
multidisciplinary knowledge, including computer science, geology, engineering,
economics, and politics. In a general POMDP formulation for carbon storage,
the states may include the geological and engineering characteristics of a
subsurface reservoir, flow properties, and joint space of well types,
controls, and configurations. The actions may include running well tests,
drilling wells, setting well controls, and even abandoning the project. The
observations can come from various sources, including local high-resolution
observations through boreholes and spatial observations such as seismic
signals through active and passive monitoring [29, 30]. The transition
function is defined by the complex simulation of flow and transport in porous
media that govern the injection and migration of the $\mathrm{CO_{2}}$ plume.
The objective is to maximize the amount of sequestrated $\mathrm{CO_{2}}$ in a
subsurface reservoir with some unknown properties and ensure safety for the
long term. In carbon storage operations, it is necessary to consider the vast
number of scenarios that arise from the high-dimensional state, action, and
observation spaces over long planning horizons. This combinatorial problem
presents a significant challenge for humans to find the optimal sequence of
actions, especially when considering the incorporation of present and future
information gathering. We propose an effective way to efficiently search for
the best sequence of actions within this complex decision-making space that
maximizes expected discounted return.

We formulate the carbon storage problem as a POMDP and solve it using
approximate belief state planning. We show the importance of automated
decision-making systems for ensuring safe operations by comparing
algorithmically generated decisions against human judgments. We demonstrate
how reasoning about future information is crucial for safety. In addition, we
exemplify how planning over complex dynamics such as multi-phase flow becomes
feasible when using a neural network surrogate model.

## 3D saline aquifer reservoir

To demonstrate the usage of the decision systems in real-world carbon capture
operations, we present a 3D saline aquifer reservoir. The aquifer is modeled
by a set of engineering parameters and a 3D geological porosity map. We have
$\mathrm{CO_{2}}$ injectors, which inject $\mathrm{CO_{2}}$ for a fixed
period, and monitoring wells, which can be used to observe $\mathrm{CO_{2}}$
saturation. After the injection period, $\mathrm{CO_{2}}$ injection ceases.
However, the $\mathrm{CO_{2}}$ still moves and spreads due to the
gravitational and buoyancy forces for a more extended time, referred to as the
post-injection period. The injected $\mathrm{CO_{2}}$ can be classified into
three categories: exited, free, and trapped. The exited $\mathrm{CO_{2}}$
indicates the injected $\mathrm{CO_{2}}$ has leaked out of the aquifer. The
free $\mathrm{CO_{2}}$ is $\mathrm{CO_{2}}$ that remains in the aquifer but is
not physically bounded and is, therefore, movable. The trapped
$\mathrm{CO_{2}}$ is already in equilibrium and chemically or physically
trapped within the aquifer.

The POMDP formulation for the 3D example is shown in Fig. 2a. The state is
defined by the joint space of well configurations, the amount of trapped,
free, and exited $\mathrm{CO_{2}}$, and the 3D porosity map. The porosity
field remains unknown. The action space is the locations where the wells may
be placed, including both $\mathrm{CO_{2}}$ injectors and monitoring wells.
All wells measure the porosity at the well locations. The monitoring wells
observe a time history of $\mathrm{CO_{2}}$ saturation. All observations
contain normally distributed noise. The transition dynamics are encoded in a
multi-phase flow simulator. It simulates how $\mathrm{CO_{2}}$ migrates and
calculates the mass of each portion of injected $\mathrm{CO_{2}}$ for the
entire injection and post-injection periods. The reward function penalizes the
total amount of $\mathrm{CO_{2}}$ escaping from the aquifer and the free
$\mathrm{CO_{2}}$ and positively rewards the total trapped $\mathrm{CO_{2}}$
mass, where $\Delta m_{\rm exited}$, $\Delta m_{\rm free}$, and $\Delta m_{\rm
trapped}$ are, respectively, the changes in exited, free, and trapped
$\mathrm{CO_{2}}$ between states $s$ and $s^{\prime}$.

(a)

(b)

Figure 2: 3D saline aquifer carbon storage POMDP demonstration and its policy.
(a) POMDP specification for individual components, emphasizing that the
proposed belief state planning approach is a general methodology, without
making specific claims about the superiority of a particular planner or
parameter set. This allows operators and stakeholders to make trade-offs
between various objectives. (b) The policy and overall performance of our
approach in a single representative case, with 3D grids illustrating the
aquifer structure.

## Reducing computational time using a surrogate model

While planning, the belief state planning approach simulates future scenarios
and reasons how future information would affect the current decision. This
process requires a large number of physical simulations. In our work, this
implies our approach needs many simulation episodes of external multi-phase
flow in porous media, which is highly computationally intensive (minutes to
hours per simulation). To make planning feasible under a multi-phase flow
environment, we use a deep learning-based surrogate model as an alternative to
numerical simulation. Once trained, the surrogate model provides significant
speed-up compared to traditional simulation while maintaining comparable
accuracy. Therefore, surrogate models enable efficient repetitive forward
simulations needed by our approach.

We use two surrogate models to predict multi-phase flow responses during the
injection and post-injection period. As demonstrated in Fig. 3, the surrogate
model consists of two components: a Fourier Neural Operator (FNO) [31] that
predicts 4D (3D space-time) CO2 gas saturation distribution, followed by a
convolutional neural network (CNN) that predicts the temporal evolution of
trapped, free, and exited CO2 given 4D gas saturation distribution.

Figure 3: Surrogate model workflow for injection and post-injection period

FNO is a type of neural operator [32] that uses the Fourier integral kernel
operator to learn empirical relationships given input-output data pairs. FNO-
typed architecture has shown superior accuracy and data efficiency for CO2
storage compared to the previous convolutional neural network-based surrogate
model problems [33, 34]. The input data to the injection period FNO model
includes a combination of well locations, drilling strategy, porosity
heterogeneity, rock compressibility, irreducible water/gas saturation,
capillary entry pressure, water viscosity, as well as spatial and temporal
encoding. Since all injection wells are shut down during the post-injection
period, the well location and drilling strategy information are not needed by
the post-injection FNO model. Instead, we provide the injection period gas
saturation map to aid the prediction of the post-injection $\mathrm{CO_{2}}$
migration.

We extend the original FNO architecture into 4D to predict the dynamic change
of 3D gas saturation over 30 years of active injection and 500 years of post-
injection monitoring. As shown in Fig. 3, the 3D spatial domain of the
training data is an irregular geological formation. This leads to complex
interactions between the flow behavior and gravity. The onset of and location
of each injection well further increases the dimensionality of this problem.
Nevertheless, the FNO architecture successfully learns this complex input-
output 4D mapping with minimal overfitting. After obtaining the gas saturation
prediction, we then use a CNN model to map the dynamic 3D gas saturation at
each time step into time series predictions of $\mathrm{CO_{2}}$ mass. The CNN
architecture consists of three channels for trapped, free, and mobile
$\mathrm{CO_{2}}$ masses (see Extended Data Table 3–6 for details of the
surrogate architecture and parameters). All training and test datasets
utilized in this work were generated employing the MATLAB Reservoir Simulation
Toolbox (MRST) [35], a versatile and widely recognized computational tool in
the domain of reservoir simulations. Using the surrogate model, we obtain
predictions with three orders of magnitude speed-up compared to traditional
numerical simulators.

## Results for the 3D case

In this work, we compare our proposed method with two baseline approaches. The
first baseline, referred to as the random policy, involves selecting injectors
randomly. The second baseline, called the expert policy, relies on experts
utilizing their knowledge of the aquifer’s geological characteristics and the
subsurface flow processes to determine injector locations. We evaluate the
performance of the random, expert, and belief state planning policies.

Our approach makes decisions based on its belief at each time step. To
illustrate the belief regarding the true porosity map, we present the variance
maps of 100 sampled porosity maps from the belief after each decision step
(posterior) in Fig. 2b. Darker colored areas signify lower variance, as the
data acquired through drilling a well and observing $\mathrm{CO_{2}}$
saturation over time via the monitoring well lead to greater certainty in the
belief. Fig.4a summarizes the results due to different approaches on a single
representative case. Our approach demonstrates superior performance, resulting
in a 16.3% increase in trapped $\mathrm{CO_{2}}$ and a 98.8% reduction in
leaked $\mathrm{CO_{2}}$ compared to the expert policy at the end of the post-
injection period.

To comprehensively evaluate and compare the three approaches, we test their
policies on ten cases with different 3D porosity maps as ground truths.
Compared to the expert policies, the policies that result from our approach
end up achieving 87% less exited $\mathrm{CO_{2}}$ and 4% more trapped
$\mathrm{CO_{2}}$ (see Fig.4b). On the other hand, the random policies have a
2.9-fold increase in the quantity of exited $\mathrm{CO_{2}}$ when compared to
the expert policies. In addition, we use both mean and standard error of
discounted return as metrics. As a result, our approach attains the highest
expected discounted return across ten cases. A larger expected discounted
return suggests a safer and more effective carbon capture policy. A smaller
standard error of the discount return indicates a more robust performance
around the expected discounted return. The results suggest our approach can
provide robust performance compared to the other two baseline methods,
ensuring safe and effective carbon storage operations.

(a)

(b)

Figure 4: Performance comparison of three solution methods on the 3D example.
(a) $\mathrm{CO_{2}}$ mass profiles at the end of the post-injection period
for a single representative case. (b) Aggregated performance metrics across
ten cases.

## Effects of Surrogate Model Fidelity Levels on Decision Quality

In this work, we use a trained surrogate model as a proxy for the complex
decision environment, raising the question of how planning with different
fidelity levels of the surrogate model might affect the performance of our
proposed approach. To address this question, we trained surrogate models at
two distinct fidelity levels. The high-fidelity model exhibits significantly
smaller root mean squared errors (RMSE) in predicting $\mathrm{CO_{2}}$ mass
over time compared to its low-fidelity counterpart. For instance, the RMSE of
trapped, free, and exited $\mathrm{CO_{2}}$ mass predictions for the high-
fidelity model are 1.1 MT, 1.8 MT, and 2.3 MT, respectively. Conversely, the
corresponding RMSE values for the low-fidelity model are 1.8 MT, 2.8 MT, and
3.6 MT.

When comparing the density plots of errors between high and low fidelity
surrogate model predictions of $\mathrm{CO_{2}}$ mass and MRST outputs at the
end of the post-injection period, we find that the error density plots of the
high fidelity model are more closely centered around 0 and exhibit smaller
variance than those of the low fidelity model (see Fig. 5a). We calculated the
JS-divergences between the density curves, which are 0.0075, 0.0126, and
0.0057, respectively. The fidelity level of the surrogate models used in
planning has a significant impact on decision quality. Fig. 5b shows that when
planning with the high fidelity model, the density curve of the discounted
return moves significantly to the right, indicating safer and more effective
policies than planning with the low fidelity model. The JS-divergence between
the two density plots on the discounted return due to high and low fidelity
models is 42.2549, highlighting the significant impact that even small changes
in the fidelity of the surrogate model can have on overall discounted returns.

In light of these results, the question naturally arises:“How accurate is
accurate enough?" While it is clear that model accuracy plays a crucial role
in the quality of the decision-making process, it is challenging to determine
an exact threshold for sufficient accuracy. Practically, perfect accuracy is
often unattainable, and operators must strike a balance between computational
cost and model fidelity. The key takeaway from our analysis is that accuracy
matters, but the surrogate model does not have to be perfect. We recommend
that, as a best practice, any plan derived from a surrogate model should be
validated using a higher fidelity model before implementation. This ensures
that the most critical decisions are based on the most accurate information
available, while still benefiting from the computational efficiency provided
by surrogate models during the planning process.

(a)

(b)

Figure 5: Comparative analysis of surrogate model fidelity. (a) Density plots
illustrating the errors between the high and low fidelity surrogate
predictions and MRST outputs. (b) Density plot on discounted return for high
and low fidelity surrogate models across ten cases.

## Demonstrating the value of the monitoring well

As shown in Fig. 2, our approach needs to make sequential decisions on the
locations of both the monitoring well and injectors. While it is common
practice to include a monitoring well in the field to obtain subsurface
information and monitor $\mathrm{CO_{2}}$ saturation, our study, for the first
time, quantifies the utility of the monitoring well in this context. Adding
the monitoring well increases the complexity of making sequential decisions
because of the alignment of the monitoring well with injectors to maximize
information gain and better inform future actions. To evaluate the value and
necessity of the monitoring well, we use the same POMDP formulation except
without the monitoring well in the action space, resulting in a different
observation space.

Table 1: Our proposed approach performance comparison under different
observation strategies

Plan scenarios | Discounted return |  $\mathrm{CO_{2}}$ mass profile (MT)a
---|---|---
(mean$\pm$SE) | Trapped | Free | Exited
No monitoring well | $-1262.9\pm 485.3$ | 30.5 | 26.9 | 1.6
Monitoring well | $-593.8\pm 367.7$ | $\mathbf{30.9}$ | 27.2 | 0.9
4D seismic | $\mathbf{-475.8}\pm\mathbf{157.8}$ | 30.4 | 27.8 | $\mathbf{0.8}$

11footnotetext: The values presented in the $\mathrm{CO_{2}}$ mass profile
represent the data collected at the conclusion of the post-injection phase.

As shown in the first two rows of Table 1, planning with the monitoring well
has a higher expected discounted reward and a lower standard error across ten
cases. This suggests that the extra information from the monitoring well can
enhance the decision quality.

## Extending the point-wise observations to spatial observations

One of the key challenges in safe carbon storage operations is the lack of
knowledge on $\mathrm{CO_{2}}$ plume migration under complex geological
regimes while planning. Geophysical monitoring (e.g., through seismic surveys)
can provide predictions on $\mathrm{CO_{2}}$ plume saturation and delineate
its boundaries [36]. To demonstrate the possibility of adding seismic
observations to the carbon storage POMDP, we replace placing the monitoring
well with the ability to make seismic observations before drilling each
injector in the action space of the 3D example. We mimic the seismic
saturation observations by adding spatial filtering with Gaussian kernel [37]
over the underlying true saturation map (see Extended Data Fig. 1 for an
example). Among all the observation strategies outlined in Table 1, employing
seismic observations in the planning process results in the safest policies
and yields the highest discounted return and the lowest standard errors
associated with it.

## Discussion and conclusion

In this study, we address a sequential decision problem in carbon storage by
employing belief state planning to ensure safe and effective operations. Our
method significantly outperforms expert decisions in terms of expected
discounted return and standard error across ten cases. For instance, in the
representative 3D case, our proposed approach reduces the amount of leaked
$\mathrm{CO_{2}}$ by over 87% when compared to the expert policy. Furthermore,
we find that incorporating information from monitoring wells and seismic
observations allows our approach to achieving safer and more effective
policies.

To make planning under a multi-phase flow environment feasible, we train an
FNO neural network as a surrogate flow simulator. We investigate the impact of
different fidelity levels of the surrogate model on the decision quality of
our proposed approach. Our results indicate that planning with a high-fidelity
model leads to substantially safer and more effective decisions for carbon
storage operations. Large-scale simulations serve as an effective means of
evaluating our methodology, as physical experiments would not only take years
to complete but also fail to adequately assess the robustness of the decision
strategies. Moreover, it would be challenging to fairly compare various
methodologies across a broad range of scenarios using physical experiments
alone.

As the development of carbon storage projects gains momentum globally,
ensuring the safe operation of these initiatives is of paramount importance.
Our research highlights the potential necessity of our proposed approach in
ensuring the long-term safety of carbon storage projects.

## Data availability

The data used for this work is made available in the repository
https://github.com/yizhengw/OptimizingCarbonStorage.

## Code availability

The code used to produce the results for this work is made available in the
repository https://github.com/yizhengw/OptimizingCarbonStorage.

## References

  * [1] United Nations Environment Programme (2022). Emissions Gap Report 2022: The closing window — climate crisis calls for rapid transformation of societies. Nairobi. See https://www.unep.org/emissions-gap-report-2022/ (accessed 10 December 2022).
  * [2] IPCC. 2018 Global Warming of 1.5°C. An IPCC special report on the impacts of global warming of 1.5°C above pre-industrial levels and related global greenhouse gas emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty. See https://www.ipcc.ch/sr15/ (accessed 10 December 2022).
  * [3] Christian LE Franzke.  Nonlinear climate change.  Nature Climate Change, 4(6):423–424, 2014.
  * [4] Robert WR Parker, Julia L Blanchard, Caleb Gardner, Bridget S Green, Klaas Hartmann, Peter H Tyedmers, and Reg A Watson.  Fuel use and greenhouse gas emissions of world fisheries.  Nature Climate Change, 8(4):333–337, 2018.
  * [5] Kevin E Trenberth, Aiguo Dai, Gerard Van Der Schrier, Philip D Jones, Jonathan Barichivich, Keith R Briffa, and Justin Sheffield.  Global warming and changes in drought.  Nature Climate Change, 4(1):17–22, 2014.
  * [6] Jon Norberg, Mark C Urban, Mark Vellend, Christopher A Klausmeier, and Nicolas Loeuille.  Eco-evolutionary responses of biodiversity to climate change.  Nature Climate Change, 2(10):747–751, 2012.
  * [7] David EHJ Gernaat, Harmen Sytze de Boer, Vassilis Daioglou, Seleshi G Yalew, Christoph Müller, and Detlef P van Vuuren.  Climate change impacts on renewable energy supply.  Nature Climate Change, 11(2):119–125, 2021.
  * [8] Shahrouz Abolhosseini, Almas Heshmati, and Jorn Altmann.  A review of renewable energy supply and energy efficiency technologies.  2014\.
  * [9] Diana Silva Siqueira, Josué de Almeida Meystre, Maicon Queiroz Hilário, Danilo Henrique Donato Rocha, Genésio José Menon, and Rogério José da Silva.  Current perspectives on nuclear energy as a global climate change mitigation option.  Mitigation and Adaptation Strategies for Global Change, 24(5):749–777, 2019.
  * [10] Ronald Wennersten, Qie Sun, and Hailong Li.  The future potential for carbon capture and storage in climate change mitigation–an overview from perspectives of technology, economy and risk.  Journal of Cleaner Production, 103:724–736, 2015.
  * [11] Milenka Mitrović and Alexandra Malone.  Carbon capture and storage (CCS) demonstration projects in canada.  Energy Procedia, 4:5685–5691, 2011.
  * [12] Bert Metz, Ogunlade Davidson, HC De Coninck, Manuela Loos, and Leo Meyer.  IPCC special report on carbon dioxide capture and storage.  2005\.
  * [13] International Energy Agency. 2020 CCUS in the transition to net-zero emissions. See https://www.iea.org/reports/ccus-in-clean-energy-transitions/ccus-in-the-transition-to-net-zero-emissions (accessed 10 December 2022).
  * [14] Glen P Peters, Robbie M Andrew, Tom Boden, Josep G Canadell, Philippe Ciais, Corinne Le Quéré, Gregg Marland, Michael R Raupach, and Charlie Wilson.  The challenge to keep global warming below 2 c.  Nature Climate Change, 3(1):4–6, 2013.
  * [15] Mouadh Addassi, Abdirizak Omar, Hussein Hoteit, Abdulkader M Afifi, Serguey Arkadakskiy, Zeyad T Ahmed, Noushad Kunnummal, Sigurdur R Gislason, and Eric H Oelkers.  Assessing the potential of solubility trapping in unconfined aquifers for subsurface carbon storage.  Scientific Reports, 12(1):1–15, 2022.
  * [16] Zaixing Huang, Christine Sednek, Michael A Urynowicz, Hongguang Guo, Qiurong Wang, Paul Fallgren, Song Jin, Yan Jin, Uche Igwe, and Shengpin Li.  Low carbon renewable natural gas production from coalbeds and implications for carbon capture and storage.  Nature Communications, 8(1):1–11, 2017.
  * [17] Karsten Michael, Alexandra Golab, Valeriya Shulakova, Jonathan Ennis-King, Guy Allinson, Sandeep Sharma, and Toby Aiken.  Geological storage of CO2 in saline aquifers—a review of the experience from existing storage operations.  International Journal of Greenhouse Gas Control, 4(4):659–667, 2010\.
  * [18] IEA (2022), CO2 storage resources and their development, IEA, Paris. https://www.iea.org/reports/co2-storage-resources-and-their-development (accessed 12 December 2022).
  * [19] Céline Scheidt, Lewis Li, and Jef Caers.  Quantifying uncertainty in subsurface systems.  John Wiley & Sons, 236, 2018.
  * [20] Daniel M Tartakovsky.  Assessment and management of risk in subsurface hydrology: A review and perspective.  Advances in Water Resources, 51:247–260, 2013.
  * [21] MA Celia, S Bachu, JM Nordbotten, and KW 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.
  * [22] Sandra Ó Snæbjörnsdóttir, Bergur Sigfússon, Chiara Marieni, David Goldberg, Sigurður R Gislason, and Eric H Oelkers.  Carbon dioxide storage through mineral carbonation.  Nature Reviews Earth & Environment, 1(2):90–102, 2020.
  * [23] Mykel J Kochenderfer, Tim A Wheeler, and Kyle H Wray.  Algorithms for decision making.  MIT Press, 2022.
  * [24] David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, et al.  Mastering the game of go without human knowledge.  Nature, 550(7676):354–359, 2017.
  * [25] Xuanchen Xiang and Simon Foo.  Recent advances in deep reinforcement learning applications for solving partially observable Markov decision processes (POMDP) problems: Part 1—fundamentals and applications in games, robotics and natural language processing.  Machine Learning and Knowledge Extraction, 3(3):554–581, 2021.
  * [26] Wenqian Zhang and Haiyan Wang.  Diagnostic policies optimization for chronic diseases based on POMDP model.  In Healthcare, volume 10, page 283. MDPI, 2022.
  * [27] Constantin Hubmann, Nils Quetschlich, Jens Schulz, Julian Bernhard, Daniel Althoff, and Christoph Stiller.  A pomdp maneuver planner for occlusions in urban scenarios.  In 2019 IEEE Intelligent Vehicles Symposium (IV), pages 2172–2179. IEEE, 2019.
  * [28] Yizheng Wang, Markus Zechner, John Michael Mern, Mykel J Kochenderfer, and Jef Karel Caers.  A sequential decision-making framework with uncertainty quantification for groundwater management.  Advances in Water Resources, 166:104266, 2022.
  * [29] Thomas M Daley, Ray D Solbau, Jonathan B Ajo-Franklin, and Sally M Benson.  Continuous active-source seismic monitoring of CO2 injection in a brine aquifer.  Geophysics, 72(5):A57–A61, 2007.
  * [30] James P Verdon, J-Michael Kendall, Don J White, Doug A Angus, Quentin J Fisher, and Ted Urbancic.  Passive seismic monitoring of carbon dioxide storage at weyburn.  The Leading Edge, 29(2):200–206, 2010.
  * [31] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier neural operator for parametric partial differential equations.  arXiv preprint arXiv:2010.08895, 2020.
  * [32] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar.  Multipole graph neural operator for parametric partial differential equations.  Advances in Neural Information Processing Systems, 33:6755–6766, 2020.
  * [33] 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.
  * [34] Gege Wen, Zongyi Li, Qirui Long, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M Benson.  Accelerating carbon capture and storage modeling using fourier neural operators.  arXiv preprint arXiv:2210.17051, 2022.
  * [35] Knut-Andreas Lie.  An introduction to reservoir simulation using matlab/gnu octave: User guide for the matlab reservoir simulation toolbox (mrst).  2019\.
  * [36] Manzar Fawad and Nazmul Haque Mondol.  Monitoring geological storage of CO2: A new approach.  Scientific Reports, 11(1):1–9, 2021.
  * [37] Chanda Bhabatosh et al.  Digital image processing and analysis.  PHI Learning Pvt. Ltd., 1977.

## Methods

### POMDPs

A POMDP is a model for sequential decision making under uncertainty. A POMDP
is defined as a tuple ($\mathcal{S},\mathcal{A},\mathcal{O},O,T,R,\gamma$),
where $\mathcal{S}$ is the set of states, $\mathcal{A}$ is the set of actions,
$\mathcal{O}$ is the set of observations, $O$ is the observation function, $T$
is the state transition function, $R$ is the reward function, and $\gamma$ is
the discount factor. At each step, an agent in state $s\in S$ takes action
$a\in\mathcal{A}$, transitions to the next state $s^{\prime}\in\mathcal{S}$
with probability $T(s^{\prime}\mid s,a)$ and receives both an observation
$o\in\mathcal{O}$ with probability $O(o\mid s^{\prime},a)$ and a reward
$r=R(s,a,s^{\prime})$. In a POMDP, the underlying state of the environment is
not fully known. Instead, the agent maintains a distribution over states
called a belief $b(s)$. The agent infers its belief from the history of
observations and actions according to Bayes’s rule,
$b^{\prime}(s^{\prime})=p(s^{\prime}\mid b,a,o)\propto O(o\mid
s^{\prime},a)\sum_{s\in\mathcal{S}}T(s^{\prime}\mid s,a)b(s)$.

A policy is a function that maps beliefs to actions. The goal is to find the
optimal policy $\pi^{*}$ that maximizes the expected sum of discounted reward,
$\sum_{t=1}^{n}\gamma^{n-1}r_{t}$. For the infinite horizon sequential problem
where $n\to\infty$, the discounted factor $\gamma\in[0,1)$ is introduced to
bound the total discounted reward.

In general, there are two approaches to generating a policy for POMDPs:
offline or online during action execution. Offline solvers perform most of the
computation prior to interacting with the environment while online solvers
perform computation between state transitions.

### 3D saline aquifer example

The 3D aquifer grids have an extension of eighty grid blocks in the x- and the
y-direction (lateral extension) and eight blocks in the z-direction (depth).
There are three injectors and one monitoring well. The monitoring well needs
to be placed before any injectors. The three injectors are placed in a time
sequence at the 1st, 11th, and 21st years. The injection period and post-
injection period, respectively, are 30 and 500 years resulting in a total of
530 years of simulation time. We assume permeability $K$ follows a
deterministic function of porosity $\phi$, with
$K=\frac{10^{-10}\phi^{3}}{58.32(1-\phi)^{2}}$. Therefore, our approach only
needs to keep track of the porosity map in its state space and observes
porosity values at each well location (see Extended Data Table 1-2 for
specific values of geological, engineering, POMDP, and planner parameters).

### Surrogate model training details

For both the injection and post-injection surrogate models, we train the
saturation and mass separately. We trained the gas situation FNO models using
300 training samples and 100 test samples for around 300 epochs. Each epoch
takes approximately 3 mins. For the $\mathrm{CO_{2}}$ mass CNN models, we
trained with 600 training samples and 100 test samples for around 300 epochs.
Each epoch takes approximately 11 s. $\mathrm{CO_{2}}$ mass requires more
training samples because it is more prone to overfitting. The relative
squared-error loss is used for both scenarios to ensure effective gradient
propagation. After training, the average gas saturation error given unseen
inputs is 0.013 during the injection period and 0.016 during the post-
injection period. The average mass error (defined as the mean absolute error
divided by the mass norm) is 0.012 during the injection period and 0.015
during the post-injection period. To obtain surrogate models with two distinct
fidelity levels, we utilize 12,000 training data for the high-fidelity model
and 1,700 training data for the low-fidelity model.

### Partially Observable Monte Carlo Planning with Observation Widening

Due to the relatively long time scales expected in CCS, we use an online
solver called Partially Observable Monte Carlo Planning with Observation
Widening (POMCPOW). POMCPOW is a tree search algorithm used for online
decision-making, which can effectively handle continuous state, action, and
observation spaces sunberg2018online. Prior to POMCPOW, state-of-the-art
online planers such as POMCP and DESPOT struggled to solve POMDPs with
continuous observation spaces somani2013despot,silver2010monte. Their planning
trees suffer from belief collapse beyond the root node and end up generating
suboptimal policies. POMCPOW, on the other hand, creates weighted particle-
based search trees, resulting in richer belief representations at each node,
allowing the search trees to grow deeper. Richer beliefs also make generated
policies contain information-gathering actions and thus yield better
strategies for the problems where information-gathering is needed.

### Ensemble Smoother with Multiple Data Assimilation

Ensemble Smoother with Multiple Data Assimilation (ES-MDA) is a commonly used
ensemble-based inference method for estimating state uncertainties
emerick2013ensemble. Unlike Bayesian inversion using Markov chain Monte Carlo
(MCMC) methods to obtain posterior distributions, ensemble-based methods can
significantly reduce the computational cost of belief updates. Among the
ensemble-based methods, the ensemble Kalman filter (EnKF) is one of the most
popular approaches evensen2003ensemble. However, the recurrent simulation
restarts required in the EnKF sequential data simulation process make this
method impractical for problems that involve expensive forward simulations. On
the other hand, in ensemble smoother (ES), all data are assimilated
simultaneously in a single update skjervheim2011ensemble. Therefore, there is
no need for simulation restarts, which makes it more attractive for practical
applications. However, because ES computes a single update, it may not yield a
good result that matches the data. In order to fix this problem, ES-MDA
introduces an inflation factor ($\alpha$) in front of the data-error
covariance matrix $C_{D}$. This factor helps the algorithm to correctly sample
from the posterior distribution emerick2016analysis.

### Belief update

The porosity map is partially observable. However, the belief over the
porosity map can be informed by both direct measurements of aquifer porosity
at each well location and $\mathrm{CO_{2}}$ saturation history. Because
observations are obtained from two different sources, the belief update is
formulated as follows:

$p(\mathbf{\phi}\mid
d_{1},\mathbf{d_{2}})=\sum_{\mathbf{\phi_{cond}}}p(\mathbf{\phi}\mid\mathbf{\phi_{cond}},\mathbf{d_{2}})p(\mathbf{\phi_{cond}}\mid
d_{1})$ (1)

where $\mathbf{\phi}$ stands for the porosity maps, $d_{1}$ is the observed
porosity value, and $\mathbf{d_{2}}$ is the observed $\mathrm{CO_{2}}$
saturation history at the monitoring well. The term $\mathbf{\phi_{cond}}$
represents the conditional porosity maps with only $d_{1}$ using conditional
Gaussian simulation hoshiya1995kriging. Therefore, we can assume
$\mathbf{\phi}$ is conditional independent with $d_{1}$, statistically
represented as $p(\mathbf{\phi}\mid
d_{1},\mathbf{d_{2}},\mathbf{\phi_{cond}})=p(\mathbf{\phi}\mid\mathbf{d_{2}},\mathbf{\phi_{cond}})$.

To reduce the computational cost of solving the inversion problem with high
dimensional data, we deploy ES-MDA to sample from the posterior of
$\mathbf{\phi}$ given both $\mathbf{d_{2}}$ and $\mathbf{\phi_{cond}}$.

unsrt references

Extended Data Fig. 1: Seismic observation example for a given 3D saturation map Extended Data Table 1: Geological and engineering parameters for the 3D aquifer Variable Description | Value
---|---
Porosity map mean | 0.2
Porosity map variogram sill | 0.001
Porosity map variogram range | 20.0
Porosity map variogram nugget | 0.0001
Rock compressibility | $4.35\text{\times}{10}^{-5}$
Irreducible water saturation | 0.27
Irreducible gas saturation | 0.20
Capillary entry pressure | 5.0
Water viscosity | $8\text{\times}{10}^{-4}$
Extended Data Table 2: POMDP and POMCPOW hyperparameters Variable Description | Value
---|---
Discount Factor ($\gamma$) | 0.99
$\lambda_{\rm exited}$ | -1000.0
$\lambda_{\rm free}$ | -1.0
$\lambda_{\rm trapped}$ | 10.0
Gaussian filtering kernel $\sigma$ | 3.0
Distribution of observation noise | $\mathcal{N}(0,0.1\times\phi_{obs})$
Exploration Coefficient for UCB | 20.0
Observation widening exponent ($\alpha_{\rm obs})$ | 0.7
Observation widening coefficient ($k_{\rm obs})$ | 5.0
Action widening exponent ($\alpha_{\rm act}$) | 0.7
Action widening coefficient ($k_{\rm act})$ | 2.0
Number of tree queries ($N_{\rm query}$) | 100
Extended Data Table 3: Model parameters for gas saturation during the injection period. The Padding denotes a padding operator that accommodates the non-periodic boundaries; Linear denotes the linear transformation to lift the input to the high dimensional space and the projection back to the original space; Fourier4d denotes the 4D Fourier operator; Conv1d denotes the bias term; Add operation adds the outputs together; GELU denotes a Gaussian Error Linear Units layer.  Layer | Operation | Output Shape
---|---|---
Input | - | (30, 80, 80, 8, 11)
Padding | Padding (4) | (38, 88, 88, 16, 11)
Lifting | Linear | (38, 88, 88, 16, 28)
Fourier 1 | Fourier4d/Conv1d/Add/GELU | (38, 88, 88, 16, 28)
Fourier 2 | Fourier4d/Conv1d/Add/GELU | (38, 88, 88, 16, 28)
Fourier 3 | Fourier4d/Conv1d/Add/GELU | (38, 88, 88, 16, 28)
Fourier 4 | Fourier4d/Conv1d/Add | (30, 80, 80, 8, 28)
De-padding | Depadding (4) | (30, 80, 80, 8, 28)
Projection 1 | Linear | (30, 80, 80, 8, 64)
Projection 2 | Linear | (30, 80, 80, 8, 1)
Extended Data Table 4: Model parameters for CO2 mass during the injection period. After the flattening layer, the model splits into three channels where each channel contains an FC1 and FC2. Layer | Operation | Output Shape
---|---|---
Input | - | (30, 80, 80, 8)
Conv1 | Conv3D/BN/ReLU | (30, 40, 40, 4)
Conv2 | Conv3D/BN/ReLU | (30, 20, 20, 2)
Flatten | - | (1, 24000)
FC1 | FC/BN/ReLU | (1, 100)
FC2 | FC/BN/ReLU | (1, 30)
Extended Data Table 5: Model parameters for gas saturation during post-injection period. The Padding denotes a padding operator that accommodates the non-periodic boundaries; Linear denotes the linear transformation to lift the input to the high dimensional space and the projection back to the original space; Fourier4d denotes the 4D Fourier operator; Conv1d denotes the bias term; Add operation adds the outputs together; GELU denotes a Gaussian Error Linear Units layer.  Layer | Operation | Output Shape
---|---|---
Input | - | (20, 80, 80, 8, 11)
Padding | Padding (4) | (28, 88, 88, 16, 11)
Lifting | Linear | (28, 88, 88, 16, 28)
Fourier 1 | Fourier4d/Conv1d/Add/GELU | (28, 88, 88, 16, 28)
Fourier 2 | Fourier4d/Conv1d/Add/GELU | (28, 88, 88, 16, 28)
Fourier 3 | Fourier4d/Conv1d/Add/GELU | (28, 88, 88, 16, 28)
Fourier 4 | Fourier4d/Conv1d/Add | (20, 80, 80, 8, 28)
De-padding | Depadding (4) | (20, 80, 80, 8, 28)
Projection 1 | Linear | (20, 80, 80, 8, 64)
Projection 2 | Linear | (20, 80, 80, 8, 1)
Extended Data Table 6: Model parameters for CO2 mass during the post-injection period. After the flattening layer, the model splits into three channels where each channel contains an FC1 and FC2. Layer | Operation | Output Shape
---|---|---
Input | - | (20, 80, 80, 8)
Conv1 | Conv3D/BN/ReLU | (20, 40, 40, 4)
Conv2 | Conv3D/BN/ReLU | (20, 20, 20, 2)
Flatten | - | (1, 16000)
FC1 | FC/BN/ReLU | (1, 100)
FC2 | FC/BN/ReLU | (1, 20)