Learning CO_2 plume migration in faulted reservoirs with Graph Neural Networks

label1,label3]Xin Jucor1
[label1]Department of Energy Science and Engineering, Stanford University, Stanford, CA, 94305

[label3]TotalEnergies EP Research and Technology USA, Houston, Texas, 77002

[cor1]Corresponding author

label3]François P. Hamon

label1]Gege Wen

label1]Rayan Kanfar

label3]Mauricio Araya-Polo

label1]Hamdi A. Tchelepi

Deep-learning-based surrogate models provide an efficient complement to numerical simulations for subsurface flow problems such as CO$_2$ geological storage.

Accurately capturing the impact of faults on CO$_2$ plume migration remains a challenge for many existing deep learning surrogate models based on Convolutional Neural Networks (CNNs) or Neural Operators.

We address this challenge with a graph-based neural model leveraging recent developments in the field of Graph Neural Networks (GNNs).

Our model combines graph-based convolution Long-Short-Term-Memory (GConvLSTM) with a one-step GNN model, MeshGraphNet (MGN), to operate on complex unstructured meshes and limit temporal error accumulation.

We demonstrate that our approach can accurately predict the temporal evolution of gas saturation and pore pressure in a synthetic reservoir with impermeable faults.

Our results exhibit a better accuracy and a reduced temporal error accumulation compared to the standard MGN model.

We also show the excellent generalizability of our algorithm to mesh configurations, boundary conditions, and heterogeneous permeability fields not included in the training set.

This work highlights the potential of GNN-based methods to accurately and rapidly model subsurface flow with complex faults and fractures.

Graph neural networks, carbon storage, two-phase flow, surrogate model, deep learning

§ INTRODUCTION

In geological carbon storage (GCS), large amounts of supercritical $\mathrm{CO}_{2}$ are injected into subsurface reservoirs for permanent storage and must be monitored over very long periods of time to ensure safe and effective storage [Agency, 2020, Pacala and Socolow, 2004].

Underground formations often exhibit high degrees of heterogeneity characterized by stratigraphic layering [Mallison et al., 2014] and the presence of faults and fractures [Rinaldi and Rutqvist, 2013, Morris et al., 2011].

These geological features critically impact the migration and the long-term behavior of $\mathrm{CO}_{2}$ plumes. Moreover, faults and fractures could potentially lead to hazards such as induced seismicity or leakage [Castelletto et al., 2013, Cusini et al., 2022].

To ensure the safety and effectiveness of $\mathrm{CO}_{2}$ injection projects, faults and fractures must be accurately modeled in high-fidelity (HF) numerical simulations.

As a result, there is a strong interest in building numerical models based on unstructured polyhedral meshes that conform to the complex geological features of the porous medium [Jenny et al., 2002, Mlacnik et al., 2006, Meng et al., 2018].

In addition, managing geological uncertainty in large-scale storage operations requires running a large number of accurate numerical predictions of $\mathrm{CO}_{2}$ plume migration over decades to centuries [Zhu and Zabaras, 2018, Kitanidis, 2015, Strandli et al., 2014]. This leads to extremely high computational costs for planning CO$_2$ storage projects with faults.

Data-driven deep-learning-based (DL) surrogate models for subsurface flow problems have shown great potential to complement HF numerical simulations and reduce the computational burden of uncertainty quantification studies.

Data-driven deep-learning (DL) models rely on data to learn the underlying physics. They approximate the input and output of interest by building statistical models from simulation data generated by HF simulators.

These methods target the minimization of the data loss between prediction fields and label data and can efficiently achieve converged solutions with satisfactory accuracy [Bischof and Kraus, 2021].

Previous DL-based models have shown excellent accuracy in predicting flow dynamics [Zhu and Zabaras, 2018] and better computational efficiency than HF reservoir simulators [Tang et al., 2020, Tang et al., 2021].

Mo et al. developed a DL surrogate model that integrates an autoregressive model with a convolutional-NNs-based (CNNs) encoder-decoder network to forecast CO$_2$ plume migration in random 2D permeability fields [Mo et al., 2019].

Tang et al. [Tang et al., 2020, Tang et al., 2021] combined a residual U-Net (R-U-Net) with convLSTM networks to predict the temporal evolution of the saturation and pressure fields in 2D and 3D oil production simulations.

Their recurrent R-U-Net model was later applied to CO$_2$ storage with coupled flow and geomechanics [Tang et al., 2022].

Wen et al. [Wen et al., 2021] developed an R-U-Net-based surrogate model for CO$_2$ plume migration, where injection durations, injection rates, and injection locations are encoded as channels of input images.

Recently, Wen et al. [Wen et al., 2022, Wen et al., 2023] combined U-Net and Fourier neural operator [Wen et al., 2022] by adding convolutional information in the Fourier layer, which yields significantly improved cost-accuracy trade-off.

However, these existing data-driven surrogate models are limited to Cartesian meshes with simple geometries, which fails to predict in unstructured meshes with stencils that vary in size and shape <cit.>.

For example, CNNs are designed for image processing and exploit the Cartesian structure of pixels, implying that these models can only efficiently operate on regular grids <cit.>.

This limitation significantly undermines the applicability of these surrogates to field-scale models with complex geological features such as faults and fractures.

To overcome these limitations, here we aim to construct a DL surrogate model based on a graph neural network (GNN) that can capture the flow dynamics in realistic subsurface flow problems modeled with complex unstructured meshes conforming to faults and fractures.

The key feature of GNNs is to represent unstructured simulation data as graphs with nodes and edges, in which the nodes represent cell-centered data (e.g., pore pressure, phase saturation). In contrast, the edges represent cell-to-cell connectivity and face-centered data (e.g., transmissibility, Darcy fluxes).

This is key to enabling the DL surrogate to operate with unstructured-mesh-based simulation data containing complex internal structures.

Recently, a class of GNNs named message-passing neural networks (MPNNs) has demonstrated its efficiency in learning forward dynamics [Sanchez-Gonzalez et al., 2020, Pfaff et al., 2020, Lino et al., 2021, Shi et al., 2022].

In MPNNs, a learnable message-passing function is designed to propagate information over the graph through a local aggregation process [Shukla et al., 2022].

Using the local aggregation process works effectively on embedding spaces and helps the model learn better representations. This aggregation operation is pivotal to enable a node to incorporate information from its neighbors, enriching the own representation of the node in the embedding space. As a result, these improved representations produce more accurate predictions and contribute to the potency of MPNNs in handling unstructured graph-based data.

Of particular relevance to this work is the MPNN methodology named MeshGraphNet (MGN) proposed by Pfaff et al. (2020) [Pfaff et al., 2020], in which the training graph is constructed from an HF simulation mesh.

The authors demonstrated that by encoding various physical quantities–depending on underlying physics–as edge features of a graph, MGN could be a fast surrogate trained from unstructured HF simulations.

This work also shows that MGN-based surrogate models have the ability to generalize to meshes unseen during training and can capture internal boundaries more accurately than CNN-based models.

Wu et al. [Wu et al., 2022] applied the MGN architecture to an oil-and-gas problem and developed a hybrid architecture to learn the dynamics of reservoir simulations on Cartesian meshes. The hybrid architecture used a U-Net (a variant of CNNs [Ronneberger et al., 2015]) for the pore pressure and MGN for the phase saturation.

Notably, these GNN surrogate models all focus on next-state predictions, i.e., they approximate the next state of a physical system from the current state and advance in time in an autoregressive manner.

However, recent work has shown that next-state models are prone to suffer from substantial temporal error accumulation [Pfaff et al., 2020, Han et al., 2022] when autoregressively rolling out for a long period.

This limitation is problematic to predict $\mathrm{CO}_{2}$ plume behavior, as $\mathrm{CO}_{2}$ operations often require the simulation of multiple decades of injection for a commercial-scale project.

Therefore, in this work, we introduce a graph-based recurrent neural network to mitigate temporal error accumulation and achieve a reliable long-term prediction.

The proposed recurrent GNN framework includes (1) a modified MGN that encodes and processes the current physical state into embedding spaces of the entire graph and (2) a graph-based recurrent convLSTM (GConvLSTM) [Seo et al., 2018] model that predicts the next state based on the embeddings computed by MGN and on recurrent memories from past states.

In comparison to the original next-state MGN predictor, the proposed algorithm, referred to as MeshGraphNet-Long Short-Term Memory (MGN-LSTM), can better mitigate temporal error accumulation and significantly improves the performance for the long-term prediction of CO$_2$ plume behavior.

Our MGN-LSTM algorithm can accurately approximate HF simulations on unstructured meshes and is generalizable to meshes, boundary conditions, and permeability fields unseen during training.

Our main contributions include: (1) using GNN to overcome the current limitations of previous surrogate models to handle complex simulation data on unstructured meshes in the context of CO$_2$ geological storage in faulted reservoirs; (2) introducing the accurate MGN-LSTM architecture reducing temporal error accumulation; and (3) demonstrating the generalizability of MGN-LSTM to unseen meshes and boundary conditions as well as its stable extrapolation to future states.

This paper proceeds as follows.

In Section <ref>, we introduce the two-phase flow equations applicable to CO$_2$ geological storage.

Section <ref> describes the proposed surrogate model (MGN-LSTM) and associated data-processing and training procedures.

In Section <ref>, the MGN-LSTM is used to predict saturation and pore pressure fields in two-phase flow (CO$_2$-brine) problems.

Section <ref> demonstrates the improved accuracy of MGN-LSTM compared to the standard MGN algorithm.

Section <ref> concludes the work and suggests future research directions.

§ PROBLEM STATEMENT

§.§ Governing equations of CO$_2$-brine flow

In this work, we consider miscible two-phase (gas and aqueous) two-component (H$_2$O and CO$_2$) flow in a compressible porous medium.

We employ a 2D domain in the $x-y$ plane for simplicity, but the model and algorithms presented here can be extended to 3D.

The H$_2$O component is only present in the aqueous phase, while the CO$_2$ component can be present in both the aqueous and the gas phases.

We denote the aqueous and the gas phases using the subscripts $a$ and $g$, respectively.

The mass conservation of each component reads:

\begin{equation}
\frac{\partial}{\partial t} \left( \phi \sum_{\ell = 1}^{2} x_{c\ell} \rho_\ell s_\ell \right) + \nabla \cdot \left( \sum_{\ell = 1}^{2} x_{c\ell} \rho_\ell \boldsymbol{v}_\ell \right) +  \sum_{\ell = 1}^{2} x_{c\ell} \rho_\ell q_\ell = 0, \qquad c = \{ \text{H}_2\text{O}, \text{CO}_2 \},
\label{eq:fluid-mass-balance}
\end{equation}

where $\phi$ is the porosity, $x_{c\ell}$ is the mass fraction of component $c$ in phase $\ell$, $\rho_\ell$ is the density of phase $\ell$, $s_\ell$ is the saturation of phase $\ell$, $\boldsymbol{v}_\ell$ is the Darcy velocity of phase $\ell$, and $q_\ell$ is the source flux for phase $\ell$.

Using the multiphase extension of Darcy's law, we write that the Darcy velocity is proportional to the gradient of the pressure:

\begin{equation}
\boldsymbol{v}_\ell = -\frac{k_{r\ell}}{\mu_\ell} \bar{\mathbf{k}} \cdot \nabla p_\ell, \qquad \ell = \{ a, g \},
\label{eq:darcy-eq}
\end{equation}

where $k_{r\ell}$ is the relative permeability of phase $\ell$, $\mu_\ell$ is the viscosity of phase $\ell$, $\bar{\mathbf{k}}$ is the permeability tensor, and $p_{\ell}$ is the phase pressure.

In this work, we assume that the permeability tensor is diagonal with an equal value for each entry.

The system is closed with the following constraints

\begin{align}
&s_g + s_a = 1,  \qquad &\text{(saturation constraint)} \\
&p_g - p_a = p_c(s_g), \qquad &\text{(capillary pressure constraint)} \\
&x_{\text{H}_2\text{O},\ell} + x_{\text{CO}_2,\ell} = 1, \qquad \ell \in \{ a, g \}, \qquad &\text{(component fraction constraints)}
\end{align}

as well as standard thermodynamics constraints on fugacities.

The partitioning of the mass of the CO$_2$ component between the gas phase and the aqueous phase is determined as a function of pressure, temperature, and salinity using the model of Duan and Sun [Duan and Sun, 2003].

The gas phase densities and viscosities are computed using the Span-Wagner [Span and Wagner, 1996] and Fenghour-Wakeham [Fenghour et al., 1998] correlations,
respectively, while the brine properties are obtained using the methodology of Phillips et al. [Phillips et al., 1981].

The relative permeabilities are computed with the Brooks-Corey model as $k_{rg}(s_g) = 0.95 ( s_g / s_{g,max} )^2$ and $k_{ra}(s_a) = ( (s_a - s_{a,min}) / (1-s_{a,min}) )^6$.

The capillary pressure $p_c$ is computed from $s_g$ using the Leverett J-function relationship.

The domain is initially saturated with brine, with an initial pressure of 10 MPa and an initial temperature of 143.76 degrees Celsius.

We use analytical (Carter-Tracy) aquifer boundary conditions.

A well injects pure supercritical CO$_2$ at a rate of 0.058 kg/s for 950 days, assuming a storage reservoir with unit meter thickness.

§.§ Finite-volume discretization on unstructured polygonal meshes

Considering a mesh with $n_{C}$ cells, the system of equations (<ref>)-(<ref>) is discretized with a cell-centered, fully implicit (backward-Euler) finite-volume scheme based on a two-point flux approximation (TPFA) and single-point upstream weighting.

The primary variables are chosen to be the gas (CO$_2$-rich) phase pressure $p = p_g$, the overall component densities $\rho_{\text{H}_2\text{O}}$, and $\rho_{\text{CO}_2}$, where an overall component density represents the mass of a given component per unit volume of mixture.

The primary variables
can be related to the variables of equation (<ref>) using the formulas given in [Voskov and Tchelepi, 2012].

At each time step, the nonlinear system of discretized equations is solved with Newton's method with damping to update all the primary variables in a fully coupled fashion.

All the simulations are performed with GEOS, an open-source multiphysics simulator targeting geological carbon storage and other subsurface energy systems [Bui et al., 2021, Camargo et al., 2022, Cusini et al., 2022].

Unstructured polygonal meshes are well suited to represent complex faults and to perform local spatial refinement around the wells.

In this work, we focus on the class of perpendicular bisector (PEBI) grids [Palagi and Aziz, 1994, Jenny et al., 2002, Mlacnik et al., 2006, Meng et al., 2018, Borio et al., 2021] generated with MRST (Matlab Reservoir Simulation Toolbox) [Lie, 2016] to mesh a 1 km x 1 km x 1 m domain containing an injector well and two straight impermeable faults in a conforming fashion (see an example mesh shown in Figs. <ref>(a)).

As PEBI meshes are orthogonal by construction, the flow simulations can be performed with the TPFA finite-volume scheme without compromising the accuracy of the solution.

§ RECURRENT GNN SURROGATE MODEL

In this section, we first describe the graph representation of unstructured mesh-based simulation data. Then, we introduce the key aspects of the deep-learning-based surrogate model for two-phase subsurface flow problems, including the learning task, model architecture, training procedure, and data preparation.

Construction of the graph representation from an unstructured mesh: (a) an example mesh with two impermeable faults and one injector; highlighted areas are zoomed-in sections of (b)  cells along fault line #2; (c) internal cells and their graph representations; and (d) input node features and edge features as well as node output (node label).

§.§ Input and output graph representation

To leverage the capabilities of graph-based machine learning, we represent the unstructured mesh-based input and output data at a given time as a directed graph with properties.

The graph is denoted by $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ where $\mathcal{V} $ and $ \mathcal{E} $ are node and edge sets, respectively.

As shown in Fig. <ref>(a), each mesh cell $i$ is represented by a graph node $i \in \mathcal{V}$.

Two adjacent cells $i$ and $j$ connected with non-zero transmissibility are represented by two graph edges, with edge $(i,j) \in \mathcal{E}$ pointing to node $j$ and edge $(j,i)$ pointing to node $i$.

The transmissibility $T_{ij}$ at each connection between cells $i$ and $j$ is computed as a function of grid geometry and rock permeability as explained in [Lie, 2016, Lie and Møyner, 2021], and is equal to zero when two cells are separated by a fault.

For a given set of mesh and permeability, the graph structure remains
fixed for all time steps.

The properties associated with node $i$ at time $t_n$ are termed as node feature, $\mathbf{x}_i^n : \mathcal{V} \rightarrow \mathbb{R}^{n_N}$.

The properties corresponding to edge $(i, j)$ are referred to as edge feature, $\mathbf{f}_{(i, j)}: \mathcal{E} \rightarrow \mathbb{R}^{n_E}$ and are independent of time.

The dimensions of the node and edge features, $n_N$ and $n_E$, are specified below.

As illustrated in Fig. <ref>(d), the node feature $\mathbf{x}_i^n$ consists of a combination of dynamic variables at time $t_n$, $y_i^{n} \in \mathbb{R}$, and static model features $\mathbf{m}_i \in \mathbb{R}^{n_M}$, where $n_M$ is the dimension of the static model features.

Specifically, $y_i^{n}$ denotes the state variable whose dynamics are learned by the neural network, such as phase saturation, $s_{g,i}$, or pore pressure, $p_{g,i}$.

Note that in this work, for each dynamic variable ($s_g$ or $p_g$), we train one individual prediction model.

The static model parameters $\mathbf{m}_i$ include the scalar permeability $k_i$ (one dimension), the cell volume $V_i$ (one dimension), the cell center coordinates $\boldsymbol{x}_i$ (two dimensions), and the cell type $\boldsymbol{n}_i$ (a one-hot vector of dimension 4), such that the dimension of the model parameters is $n_M = 1 + 1 + 2 + 4 = 8$ and the dimension of the node features is $n_N = n_M + 1 = 9$.

As shown in Figs. <ref> (b) to (c), the cell type is used to identify cells playing a key role in the simulation, such as the cells along the faults and the cells where source terms (well) or boundary conditions are imposed.

The fourth cell type includes the remaining cells (not along faults, and not where source and boundary conditions are imposed).

The edge feature $\mathbf{f}_{(i,j)}$ of edge $(i, j)$ is constructed to enrich the graph connectivity information with the (signed) distance between cell centers, $\boldsymbol{x}_j – \boldsymbol{x}_i \in \mathbb{R}^2$, (see Fig. <ref>(d)) and its absolute value, such that $n_E = 3$.

The schematic expressions of node and edge features used in this study are given in Fig. <ref>(d)). We explored four configurations of node and edge inputs; the corresponding variables are detailed in Table <ref>.

§.§ MeshGraphNet-Long Short-Term Memory (MGN-LSTM)

The proposed MGN-LSTM model is designed to learn the spatio-temporal evolution of the selected dynamic variable (pressure or CO$_2$ saturation) of the two-phase flow problem defined in Section <ref>.

Given the initial state $\mathbf{Y}^{0} = [y^0_1, \dots, y^0_{n_{C}}]^T$ and the static model features $\mathbf{M} = [(\mathbf{m}_1)^T, \dots, (\mathbf{m}_{n_{C}})^T]^T$ defined at the nodes of $\mathcal{G}$, and given the features at all edges $\mathbf{F}$, we compute the sequence of dynamic variables $(\hat{\mathbf{Y}}^{1}, \dots ,\hat{\mathbf{Y}}^{n_{T}})$ in an autoregressive way as follows:

\begin{align}
& \hat{\mathbf{Y}}^{0} = \mathbf{Y}^{0}, \\
&\hat{\mathbf{Y}}^{n+1} = f_{\text{MGN-LSTM}, \theta} ( \mathcal{G}, \hat{\mathbf{Y}}^{n}, \mathbf{M}, \mathbf{F}), \qquad n \in \{1, \dots, n_{T}\},
\label{eq-RMGN}
\end{align}

where  $f_{\text{MGN-LSTM}, \theta}$ is the surrogate model parameterized by the model weights $\theta$.

The number of cells in the mesh is denoted by $n_C$ and the number of temporal snapshots by $n_T$.

Once the models are trained, the inference is identical to the training process of Eq. (<ref>) and  produces a surrogate model with all the learnable parameters being fixed.

The training process will be discussed in detail in the next section.

Workflow schematic for the recurrent GNN model proposed in this article, referred to as MGN-LSTM. Given a domain in which supercritical CO$_2$ is injected in the center, the goal of MGN-LSTM is to learn a sequence of dynamic variables. The domain is discretized using an unstructured mesh that conforms to two fault lines and an injection well. A GNN framework based on MGN is used to predict physical quantities from time step $n$ to time step $n+1$ until the end of the simulation. A graph-based recurrent neural network framework based on Graph ConvLSTM  [Seo et al., 2018] is employed to better capture and represent temporal dynamics.

In MGN-LSTM, the temporal evolution of the dynamic variable is captured by learning directly on the graph using latent node and edge features derived from the physical node and edge features reviewed in Section <ref> as in the recent work of Pfaff et al. [Pfaff et al., 2020].

Specifically, our algorithm, sketched in Fig. <ref> and described in the next sections, combines the encoder-processor-decoder procedure of MGN [Pfaff et al., 2020] with a graph-based sequence model named Graph ConvLSTM [Seo et al., 2018] to perform the learning task.

Let us consider a one-step prediction, i.e., the computation of $\hat{\mathbf{Y}}^{n+1}$ using $\hat{\mathbf{Y}}^{n}$ in Eq. (<ref>).

First, the encoder-processor steps detailed in Section <ref> map the physical node/edge features to the latent node/edge features.

Then, as described in Section <ref>, the latent space variables are used as input to the Graph ConvLSTM algorithm, which aims at retaining spatial-temporal information encoded in the recurrent memories.

Finally, the output of Graph ConvLSTM is decoded and mapped to the physical space using the procedure of Section <ref>.

At that point, the one-step prediction of the dynamic variable is complete.

In the following subsections, we define the main components of the proposed architecture.

§.§.§ Encoder and processor

The encoder is the first step of the prediction in MGN-LSTM.

Using the graph representation of Section <ref>, we compute the initial latent feature vectors at time $t_n$,
$\mathbf{v}_{i}^{n, 0} \in \mathbb{R}^{n_H}$ and $\mathbf{e}_{i}^{n,0} \in \mathbb{R}^{n_H}$ from the physical feature vectors, $\mathbf{x}_{i}^{n}$ and $\mathbf{f}_{(i,j)}$.

The hyperparameter $n_H$ denotes the size of the latent vectors.

The computation of the latent vectors is done using the node and edge multilayer perceptrons (MLPs), denoted respectively by $\operatorname{mlp}_{v}^{0}$ and $\operatorname{mlp}_{e}^{0}$, as follows:

\begin{equation}
\begin{aligned}
\mathbf{v}_{i}^{n, 0} &=\operatorname{mlp}_{v}^{0}\left(\mathbf{x}_{i}^{n}\right), \\
\mathbf{e}_{(i,j)}^{n, 0} &=\operatorname{mlp}_{e}^{0}\left(\mathbf{f}_{(i,j)}\right).
\end{aligned}
\label{encoder}
\end{equation}

The graph equipped with the initial latent node and edge features computed by the encoder at time $t$ in Eq. (<ref>) is the input to the processor.

The processor consists of $m$ message-passing steps computed in sequence.

At step $\ell \in \{ 1, \dots, m \}$ in the sequence, each graph edge feature $\mathbf{e}_{(i,j)}^{n, \ell}$ is updated using its value at the previous message-passing step $\ell-1$ and the values of the adjacent node features at step $\ell-1$, as follows

\begin{equation}
\label{processor1}
\mathbf{e}_{(i,j)}^{n, \ell}=\operatorname{mlp}_{e}^{\ell}\left([\mathbf{e}_{(i,j)}^{n,\ell-1}, \mathbf{v}_{i}^{n, \ell-1}, \mathbf{v}_{j}^{n, \ell-1}]\right), \qquad
\ell \in \{ 1, \dots, m \}
\end{equation}

to obtain the updated value.

In Eq. (<ref>), the operator $\operatorname{[\cdot]}$ concatenates the given arguments on the feature dimension.

The mapping in each step in the message passing is computed using MLP with residual connection and Rectified Linear Unit (ReLu) as the non-linear activation function.

Then, each graph node $\mathbf{v}_{i}^{n, \ell}$ is updated using its value at the previous message-passing step, $\ell-1$, and the aggregation of its incident edge features at step $\ell$:

\begin{equation}
\label{processor2}
\quad \mathbf{v}_{i}^{n, \ell}=\operatorname{mlp}_{v}^{\ell} \left( \operatorname{[} \mathbf{v}_{i}^{n,\ell-1}, \sum_{j \in \text{adj}(i)} \mathbf{e}_{(i,j)}^{n,\ell} \operatorname{]} \right), \qquad \ell \in \{ 1, \dots, m \}
\end{equation}

where $\text{adj}(i)$ is the set of nodes connected to node $i$.

Using Eqs. (<ref>)-(<ref>), the processor computes an updated set of node features that are then used by the Graph ConvLSTM to produce recurrent memories.

The update of the edge-based messages, $\mathbf{e}_{(i,j)}^{n, \ell}$, is key to the accuracy of the MGN flow predictions as it propagates information between neighboring graph nodes (i.e., between neighboring mesh cells).

This design choice differentiates MGN from other classical GNN frameworks relying only on node features (see [Zhou et al., 2020]), such as GCN and GraphSAGE.

Moreover, leveraging edge information makes it possible to capture nontrivial topological information regarding connectivity and transmissibility, which play an important role in HF simulations and could be used to infuse more physics into the data-driven model (see Section <ref>).

§.§.§ Graph-based convolutional LSTM model

To limit the temporal error accumulation and improve prediction accuracy, we complement MGN with a variant of Convolutional Long Short-Term Memory (ConvLSTM) that operates on graph data named Graph Convolutional LSTM (GConvLSTM) [Seo et al., 2018].

The latter is obtained by replacing the convolutional operator in ConvLSTM with a graph operator.

Specifically, we follow the choice made in [Seo et al., 2018] and replace the Euclidean 2D convolution kernel with the Chebyshev spectral convolutional kernel [Defferrard et al., 2016], whose hyperparameters are given in Appendix <ref>.

Graph spectral filters are known to perform effectively on graph-based data with a small number of parameters thanks to their isotropic nature [Seo et al., 2018].

The goal of the Graph ConvLSTM step at time $t_n$ is to compute the cell state $\mathbf{C}^n$ and the hidden state $\mathbf{H}^n$.

This is done as a function of the latent representation of node features computed  by the processor of MGN at time $t_n$ (Section <ref>) and of the recurrent memories $\mathbf{C}^{n-1}$ and $\mathbf{H}^{n-1}$.

Using the terminology of ConvLSTM, the GConvLSTM architecture involves a set of memory cells, namely the cell state $\mathbf{C}^n$ and the hidden state $\mathbf{H}^n$.

GConvLSTM also relies on input gates $\mathbf{i}^n$, output gates $\mathbf{o}^n$, and forget gates $\mathbf{f}^n$ defined in [Seo et al., 2018].

These gates are based on a graph convolution operator and are used to control the flow of information into, out of, and within the memory cells.

By construction, the cell state $\mathbf{C}^n \in \mathbb{R}^{n_C \times n_H}$ and the hidden state $\mathbf{H}^n \in \mathbb{R}^{n_C \times n_H}$ exhibit temporal dynamics and can contain spatial structural information of the graph-based input $\boldsymbol{V}^{n,m} = [(\mathbf{v}^{n,m}_1)^T, \dots, (\mathbf{v}^{n,m}_{n_C})^T]^T$ at time $t_n$.

The functions of GConvLSTM are computed as follows:

\begin{equation}
\begin{aligned}
\mathbf{i}^n &=\sigma\left( \mathbf{W}_{x i} *_{\mathcal{G}} \boldsymbol{V}^{n,m}+\mathbf{W}_{h i} *_{\mathcal{G}} \mathbf{H}^{n-1} + \mathbf{b}_i\right), \\
\mathbf{f}^n &=\sigma\left(\mathbf{W}_{x f} *_{\mathcal{G}} \boldsymbol{V}^{n,m}+\mathbf{W}_{h f} *_{\mathcal{G}} \mathbf{H}^{n-1} + \mathbf{b}_f\right), \\
\mathbf{C}^n &= \mathbf{f}^n \odot \mathbf{C}^{n-1} + \mathbf{i}^n \odot \tanh \left( \mathbf{W}_{x c} *_{\mathcal{G}} \boldsymbol{V}^{n,m} +\mathbf{W}_{h c} *_{\mathcal{G}} \mathbf{H}^{n-1}+\mathbf{b}_c\right), \\
\mathbf{o}^n &=\sigma\left(\mathbf{W}_{x o} *_{\mathcal{G}} \boldsymbol{V}^{n,m} +\mathbf{W}_{h o} *_{\mathcal{G}} \mathbf{H}^{n-1} +\mathbf{W}_{c o} \mathbf{C}^n + \mathbf{b}_o\right), \\
\mathbf{H}^n &= \mathbf{o}^n \odot \tanh \left(\mathbf{C}^n\right),
\end{aligned}
\label{eq:lstm}
\end{equation}

where $*_{\mathcal{G}}$ denotes the graph convolution operator, $\odot$ denotes the Hadamard product, and $\sigma(\cdot)$ is the sigmoid activation function.

$\mathbf{W}$ and $\mathbf{b}$ are respectively the weights of the graph convolutional kernel and the bias term.

$\mathbf{i}^n$ determines how much of the new input is incorporated into the cell state. $\mathbf{f}^n$ controls the information to eliminate from the previous cell state. $\mathbf{o}^n$ determines how much of the cell state is output to the next time step.

In MGN-LSTM, the input vector $\boldsymbol{V}^{n,m} \in \mathbb{R}^{n_C \times n_H} $ is the latent representation of node features computed by the processor of MGN at time $t_n$ (Section <ref>).

In Eq. (<ref>), the vectors $\mathbf{C}^{n-1}$ and $\mathbf{H}^{n-1}$ are the recurrent memories obtained from the previous time step.

After evaluating Eq. (<ref>), the updated $\mathbf{H}^{n}$ is decoded into the next-step physical state as explained in Section <ref>.

More details regarding GConvLSTM gates and states can be found in [Seo et al., 2018].

§.§.§ Decoder

The decoder maps the updated hidden state $\mathbf{H}^{n}$ computed by GConvLSTM to the dynamic node-based properties in physical state using MLP as follows:

\begin{equation}
\label{decoder}
\hat{y}_{i}^{n+1}=\operatorname{mlp}_{v}^{m+1}\left(\mathbf{h}_{i}^{n}\right)
\end{equation}

where $\mathbf{h}_{i}^{n}$ contains the rows of the hidden state vector $\mathbf{H}^{n}$ corresponding to the updated latent vector of graph node $i$, and $m$ is the number of steps performed by the processor.

Detailed illustrations of the encoder, processor, LSTM cell, and decoder are given in Appendices <ref>, <ref>, and <ref>, respectively.

§.§ Loss function and Optimizer

We train MGN-LSTM on training time steps by minimizing the misfit between the true node label (HF simulation results) and the predicted node label.

We use the per-node root mean square error (RMSE) loss to quantify the data mismatch for each time step.

The loss function reads:

\begin{equation}
\mathcal{L}_{\text {MGN-LSTM}} =  \sum_{n=1}^{n_T} \sqrt{ \frac{1}{n_{B}} \sum_{i=1}^{n_{B}} \left\|y^{n}_i-\hat{y}^{n}_i\right\|_2^2}
\label{loss_fn}
\end{equation}

where $n_{B}$ denotes the number of nodes in a batch of training meshes, $y^{n}_i$ denotes the true output in the data set, $\hat{y}^{n}_i$ is the output predicted by MGN-LSTM, as formalized in Eq. (<ref>).

During training, the learning weights of $f_{\text{MGN-LSTM}, \theta}$ are updated based on the gradient to the loss function through back-propagation.

Unlike next-step models such as MeshGraphNet [Pfaff et al., 2020], MGN-LSTM can propagate gradients throughout the entire sequence, allowing the model to utilize information from previous steps to produce stable predictions with small error accumulation.

Furthermore, recurrent memories output by LSTM in MGN-LSTM retain information from previous inputs, which can be utilized to inform future predictions.

An adaptive moment estimation (ADAM) optimizer is used and the learning rate is gradually decreased from $10^{-4}$ to $10^{-6}$.

With the dataset and a fixed set hyperparameters at hand, each epoch takes 40s to 100s on an NVIDIA A100 GPU depending on the architecture dimension. The value of specific training hyperparameters is given in Appendix <ref>.

§ RESULTS

In this section, we first describe the dataset considered to train and test the MGN-LSTM model.

We illustrate that MGN-LSTM can accurately predict the CO$_2$ saturation plume and pore pressure evolution in the presence of impermeable faults.

At the end of this section, we also show how the saturation prediction accuracy can be improved by incorporating more physical properties (for instance, relative permeability) in the node and edge features of the graph.

The results of MGN-LSTM will be compared with those obtained with standard LSTM in Section <ref>.

§.§ Data description and MGN-LSTM training setup

We generate a total of 500 realizations of the synthetic geological models of size 1 km x 1 km x 1 m.

The domain shape as well as the position of the two impermeable faults are fixed across all realizations.

The coordinates of the endpoints of fault line 1 are (100m, 300m) and (400m, 600m), and the coordinates of the endpoints of fault line 2 are (400m, 500m) and (800m, 800m).

The synthetic models differ in their geological parameters (permeability), mesh configuration, and well location.

In each synthetic model, we first randomly generate the location of the injection well constrained within a prescribed 200 x 200 m box in the center of the domain to ensure that the injector is not placed too close to the boundary.

Then, a PEBI mesh conforming exactly to the specified faults and refined around the well is generated.

After that, we create a geomodel in which the permeability values are assigned to each cell according to a randomly generated Von Karman distribution [Carpentier and Roy-Chowdhury, 2009] using SGeMS [Remy et al., 2009]. Specifically, the mean and standard deviation of log-permeability are 3.912 ln(mD) and 0.5 ln(mD), respectively, which results in an average permeability of 50 mD in the reservoir. A constant porosity of 0.2 is assigned to all cells.

Fig. <ref> shows three sampled geomodels and meshes used for training.

High-fidelity numerical simulation is then performed for each model using the GEOS simulator.

MGN-LSTM is trained with 450 input meshes, in which each mesh has an 11-step rollout of simulation data, representing 550 days of CO$_2$ injection.

The rollout division is linear with 50 days per step.

The trained model is then tested on 50 unseen meshes for 11-step rollouts and 19-step rollouts with temporal extrapolation to 950 days.

The average number of graph nodes and graph edges in the dataset are 1885 and 4500 respectively.

We use two separate MGN-LSTM instances of the same architecture to predict the two different dynamical quantities, namely gas saturation and pore pressure.

The only difference between these two models is the dynamical quantity used to form node features and node labels during the training and inference stages.

In this study, we use `detrending' scaling [Tang et al., 2022] for all fields of node/edge features.

Specifically, the preprocessing can be expressed as:

\begin{equation}
\tilde{\mathbf{x}}^{n}_i = \frac{\mathbf{x}^{n}_i -  \text{mean}([\mathbf{x}_1^n, \ldots,\mathbf{x}_{n_{S}}^n ])}{\text{std}([\mathbf{x}_1^n, \ldots,\mathbf{x}_{n_{S}}^n ])}, \quad i=1,\ldots, n_{S}, \ \ n = 1, \ldots, n_T,
\label{eq:p-0}
\end{equation}

where $n_{S}$ is the number of training samples. With this approach, a field of node features at time $t_n$, $\mathbf{x}^{n}_i$, can be normalized by subtracting the mean field (over all $n_{S}$ samples) and dividing each element by its standard deviation at that time step.

Heterogeneous permeability realizations with two fixed impermeable faults and one injection well for three cases. The coordinates of well location in each sampled case is shown on the top, and the insets show the enlarged vicinity of each well.

§.§ Evaluation metrics

To quantify the prediction accuracy for gas saturation, we use the plume saturation error, $\delta^{s_{g}}$, introduced in [Wen et al., 2023] and defined as:

\begin{equation}
\begin{aligned}
& \delta^{s_{g}}=\frac{1}{\sum_{i,n} I^{n}_{i}} \sum_{n = 1}^{n_T} \sum_{i = 1}^{n_C} I^{n}_{i}\left|s^{n}_{g, i} - \hat{s}_{g, i}^n\right|, \\
& I_{i}^n =1 \quad \text { if } \quad\left(s_{g,i}^{n}>0.01\right) \cup\left(\left|\hat{s}^n_{g,i}\right|>0.01\right),
\end{aligned}
\label{gs_error}
\end{equation}

where $I^{n}_{i}=1$ indicates that a mesh cell has a non-zero gas saturation in either the ground truth or the prediction, $s^n_g$ denotes the gas saturation from HF simulations (ground truth), $\hat{s}^n_g$ is the predicted gas saturation, $n_T$ is the number of temporal snapshots, including training ranges (11 steps, 550 days) and extrapolated ranges (8 steps, 400 days), and $n_C$ is the number of cells in the mesh, which can vary between simulation models.

As discussed in [Wen et al., 2023], $\delta^{s_g}$ provides a strict metric to evaluate CO$_2$ gas saturation because it focuses on the accuracy within the separate phase plume.

We use the relative error $\delta^{p_g}$ defined below to evaluate the prediction accuracy for pore pressure:

\begin{equation}
\begin{aligned}
\delta^{p_g}=\frac{1}{n_{C} n_T} \sum_{n = 1}^{n_T} \sum_{i = 1}^{n_C} \frac{\left|p_{g,i}^{n} -\hat{p}_{g,i}^n\right|}{p_{ g,\textit{init}}},
\end{aligned}
\label{pp_error}
\end{equation}

where $p^n_g$ denotes the ground truth pore pressure given by the HF simulation, $\hat{p}^n_g$ is the predicted pore pressure, and $p_{ g,\textit{init} }$ is the initial reservoir pressure, which remains identical ($10$ MPa) for all test/training cases.

In the following sections, we use these metrics to illustrate the accuracy of MGN-LSTM in two steps.

First, in Section <ref>, we consider a representative mesh in the test set and demonstrate the ability of MGN-LSTM to capture the complex plume dynamics for time steps beyond the training period.

Then, in Section <ref>, we consider 10 representative test meshes to illustrate the reliability of MGN-LSTM predictions for unseen well locations, permeability fields, and meshes yielding very different CO$_2$ plume shapes.

§.§ Predicting complex spatio-temporal dynamics beyond the training period

In this section, we consider mesh 468 from the test set as an example.

The mesh and location of the well are highlighted in the insets of Fig. <ref>.

This test case has a CO$_2$ plume saturation error (Eq. (<ref>)) within the interquartile range of the test ensemble and can therefore be considered as representative of the accuracy of the surrogate model predictions (see Fig. <ref>).

This particular configuration is illustrated here because its well location is close to the faults and therefore produces an interesting saturation plume.

Figure <ref> compares the MGN-LSTM prediction of the CO$_2$ saturation fields with the HF simulation results at five snapshots ($t$ = 200, 400 days in the training period, $t$ = 600, 800, 950 days beyond the training period).

MGN-LSTM and HF simulations exhibit an excellent match.

Despite the presence of impermeable faults near the injector, we observe that (1) MGN-LSTM can accurately capture the complex temporal evolution of the CO$_2$ plume in the presence of impermeable faults, and (2) MGN-LSTM can extrapolate beyond the training horizon with a mild error accumulation at the saturation front.

The MGN-LSTM pore pressure predictions are presented in Fig <ref>.

They also exhibit a very good agreement with the HF simulation results during and after the training period.

Unlike for the saturation variable, the prediction errors are distributed over the entire domain due to the elliptic behavior of the pressure variable.

Temporal evolution of the CO$_2$ saturation plume. The first and second rows respectively show the CO$_2$ saturation fields from HF and MGN-LSTM for a representative test case (mesh 468) at five different times. The third row shows the saturation error between MGN-LSTM and HF. The first three columns represent predictions within the training period and the last two columns represent the model predictions for time steps beyond the training period.

Temporal evolution of the pore pressure. The first and second rows respectively show the pore pressure fields from HF and MGN-LSTM for a representative test case (mesh 468) at five different times. The third row shows the pressure error between MGN-LSTM and HF. The first two columns show predictions within the training period and the last three columns show the model predictions for time steps beyond the training period.

Accuracy of the MGN-LSTM prediction of (a) gas saturation and (b) pore pressure after rolling out 550 and 950 days for all test meshes. The meshes discussed in Section <ref> are highlighted as triangles and shown on top of the boxplot of the 950-day rollout results.

§.§ Generalizability to meshes, well locations, and permeabilities not seen during training

In this section, we consider 10 representative meshes from the test set (meshes 484, 494, 469, 490, 497 for the saturation and meshes 467, 451, 464, 450, 472 for pressure) to demonstrate that MGN-LSTM generalizes well for different meshes, boundary conditions (well locations) and permeability fields that are not included in the training set.

As shown in Fig. <ref>, the 10 test cases have a saturation error within the interquartile range of the test ensemble and are therefore representative of the surrogate accuracy.

CO$_2$ saturation predictions at 950 days for five distinct test meshes. From top to bottom, the rows show reservoir permeability, HF simulation of CO$_2$ saturation, MGN-LSTM prediction of CO$_2$ saturation, and MGN-LSTM prediction errors. The coordinates of well location in each test case is shown on the top of first row and the insets show the enlarged vicinity of each well.

The results after 950 days are presented in Figs. <ref>-<ref> for the CO$_2$ saturation and the pore pressure, respectively.

The figures show that the permeability field and well location relative to the faults vary drastically from one case to the other.

Due to the complex interplay between the saturation front and the faults, the differences in initial setup yield very different plume shapes.

Still, MGN-LSTM achieves an excellent agreement with the HF simulation results after 950 days for both saturation and pressure.

MGN-LSTM displays a remarkable generalizability considering the high dimensionality of the problem and the training data size.

The median saturation errors in the CO$_{2}$ plume prediction for the training set and the extrapolated ranges in the testing set are only $1.0\%$ and $1.3\%$, respectively.

This degree of accuracy in the prediction of CO$_2$ plume migration is sufficient for many practical applications, including the estimation of sweep efficiencies  [Wen et al., 2023].

Similarly, the pressure prediction exhibits an excellent accuracy.

The median pore pressure errors for the training set and the extrapolated ranges in the testing set stand at $0.16\%$ and $0.19\%$.

This confirms the reliability of MGN-LSTM predictions for complex configurations unseen during training.

Pore pressure predictions at 950 days for five distinct test meshes. From top to bottom, the rows show reservoir permeability, HF simulation of pore pressure, MGN-LSTM prediction of pore pressure, and MGN-LSTM prediction errors. The coordinates of well location in each test case is shown on the top of first row and the insets show the enlarged vicinity of each well.

§.§ Improving performance using augmented physics-based graph node/edge features

In this section, we describe a methodology to improve the prediction accuracy of MGN-LSTM.

In the previous sections, the edge feature of MGN-LSTM only accounted for distance-related information (see Section <ref>).

But, more physics insights can be infused into the algorithm by incorporating more information in the edge and node features.

Specifically, we explore three modifications of the gas saturation model discussed in earlier sections, focusing on the incorporation of additional features.

The first modification involves the inclusion of static transmissibility ($T_{ij}$) as an additional edge feature, the second variation introduces phase relative permeability ($k_{r,i}$) as an additional node feature, and the third variation combines both static transmissibility and phase relative permeability as extra edge and node features, respectively.

The transmissibility  remains fixed throughout rollout steps.

The relative permeability of gas phase $k_{r,i}$ is computed as a function of the predicted gas saturation, $s_{g,i}^{n-1}$.

The input and output variables for each case are summarized in Table <ref>.

Input and output of the surrogate model: $s_{g,i}^n, V_i, k_i, \boldsymbol{n}_i, \boldsymbol{x}_i, k_{r,i}(s_{g,i}^{n-1})$ denote gas saturation, cell volume, scalar permeability, cell type, cell center, cell relative permeability, respectively, for a given cell.

Note that $k_{r,i}$ takes the gas saturation at the previous prediction step, $s_{g,i}^{n-1}$, to update its value.

Here, $\boldsymbol{n}_i$ is a one-hot vector of size 4, encoding whether a cell is an internal cell, injector, cell along fault lines, and boundary cell.

$T_{i,j}$ denotes the static transmissibility between cells $i$ and $j$.

Case     Node
input     Edge
input     Node

Baseline     $s_{g,i}^n, V_i, k_i, \boldsymbol{n}_i, \boldsymbol{x}_i $     $\boldsymbol{x}_j-\boldsymbol{x}_i$,

$\left|\boldsymbol{x}_j-\boldsymbol{x}_i\right|$     $s_{g,i}^{n+1}$

Static transmissibility     $s_{g,i}^n, V_i, k_i, \boldsymbol{n}_i, \boldsymbol{x}_i $     $\boldsymbol{x}_i-\boldsymbol{x}_j$,

$T_{i,j}$     $s_{g,i}^{n+1}$

Relative permeability     $s_{g,i}^n, V_i, k_i, \boldsymbol{n}_i, \boldsymbol{x}_i, k_{r,i}(s_{g,i}^{n-1})$     $\boldsymbol{x}_i-\boldsymbol{x}_j$,

$\left|\boldsymbol{x}_i-\boldsymbol{x}_j\right|$      $s_{g,i}^{n+1}$

Static transmissibility

Relative permeability     $s_{g,i}^n, V_i, k_i, \boldsymbol{n}_i, \boldsymbol{x}_i, k_{r,i}(s_{g,i}^{n-1})$     $\boldsymbol{x}_i-\boldsymbol{x}_j$,

$T_{i,j}$       $s_{g,i}^{n+1}$

These three modifications are evaluated separately and incrementally, and are only applied to the prediction of CO$_2$ saturation.

We train each experiment with different random seeds.

We evaluate the effect of augmented features on model performance by measuring the 11-step rollout plume error, which only includes the training range, namely the first 550 days.

Figure <ref> confirms that adding more physical information in the node/edge features clearly improves CO$_2$ saturation prediction accuracy.

Using the transmissibility as an edge yields a mild error reduction in $\delta^{s_{g}}$, while the largest improvement is observed for the addition of relative permeability in the node feature with a near 10%-reduction in $\delta^{s_{g}}$. This experiment demonstrated that incorporating more physics into the architecture can improve prediction accuracy even in a purely data-driven framework.

Effects of injecting more physics information on the 11-step rollout (550 days of injection) performances of MGN-LSTM at predicting CO$_2$ saturation. The error bars correspond to min/median/max performance across 12 random seeds.
"Trans" denotes the case with transmissibility included in the edge features, "Rela. Perm." denotes the case with relative permeability included in the node features, and "Trans + Rela. Perm." combines the two approaches.

§ DISCUSSION

§.§ Comparison between MGN-LSTM and standard MGN

The standard MGN approach [Pfaff et al., 2020] suffers from temporal error accumulation and often requires mitigation strategies to maintain a stable long-term prediction.

In the present study, we rely on a noise injection technique [Han et al., 2022, Pfaff et al., 2020] to enable a stable rollout of standard MGN.

We compare standard MGN and MGN-LSTM with the same encoder-processor-decoder architecture.

In MGN-LSTM, two separate standard MGN models are trained for predicting gas saturation and pore pressure.

The accuracy improvements discussed in Section <ref> are not used here.

The noise injection strategy and training procedure for MGN are detailed in Appendix <ref>, as well as the model parameters used in the study.

The goal of this ablation-nature comparison is to demonstrate the ability of LSTM to constrain temporal errors.

Prediction accuracy for (a) gas saturation and (b) pore pressure after rolling out 550 and 950 days for mesh 468 (see Section <ref>).

Next, we briefly illustrate the difference in prediction accuracy between MGN-LSTM and standard MGN for the meshes considered in Sections <ref> and <ref>.

Considering mesh 468 of Section <ref> first, we compare the predictions of standard MGN and MGN-LSTM after the training period in Fig. <ref>.

It is clear that for both saturation and pressure, the standard MGN predictions start deviating significantly from the HF data after 800 days.

The standard MGN error accumulates at the front for CO$_2$ saturation and is more diffused for pressure.

In both cases however, the MGN-LSTM error remains constrained thanks to the addition of the LSTM cell.

Considering now the 10 meshes of Section <ref>, Fig. <ref> confirms the poor prediction accuracy of standard MGN after 950 days of injection for both CO$_2$ saturation and pressure compared to MGN-LSTM.

Prediction accuracy for (a) gas saturation and (b) pore pressure after rolling out 550 and 950 days for all test meshes of Section <ref>.

We summarize the respective accuracy of MGN-LSTM and standard MGN by comparing their ensemble results.

Figure <ref> shows the boxplots of error metrics $\delta^{S_g}$ and $\delta^{P_g}$ evaluated at (1)  550 days (the end of training) and (2) 950 days (400 days of extrapolation).

We observe that in both cases, MGN-LSTM outperforms standard MGN in terms of accuracy for the prediction of CO$_2$ saturation as well as pressure.

This is particularly the case after the end of the training period (i.e., between 550 days and 950 days), as the standard MGN saturation prediction accuracy deteriorates significantly due to error accumulation.

Boxplots of the test ensembles for plume saturation and pressure errors of MGN-LSTM and standard MGN. Panels (a) and (c) compare the CO$_2$ saturation predictions of MGN-LSTM and standard MGN after rolling out 11 steps (550 days) and 19 steps (950 days), respectively, on the test ensembles. Panel (b) and (d) compare the pore pressure predictions of MGN-LSTM and standard MGN after rolling out 11 steps (550 days) and 19 steps (950 days), respectively, on the test ensembles. MGN-LSTM is in orange and standard MGN is in blue.

§.§ Computational efficiency

To analyze the computational efficiency, we compare the inference times of MGN-LSTM, standard MGN, and the high-fidelity simulator, GEOS, as demonstrated in Table <ref>. Utilizing the previously discussed dataset with 1,885 cells per mesh on average, MGN-LSTM requires an average of 0.31 seconds for a 19-step rollout on an NVIDIA Tesla A100 GPU to process a single batch. This contrasts with standard MGN's shorter average time requirement of 0.07 seconds on the same hardware. The inference time for MGN-LSTM is around four times that of standard MGN, primarily due to the auto-regressive prediction of dynamic predictions in LSTM. The reduction of this overhead is an area of focus for future work.

Comparing these surrogate models with GEOS demonstrates a significant performance gain. Specifically, on the same dataset, MGN-LSTM exhibits a nearly 160-fold reduction in execution time in comparison to GEOS, which operates on a CPU Intel(R) Xeon(R) E5-2680 v4 2.10GHz. We expect that this performance gain will be even more prominent with a larger mesh size.

Average inference times for MGN-LSTM and standard MGN with 11-step rollout (550 days) and 19-step rollout (950 days). The inference times are compared with the GEOS run time after 550 days and 950 days.

MGN-LSTM avg.
inference time (s)a     Standard MGN avg.
inference time (s)b     GEOS
run time (s)c

11-step rollout (550 days)      0.18     0.04     22.12

19-step rollout (950 days)      0.31     0.07     49.02

a On an NVIDIA Tesla A100 GPU, single-batch inference run

b On an NVIDIA Tesla A100 GPU, single-batch inference run

c On an Intel Xeon E5–2695 v4, single-core serial run

§ CONCLUDING REMARKS

This paper aims at demonstrating the applicability of deep-learning-based surrogate models for subsurface flow simulations with complex geological structures.

We present a graph-based neural surrogate model that can operate on unstructured meshes and naturally handle geological fault structures.

Our model combines a graph-based Long-Short-Term-Memory (LSTM) cell with a one-step graph neural network model, MeshGraphNet (MGN), to control temporal error accumulation.

The model is trained on 450 high-fidelity, unstructured simulation results.

The mesh configuration, well location, and reservoir permeability field vary from one simulation to the other.

The accuracy and performance of MGN-LSTM were analyzed on a set of 50 test cases.

Our results demonstrate that MGN-LSTM can accurately predict the temporal evolution of gas saturation and pore pressure in a CO$_2$ injection case with two impermeable faults.

The model exhibits excellent generalizability to mesh configurations, well locations, and permeability fields not included in the training set.

Furthermore, our comparison study shows that MGN-LSTM is more accurate than standard MGN and predicts dynamic fields with a smaller temporal error accumulation thanks to the LSTM cell.

These conclusions have key implications in the field of CO$_2$ storage since accurately capturing uncertainties related to geological structures is critical.

By enabling the construction of surrogates on unstructured meshes, MGN-LSTM has the potential to accelerate the quantification of uncertainties in CO$_2$ injection cases with complex geological structures such as faults and fractures.

Future work in this area should address several topics.

We observed that using full-graph embeddings to retain recurrent memories in LSTM could require a significant amount of memory, limiting the size of the model that could be trained.

Future work could explore alternative architecture designs to avoid using full-graph embeddings.

Furthermore, the model training is currently limited to 2D cases, and scaling up to larger domains remains an open challenge.

Finally, combining data assimilation techniques with the proposed surrogate model is also of interest.

§ CODE AVAILABILITY

The MGN-LSTM model architecture code and training dataset used in this study can be accessed at <https://github.com/IsaacJu-debug/gnn_flow> upon the publication of this manuscript.

§ ACKNOWLEDGEMENTS

We are grateful for the funding provided by TotalEnergies through the FC-MAELSTROM project. We also acknowledge the Livermore Computing Center and TotalEnergies for providing the computational resources used in this work. We thank Nicola Castelletto for his insights, fruitful discussions, and help with PEBI mesh generation and support in GEOS. Portions of this work were performed during a summer internship of the first author with TotalEnergies in 2022.

§ APPENDIX

§.§ Details of MGN-LSTM

§.§.§ Encoder

This subsection gives the detailed structure for the node level MLP ($\operatorname{mlp}_{v}$) and edge level MLP ($\operatorname{mlp}_{e}$), as shown in Table <ref>.

Details of node and edge encoders

Node MLP $\left(\operatorname{mlp}_{v}\right)$     Edge MLP $\left(\operatorname{mlp}_{e}\right)$

Input: $x$     Input: $x$

$x=\operatorname{Linear}(8,100)(x)$     $x=\operatorname{Linear}(3,100)(x)$

$x=\operatorname{RELU}(x)$     $x=\operatorname{RELU}(x)$

$x=\operatorname{Linear}(100,100)(x)$     $x=\operatorname{Linear}(100,100)(x)$

$x=\operatorname{LayerNorm}(100)(x)$     $x=\operatorname{LayerNorm}(100)(x)$

§.§.§ Processor

As discussed in Section <ref>, the processor is a stack of $m$ graph neural network (GNN) layers, each one performing a round of message passing. The hyperparameter values used for the processor are detailed in Table <ref>.

Except for the input dimension, the structure of $\operatorname{mlp}_{v}$ and $\operatorname{mlp}_{e}$ employed in the processor are the same with the ones used in the encoder (Table <ref>).

Hyperparameters used for processor

Parameter name     Value

Number of GNN layer for the processor     10

Latent size for the processor     100

Activation     Relu

Type of normalization     Layer normalization

Input feature size of Node MLP for one GNN layer     200

Input feature size of Edge MLP for one GNN layer     300

§.§.§ Graph Convolutional LSTM cell

The GConvLSTM cell employs identical input and output channel dimensions, both set to the processor's latent size of 100. The  filter size of Chebyshev spectral convolutional kernel is set to 8, with symmetric normalization being applied to the graph Laplacian [Seo et al., 2018].

§.§.§ Training hyperparameters

All hyperparameters employed in the training process are listed in Table <ref>.

Hyperparameters used for training

Parameter name     value

Loss function     RMSE

Number of epochs     500

Batch size     10

Learning rate     0.001

Weight decay     $5 \mathrm{e}-4$

Optimizer     Adam

Optimizer scheduler     cos

§.§ Details of MeshGraphNet (MGN) baseline

Unlike MGN-LSTM, where entire sequences of training steps are trained, MGN training involves only minimizing data mismatch between predicted node label and true node label for a single time step.

Moreover, MGN utilizes the same encoder-processor-decoder components and training hyperparameters as those used in MGN-LSTM.

§.§.§ Noise injection in MGN

Following the noise injection technique given in [Pfaff et al., 2020], we apply a zero-mean Gaussian noise to the field of dynamical quantities, such as gas saturation, during training.

The standard deviation, $s_{n}$, of the noise is a hyperparameter requiring tuning to achieve a stable rollout.

In this work, we have found $s_{n} = 0.05\cdot s $ to be a good choice, where $s$ is the standard deviation of dynamical quantities in the dataset.

We remark that the proposed MGN-LSTM model does not need perturbing corresponding dynamical quantities in the training dataset during the training.

[Agency, 2020]
authorI. E. Agency, titleEnergy Technology Perspectives
2020-Special Report on Carbon Capture Utilisation and Storage: CCUS in Clean
Energy Transitions, publisherOECD Publishing,

[Pacala and Socolow, 2004]
authorS. Pacala, authorR. Socolow,
titleStabilization wedges: solving the climate problem for the
next 50 years with current technologies, journalscience
volume305 (number5686) (year2004)

[Mallison et al., 2014]
authorB. Mallison, authorC. Sword,
authorT. Viard, authorW. Milliken,
authorA. Cheng, titleUnstructured cut-cell grids for
modeling complex reservoirs, journalSpe Journal
volume19 (number02) (year2014)

[Rinaldi and Rutqvist, 2013]
authorA. P. Rinaldi, authorJ. Rutqvist,
titleModeling of deep fracture zone opening and transient ground
surface uplift at KB-502 CO2 injection well, In Salah, Algeria,
journalInternational Journal of Greenhouse Gas Control
volume12 (year2013) pages155–167.

[Morris et al., 2011]
authorJ. P. Morris, authorY. Hao,
authorW. Foxall, authorW. McNab, titleA
study of injection-induced mechanical deformation at the In Salah CO2 storage
project, journalInternational Journal of Greenhouse Gas Control
volume5 (number2) (year2011)

[Castelletto et al., 2013]
authorN. Castelletto, authorP. Teatini,
authorG. Gambolati, authorD. Bossie-Codreanu,
authorO. Vincké, authorJ.-M. Daniel,
authorA. Battistelli, authorM. Marcolini,
authorF. Donda, authorV. Volpi,
titleMultiphysics modeling of CO2 sequestration in a faulted
saline formation in Italy, journalAdvances in Water Resources
volume62 (year2013) pages570–587.

[Cusini et al., 2022]
authorM. Cusini, authorA. Franceschini,
authorL. Gazzola, authorT. Gazzola,
authorJ. Huang, authorF. Hamon,
authorR. Settgast, authorN. Castelletto,
authorJ. White, titleField-scale simulation of
geologic carbon sequestration in faulted and fractured natural formations,
typeTech. Rep., institutionLawrence Livermore National
Lab. (LLNL), Livermore, CA (United States), year2022.

[Jenny et al., 2002]
authorP. Jenny, authorC. Wolfsteiner,
authorS. H. Lee, authorL. J. Durlofsky,
titleModeling flow in geometrically complex reservoirs using
hexahedral multiblock grids, journalSPE Journal
volume7 (number02) (year2002)

[Mlacnik et al., 2006]
authorM. J. Mlacnik, authorL. J. Durlofsky,
authorZ. E. Heinemann, titleSequentially adapted
flow-based PEBI grids for reservoir simulation, journalSPE
Journal volume11 (number03) (year2006)

[Meng et al., 2018]
authorX. Meng, authorZ. Duan,
authorQ. Yang, authorX. Liang, titleLocal
PEBI grid generation method for reverse faults, journalComputers
& Geosciences volume110 (year2018)

[Zhu and Zabaras, 2018]
authorY. Zhu, authorN. Zabaras,
titleBayesian deep convolutional encoder–decoder networks for
surrogate modeling and uncertainty quantification, journalJournal
of Computational Physics volume366 (year2018)

[Kitanidis, 2015]
authorP. K. Kitanidis, titlePersistent questions of
heterogeneity, uncertainty, and scale in subsurface flow and transport,
journalWater Resources Research
volume51 (number8) (year2015)

[Strandli et al., 2014]
authorC. W. Strandli, authorE. Mehnert,
authorS. M. Benson, titleCO2 plume tracking and
history matching using multilevel pressure monitoring at the Illinois
Basin–Decatur Project, journalEnergy Procedia
volume63 (year2014) pages4473–4484.

[Bischof and Kraus, 2021]
authorR. Bischof, authorM. Kraus,
titleMulti-objective loss balancing for physics-informed deep
learning, journalarXiv preprint arXiv:2110.09813 .

[Tang et al., 2020]
authorM. Tang, authorY. Liu, authorL. J.
Durlofsky, titleA deep-learning-based surrogate model for data
assimilation in dynamic subsurface flow problems, journalJournal
of Computational Physics volume413 (year2020)

[Tang et al., 2021]
authorH. Tang, authorP. Fu, authorC. S.
Sherman, authorJ. Zhang, authorX. Ju,
authorF. Hamon, authorN. A. Azzolina,
authorM. Burton-Kelly, authorJ. P. Morris,
titleA deep learning-accelerated data assimilation and forecasting
workflow for commercial-scale geologic carbon storage,
journalInternational Journal of Greenhouse Gas Control
volume112 (year2021) pages103488.

[Mo et al., 2019]
authorS. Mo, authorY. Zhu,
authorN. Zabaras, authorX. Shi,
authorJ. Wu, titleDeep convolutional encoder-decoder
networks for uncertainty quantification of dynamic multiphase flow in
heterogeneous media, journalWater Resources Research
volume55 (number1) (year2019)

[Tang et al., 2022]
authorM. Tang, authorX. Ju, authorL. J.
Durlofsky, titleDeep-learning-based coupled flow-geomechanics
surrogate model for CO2 sequestration, journalInternational
Journal of Greenhouse Gas Control volume118
(year2022) pages103692.

[Wen et al., 2021]
authorG. Wen, authorM. Tang, authorS. M.
Benson, titleTowards a predictor for CO$_2$ plume migration using
deep neural networks, journalInternational Journal of Greenhouse
Gas Control volume105 (year2021)

[Wen et al., 2022]
authorG. Wen, authorZ. Li,
authorK. Azizzadenesheli, authorA. Anandkumar,
authorS. M. Benson, titleU-FNO—An enhanced Fourier
neural operator-based deep-learning model for multiphase flow,
journalAdvances in Water Resources volume163
(year2022) pages104180.

[Wen et al., 2023]
authorG. Wen, authorZ. Li, authorQ. Long,
authorK. Azizzadenesheli, authorA. Anandkumar,
authorS. M. Benson, titleReal-time high-resolution CO
2 geological storage prediction using nested Fourier neural operators,
journalEnergy & Environmental Science
volume16 (number4) (year2023)

[Gao et al., 2022]
authorH. Gao, authorM. J. Zahr, authorJ.-X.
Wang, titlePhysics-informed graph neural Galerkin networks: A
unified framework for solving PDE-governed forward and inverse problems,
journalComputer Methods in Applied Mechanics and Engineering
volume390 (year2022) pages114502.

[Shukla et al., 2022]
authorK. Shukla, authorM. Xu,
authorN. Trask, authorG. E. Karniadakis,
titleScalable algorithms for physics-informed neural and graph
networks, journalData-Centric Engineering volume3
(year2022) pagese24.

[Sanchez-Gonzalez et al., 2020]
authorA. Sanchez-Gonzalez, authorJ. Godwin,
authorT. Pfaff, authorR. Ying,
authorJ. Leskovec, authorP. Battaglia,
titleLearning to simulate complex physics with graph networks,
in: booktitleInternational Conference on Machine Learning,
organizationPMLR, pages8459–8468,

[Pfaff et al., 2020]
authorT. Pfaff, authorM. Fortunato,
authorA. Sanchez-Gonzalez, authorP. W. Battaglia,
titleLearning mesh-based simulation with graph networks,
journalarXiv preprint arXiv:2010.03409 .

[Lino et al., 2021]
authorM. Lino, authorC. Cantwell,
authorA. A. Bharath, authorS. Fotiadis,
titleSimulating continuum mechanics with multi-scale graph neural
networks, journalarXiv preprint arXiv:2106.04900 .

[Shi et al., 2022]
authorN. Shi, authorJ. Xu, authorS. W.
Wurster, authorH. Guo, authorJ. Woodring,
authorL. P. Van Roekel, authorH.-W. Shen,
titleGNN-Surrogate: A Hierarchical and Adaptive Graph Neural
Network for Parameter Space Exploration of Unstructured-Mesh Ocean
Simulations, journalIEEE Transactions on Visualization and
Computer Graphics volume28 (number6)
(year2022) pages2301–2313.

[Wu et al., 2022]
authorT. Wu, authorQ. Wang, authorY. Zhang,
authorR. Ying, authorK. Cao,
authorR. Sosic, authorR. Jalali,
authorH. Hamam, authorM. Maucec,
authorJ. Leskovec, titleLearning Large-scale
Subsurface Simulations with a Hybrid Graph Network Simulator, in:
booktitleProceedings of the 28th ACM SIGKDD Conference on
Knowledge Discovery and Data Mining, pages4184–4194,

[Ronneberger et al., 2015]
authorO. Ronneberger, authorP. Fischer,
authorT. Brox, titleU-net: Convolutional networks for
biomedical image segmentation, in: booktitleMedical Image
Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International
Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18,
organizationSpringer, pages234–241,

[Han et al., 2022]
authorX. Han, authorH. Gao, authorT. Pffaf,
authorJ.-X. Wang, authorL.-P. Liu,
titlePredicting Physics in Mesh-reduced Space with Temporal
Attention, journalarXiv preprint arXiv:2201.09113 .

[Seo et al., 2018]
authorY. Seo, authorM. Defferrard,
authorP. Vandergheynst, authorX. Bresson,
titleStructured sequence modeling with graph convolutional
recurrent networks, in: booktitleInternational Conference on
Neural Information Processing, organizationSpringer,
pages362–373, year2018.

[Duan and Sun, 2003]
authorZ. Duan, authorR. Sun, titleAn
improved model calculating CO2 solubility in pure water and aqueous
NaCl solutions from 273 to 533 K and from 0 to 2000 bar,
journalChemical Geology
volume193 (number3-4) (year2003)

[Span and Wagner, 1996]
authorR. Span, authorW. Wagner, titleA new
equation of state for carbon dioxide covering the fluid region from the
triple-point temperature to 1100 K at pressures up to 800 MPa,
journalJournal of Physical and Chemical Reference Data
volume25 (number6) (year1996)

[Fenghour et al., 1998]
authorA. Fenghour, authorW. A. Wakeham,
authorV. Vesovic, titleThe viscosity of carbon
dioxide, journalJournal of Physical and Chemical Reference Data
volume27 (number1) (year1998)

[Phillips et al., 1981]
authorS. L. Phillips, authorA. Igbene,
authorJ. A. Fair, authorH. Ozbek,
authorM. Tavana, titleTechnical databook for
geothermal energy utilization, typeTech. Rep.,
institutionLawrence Berkeley National Laboratory, CA, USA,

[Voskov and Tchelepi, 2012]
authorD. V. Voskov, authorH. A. Tchelepi,
titleComparison of nonlinear formulations for two-phase
multi-component EoS based simulation, journalJournal of Petroleum
Science and Engineering volume82 (year2012)

[Bui et al., 2021]
authorQ. M. Bui, authorF. P. Hamon,
authorN. Castelletto, authorD. Osei-Kuffuor,
authorR. R. Settgast, authorJ. A. White,
titleMultigrid reduction preconditioning framework for coupled
processes in porous and fractured media, journalComputer Methods
in Applied Mechanics and Engineering volume387
(year2021) pages114111.

[Camargo et al., 2022]
authorJ. Camargo, authorF. Hamon,
authorA. Mazuyer, authorT. Meckel,
authorN. Castelletto, authorJ. White,
titleDeformation Monitoring Feasibility for Offshore Carbon
Storage in the Gulf-of-Mexico, in: booktitleProceedings of
the 16th Greenhouse Gas Control Technologies Conference (GHGT-16), October
2022, pages1–11, year2022.

[Palagi and Aziz, 1994]
authorC. L. Palagi, authorK. Aziz, titleUse
of Voronoi grid in reservoir simulation, journalSPE Advanced
Technology Series volume2 (number02)
(year1994) pages69–77.

[Borio et al., 2021]
authorA. Borio, authorF. P. Hamon,
authorN. Castelletto, authorJ. A. White,
authorR. R. Settgast, titleHybrid mimetic
finite-difference and virtual element formulation for coupled poromechanics,
journalComputer Methods in Applied Mechanics and Engineering
volume383 (year2021) pages113917.

[Lie, 2016]
authorK.-A. Lie, titleAn introduction to reservoir
simulation using MATLAB: user guide for the Matlab Reservoir Simulation
Toolbox (MRST), journalSintef Ict volume118.

[Lie and Møyner, 2021]
authorK.-A. Lie, authorO. Møyner,
titleAdvanced modelling with the MATLAB reservoir simulation
toolbox, publisherCambridge University Press,

[Zhou et al., 2020]
authorJ. Zhou, authorG. Cui, authorS. Hu,
authorZ. Zhang, authorC. Yang,
authorZ. Liu, authorL. Wang, authorC. Li,
authorM. Sun, titleGraph neural networks: A review of
methods and applications, journalAI Open volume1
(year2020) pages57–81, ISSN

[Defferrard et al., 2016]
authorM. Defferrard, authorX. Bresson,
authorP. Vandergheynst, titleConvolutional neural
networks on graphs with fast localized spectral filtering,
journalAdvances in Neural Information Processing Systems

[Carpentier and Roy-Chowdhury, 2009]
authorS. Carpentier, authorK. Roy-Chowdhury,
titleConservation of lateral stochastic structure of a medium in
its simulated seismic response, journalJournal of Geophysical
Research: Solid Earth volume114 (numberB10).

[Remy et al., 2009]
authorN. Remy, authorA. Boucher,
authorJ. Wu, titleApplied geostatistics with SGeMS: A
user's guide, publisherCambridge University Press,