# Hybrid Newton method for the acceleration of well event handling in the simulation of CO2 storage using supervised learning

Antoine Lechevallier, Sylvain Desroziers, Thibault Faney, Eric Flauraud, Frédéric Nataf

► To cite this version:

Antoine Lechevallier, Sylvain Desroziers, Thibault Faney, Eric Flauraud, Frédéric Nataf. Hybrid Newton method for the acceleration of well event handling in the simulation of CO2 storage using supervised learning. Computers & Geosciences, 2025, 197, pp.105872. 10.1016/j.cageo.2025.105872. hal-04988094v2

HAL Id: hal-04988094

https://hal.science/hal-04988094v2

Submitted on 2 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
1

# Highlights

3

## Hybrid Newton method for the acceleration of well event handling in the simulation of CO₂ storage using supervised learning

4

Antoine Lechevallier, Sylvain Desroziers, Thibault Faney, Eric Flauraud, Frédéric Nataf

5

6

* **Innovative Approach:** We develop a machine-learning based Newton preconditioning method, showcasing its potential through a proof of concept.

7

8

* **Hybrid Initialization:** Our approach integrates Newton's method with a hybrid initialization strategy, enabling rapid convergence within the domain of quadratic convergence while preserving numerical guarantees.

9

10

* **Application in CO₂ Storage Simulation:** Our method is applied to accelerate the numerical simulation of CO₂ storage, using a open-source Python framework developed during this work.

11

12

* **Initialization Strategy:** We initialize the pressure using a linear approximation and the saturation using a Fourier Neural Operator neural network trained in a supervised way.

13

14

* **Performance Boost:** Our method significantly reduces the number of Newton iterations across various well events in two test cases, showcasing its efficiency.
# Hybrid Newton method for the acceleration of well event handling in the simulation of CO₂ storage using supervised learning

Antoine Lechevallier<sup>a,b</sup>, Sylvain Desroziers<sup>c</sup>, Thibault Faney<sup>a</sup>, Eric Flauraud<sup>a</sup> and Frédéric Nataf<sup>d</sup>

<sup>a</sup>IFP Energies nouvelles

<sup>b</sup>Sorbonne-Université

<sup>c</sup>Michelin

<sup>d</sup>Laboratoire J.L. Lions

## ARTICLE INFO

Keywords:

CO₂ storage

Reservoir simulation

Newton method

Non-linear preconditioning

Machine learning

Fourier Neural Operator

## ABSTRACT

Geological storage of CO₂ is an essential instrument for efficient Carbon Capture and Storage policies. Numerical simulations provide the solution to the multi-phase flow equations that model the behavior of the CO₂ injection site. However, numerical simulations of fluid flow in porous media are computationally demanding: it can take up to several hours on a HPC cluster in order to simulate one injection scenario for a large CO₂ reservoir if we want to accurately model the complex physical processes involved. This becomes a limiting issue when performing a large number of simulations, e.g. in the process of 'history matching'. Well events, such as opening and closure, cause important numerical difficulties due to their instant impact on the pressure and saturation unknowns. This often forces a drastic reduction of the time step size to be able to solve the non-linear system of equations resulting from the discretization of the continuous mathematical model. However, these specific well events in a simulation have a relatively similar impact across space and time. We propose a proof of concept methodology to alleviate the impact of well events during the numerical simulation of CO₂ storage in the underground by using a machine-learning based non-linear preconditioning. We complement the standard fully implicit solver by predicting an initialization of Newton's method directly in the domain of quadratic convergence using supervised learning. More specifically, we replace the initialization in pressure by a linear approximation obtained through an implicit solver and we use a Fourier Neural Operator (FNO) to predict the saturation initialization. Furthermore, we present an open-source Python framework for conducting reservoir simulations and integrating machine-learning models. We apply our methodology to two test cases derived from a realistic CO₂ storage in saline aquifer benchmark. We reduce the required number of Newton iterations to handle a well opening by 53% for the first test case, i.e required number of linear system to solve and by 38% for the second test case.

## CRediT authorship contribution statement

Antoine Lechevallier: PhD student: Redaction and code developement. Sylvain Desroziers: Supervising. Thibault Faney: Supervising. Eric Flauraud: Supervising. Frédéric Nataf: Supervising.

## 1. Introduction

The reduction of CO₂ emissions into the atmosphere is mandatory to achieve ecological transition. CO₂ geological storage is an essential instrument for efficient Carbon Capture and Storage (CCS) policies. Numerical simulations provide the solution to the multi-phase flow equations that model the behavior of the CO₂ injection site. They are an important tool to decide whether or not to operate a potential carbon storage site and to monitor these operations (long term storage of injected CO₂, potential gas leakage, optimal positioning of CO₂ injection wells, etc.).

ORCID(s): 0009-0004-4605-4140 (A. Lechevallier)
However, numerical simulations of fluid flow in porous media are computationally demanding: it can take up to
several hours on a HPC cluster in order to simulate one injection scenario for a large CO₂ reservoir if we want to
accurately model the complex physical processes involved. This becomes a limiting issue when performing a large
number of simulations.

Industrial solvers often rely on fully implicit methods (FIM) (Dawson et al., 1997; Palomino Monteagudo and Firoozabadi, 2007; Liu et al., 2013; Eymard et al., 2000; LeVeque, 2002), which are unconditionally stable, enabling the solution of problems across various time scales. However, this stability comes at the expense of solving a large system of equations with mixed elliptic/hyperbolic characteristics. The resulting nonlinear system of equations is typically tackled using Newton's method which is not unconditionally convergent (Cai et al., 2002; Skogestad et al., 2013). Indeed, to attain a solution within a defined tolerance for a desired time-step, a series of Newton iterations is executed until convergence is achieved for the nonlinear algebraic equations. Nevertheless, when confronted with poor initial guesses and/or large time-steps, the standard Newton method for nonlinear solvers may fail to converge. In such instances, a common recourse is to retry with a smaller time step, incurring significant computational time.

Over the past few decades, efforts have been directed towards enhancing the convergence of Newton's method for two-phase flow problems. This has been pursued through adjustments to the time-step size or the utilization of nonlinear preconditioners (Yang et al., 2016). These preconditioners serve to enhance the initial guess, employing either explicit (IMPES) (Chen et al., 2004, 2006) or implicit sequential methods (IMPIMS) Quandalle and Savary (1989); Ouaki (2013), or employing domain decomposition preconditioning methods such as the Additive Schwarz Preconditioned Inexact Newton method (ASPIN) (Cai and Keyes, 2001; Hwang and Cai, 2005; Liu and Keyes, 2015). Another preconditioning strategy at the Newton iteration level involves damping, where a damping factor is applied cell-wise to ensure that absolute variations in saturations remain below a certain threshold. This threshold may be heuristic (Schlumberger, 2013) or physics-based, such as Trust regions (Pour et al., 2023; Klemetsdal et al., 2019), which determine safe saturation region updates using inflection points.

Well events, such as opening and closure, cause important numerical difficulties due to their large impact on the pressure and saturation unknowns. This often forces a drastic reduction of the time step size to be able to solve the non-linear system of equations resulting from the discretization of the continuous mathematical model.

However, these specific well events in a simulation are relatively similar across space and time: the degree of similarity between two well events depends on a few parameters such as the injection condition, the state of the reservoir at the time of the event, the boundary conditions or the porous media parameters (permeability and porosity) around each well. Therefore, the potential of machine learning algorithms to alleviate numerical difficulties related to well events is worth investigating.

Recent interest in machine learning applied to the prediction of physical processes has fueled the development of
92

93

94

"Physics Informed Deep Learning", wherein machine learning models either substitute or enhance traditional numerical algorithms while preserving the inherent constraints from the physical model. These models can be trained in a supervised or unsupervised manner.

95

96

97

98

99

100

In unsupervised learning, no labeled data is available and the objective is to directly enforce physical constraints. This can be achieved, for instance, by minimizing the residual of the partial differential equations governing the evolution of the solution or by penalizing deviations to mass conservation, etc. Physics Informed Neural Networks (PINN), as introduced by Lagaris et al. (1998) typically is an example of unsupervised learning applied to the resolution of differential equations. More recently, Raissi et al. (2019) remastered the PINN approach in a semi-supervised manner in order to solve PDE problems.

101

102

103

104

105

106

107

108

109

110

In supervised learning, the objective is to match the labeled data available from experiment or previous simulations. In particular, the supervised learning of mappings between infinite dimensional function space through neural operators ((Kovachki et al., 2021; Lu et al., 2019)) has yielded notable outcomes in learning PDE solutions (Takamoto et al. (2022)). Especially Fourier Neural Operator (FNO) Li et al. (2020) has demonstrated great performances and a lot of recent architectures based on it have emerged (Wen et al., 2022; Li et al., 2021; Raonić et al., 2023)). Physics Informed Deep Learning has been applied to geological sequestration of CO₂ using a wide range of machine learning techniques (Amini and Mohaghegh, 2019; You et al., 2020; Chu et al., 2023; Chen et al., 2018; Tariq et al., 2023). Numerous numerical simulations are required to validate an injection scenario. These simulations are costly, and machine learning has the potential to accelerate them, facilitating decision-making and enabling a faster and scalable deployment of CCS.

111

112

113

114

115

116

For example, in Chu et al. (2023), simulations were accelerated by a factor of 8000 compared to the traditional numerical solvers, while maintaining notable accuracy. This was achieved using a novel architecture based on FNO. In Wen et al. (2023), yet another innovative architecture, also based on the FNO achieves a remarkable speedup compared to existing methods in the numerical simulation of CO₂ storage. Several other notable articles also demonstrate a substantial reduction in the computational cost of CO₂ storage numerical simulations (Witte et al., 2023, 2022; Tang et al., 2023; Jiang et al., 2023; Yan et al., 2022; Shokouhi et al., 2021; Wu et al., 2021; Wen et al., 2021).

117

118

119

120

121

122

123

In most of the cases, the resulting machine learning model is used as a surrogate model that supplant the traditional solver. This allows for faster evaluation of new input parameters in contrast to the original numerical solver. This becomes especially valuable when the numerical solver is computationally expensive and time-consuming, as surrogate models can provide quick approximations. However, this 'black-box' approach lacks a guarantee of predicting solutions accurately, in contrast to the reliability of a numerical solver. Consequently, hybrid approaches that complement the numerical solver with a machine learning model while preserving the numerical guarantees are worth considering for evaluation.
We propose a proof of concept methodology to alleviate the impact of well events during the numerical simulation of CO₂ storage in the subsurface. First, we develop an open-source framework that allows to run reservoir simulation and include machine learning models at a low-level; which, to the best of our knowledge, is not freely available in the literature. The physical model that we implemented is simple (incompressible two-phase flow), yet exhibits typical numerical issues encountered during reservoir simulation. Moreover, we complement the fully implicit solver by predicting an initialization of Newton's method directly in the domain of quadratic convergence using a supervised learning approach based on recently developed Fourier Neural Operators. This approach can be seen as nonlinear preconditioning/elimination with a prediction-correction strategy. Our results show a significant decrease in the number of Newton iterations required for convergence, while ensuring the convergence to the correct solution and no time-step reduction. Moreover, this type of preconditioning may allow to run simulations with larger time-steps compared to the standard strategy which may lead to significant accelerations. In practice, this approach requires further development to be deployed in cutting-edge reservoir simulation software such as OPM Flow (Rasmussen et al., 2021), JutulDarcy (Møyner, 2024) or MRST (Lie, 2019). It requires the integration of state-of-the-art neural operator architectures within the nonlinear solver of a low-level code base, such as C++, and must also be applicable to more realistic physical models and meshes. Nonetheless, we believe it can serve as a stepping stone for machine-learning-based nonlinear preconditioning.

This paper is organized into 4 parts. In the first part (2), we introduce the mathematical problem, its numerical resolution, and an example of numerical behavior related to well events. In the second part (3), we detail the construction of the Hybrid Newton's method initialization guess. Then, in the third section (4), we detail the test cases on which we apply the methodology and the architectures of the neural networks. Finally, in the last part (5), we present the results, first on a single example, then on the entire dataset.

# 2. Continuous model and numerical resolution

## 2.1. Continuous model

Let us consider a two-phase fluid composed of an aqueous phase noted $w$ and a gaseous phase noted $g$ flowing in a porous medium. We are interested in determining the evolution of the pressure $P_\alpha$ and saturation $S_\alpha$ of each phase.

The standard approach of modeling a two-phase flow in a porous medium is based on the application of Darcy's law (Hubbert (1956)) for each phase $\alpha \in \{w, g\}$ that relates the phase velocity $v_\alpha$ with the pressure gradient:

$$v_{\alpha} = -\frac{\overline{K}k_{r_{\alpha}}(S_{\alpha})}{\mu_{\alpha}}\nabla P_{\alpha}, \qquad (1)$$

where $\overline{K}$ is the permeability tensor, $k_{r_\alpha}$ is the relative permeability that depends on the saturation $S_\alpha$, and $\mu_\alpha$ is the
dynamic viscosity. We consider the medium as isotropic, hence the tensor of permeability $\overline{K}$ can be considered as a scalar field $K$. We consider reservoir geometries such that gravity can be neglected. Capillary effects are neglected, therefore $P_w = P_g = P$. Finally, we use a quadratic relative permeability law, i.e $kr_\alpha(S_\alpha) = S_\alpha^2$. For incompressible fluids, the conservation equations describing the evolution of the pressure $P$ and saturations $S_w$ and $S_g$ are then:

$$ \phi \frac{\partial}{\partial t}(S_w) + \text{div}(v_w) = q_w, \quad \phi \frac{\partial}{\partial t}(S_g) + \text{div}(v_g) = q_g, \qquad (2) $$

$$ v_w = -\frac{K k r_w(S_w)}{\mu_w} \nabla P, \quad v_g = -\frac{K k r_g(S_g)}{\mu_g} \nabla P, \qquad (3) $$

$$ S_g + S_w = 1, \qquad (4) $$

where $\phi$ is the medium porosity field and $q_g$ and $q_w$ are source/sink terms representing the gas and water flow rates injected into the reservoir by wells, respectively. We suppose that the permeability $K$ and the porosity $\phi$ may vary in space but are constant through time. The viscosities $\mu_\alpha$ are constant.

## 2.2. Numerical resolution

The continuous model is discretized using a two-points finite volume spatial scheme on a cartesian mesh and an implicit Euler time scheme. The pressure and saturation are solved simultaneously on a bounded domain through a fully implicit scheme (see appendix 8.1, equation (10)).

We choose $P$ and $S = S_g$ as primary unknowns and $S_w = 1 - S$ as secondary unknown.

After discretization, described in appendix 8.1, for each cell $k$, the system (2) (3) (4) can be written in terms of residuals as:

$$ \begin{cases} R_{g_k}(P^{n+1}, S^{n+1}) = 0, \\ R_{w_k}(P^{n+1}, 1 - S^{n+1}) = 0. \end{cases} \qquad (5) $$

Then, the resulting non-linear system of equations is solved using Newton's method as detailed in Algorithm 1. In practice we select the maximum variation of saturation $\delta S_{max} = 0.2$, which is a common choice for Appleyard
damping, and clipping to ensure that the saturation remain physical i.e between 0 and 1. Regarding the stopping criterion, at the Newton iteration number *I*, the instantaneous balance for a phase *α* in a cell *k* is:

$$B_{\alpha,k}^{I} = |R_{\alpha,k}^{I}| \frac{dt}{|k|}, \qquad (6)$$

where $|k|$ aims to put into perspective the error between the large and small cells. Then, the stopping criterion is:

$$\max_{\alpha \in \{w,g\}} \max_{k=1,N} B_{\alpha,k}^{I} < \epsilon, \qquad (7)$$

with $\epsilon$ the stopping criterion value and *N* the total number of cells. This criterion can be seen as a test on the residual infinity norm. In this work, we select $\epsilon = 10^{-6}$.

Finally, we wished to use a framework that would allow us to run reservoir simulation and include state-of-the-art machine learning algorithms but we could not find any framework meeting this requirements. Therefore we develop our own open-source framework, YADS: Yet Another Darcy Solver (https://github.com/PINN-Well-opening-and-closing-Yads) python reservoir simulation library to run the simulations. This solution is not ideal, and the methodology needs to be tested in state-of-the-art open-source numerical solvers such as OPM Flow, the MATLAB Reservoir Simulation Toolbox (MRST), or JutulDarcy. This work humbly assesses the potential of the methodology within a Python environment that supports state-of-the-art neural network architectures.

### Newton initialization

In the 'standard' or 'classic' Newton's method, the initialization for step *n* + 1 is derived from the solution $X_{sol}^n$ obtained in the preceding step *n*. Therefore, we denote $X_{init}^{n+1} = X_{sol}^n$ the standard initialization of the step *n* + 1.

## 2.3. Impact of well events during the numerical resolution

Well events generate pressure and saturation discontinuities (in time) inside the reservoir and lead to nonlinear convergence problems as they act as singular point sources that are tightly coupled to the reservoir model Ahmed et al. (2022). These discontinuities often prevent Newton's method from converging while attempting to solve the system as the initial guess may be far from the solution. In such situations, we attempt to solve the system again with a reduced time step and continue this process until convergence is achieved. Well events can thus account for a significant portion of the actual simulation time. An illustrative scenario, inspired by the test case 2 discussed in Section 4, is depicted in Figure 1. This scenario involves a simulation of 40 time steps, each spanning 2 years. At 4 years, an injection well is opened using an well injection flow rate of $10^{-3.1} s^{-1}$ and closed at 24 years. We observe that immediately following the well opening, Newton's method necessitates 40 iterations to reach convergence, with a minimum of 5 iterations required at each subsequent step throughout the entire well opening process. In contrast, in the absence of well events,
**Algorithm 1** Newton algorithm for Fully Implicit scheme

**Inputs:**

$(P^n S^n)$: Actual Pressure and gas saturation in the reservoir

$\delta S_{max}$: Maximum saturation variation per Newton iteration

dt: Time-step

NMAXITER: Maximum number of Newton iterations

**Outputs:**

$(P^{n+1} S^{n+1})$: Solution Pressure and gas saturation in the reservoir

**Initialization:** $P^{I=0} = P^n, S^{I=0} = S^n$ (and $S_{w}^{I=0} = S_{w}^{n}$)

**Iteration:**

for $I \leftarrow 1$ to NMAXITER do

Solve the linear system:

$$ \left( \begin{matrix} \frac{\partial R_g^I}{\partial P} & \frac{\partial R_g^I}{\partial S} \\ \frac{\partial R_w^I}{\partial P} & \frac{\partial R_w^I}{\partial S} \end{matrix} \right) \left( \begin{matrix} \delta P \\ \delta S \end{matrix} \right) = - \left( \begin{matrix} R_g(P^I, S^I) \\ R_w(P^I, 1-S^I) \end{matrix} \right) \quad (8) $$

Update unknowns

$$ \alpha = \min(1, \delta S_{max}/\max(\delta S)) \triangleright \text{Compute relaxation/damping coefficient} $$

$$ P^{I+1} = P^I + \alpha \delta P, \qquad S^{I+1} = S^I + \alpha \delta S \quad (\text{and } S_{w}^{I+1} = 1 - S^{I+1}) $$

$$ S^{I+1} = \text{clipping}(S^{I+1}, \min = 0, \max = 1) \triangleright \text{clip the saturation to } [0, 1] $$

$$ \text{if convergence then} \quad \triangleright \text{Convergence test on the residual } L_{\infty} \text{ norm} $$

$$ P^{n+1} = P^{I+1}, \qquad S^{n+1} = S^{I+1} \quad (\text{and } S_{w}^{n+1} = 1 - S^{I+1}) $$

break

$$ \text{else if } I == \text{NMAXITER then} \quad \triangleright \text{Max number of iterations reached} $$

restart algorithm with $dt = dt/2$

end if

end for

End of Newton's iterations

the maximum number of Newton iterations is 5 immediately after well closure and 3 otherwise. Consequently, the well event-induced Newton iterations constitute the majority of the total Newton iterations throughout the entire scenario, accounting for 54%. Furthermore, if the maximum allowable iteration count were less than 40, the time step would be halved.

## 3. Methodology

We propose a methodology based on the Hybrid Newton algorithm (Aghili et al., 2024; Odot et al., 2021; Huang et al., 2020; Yahia et al., 2021; Bhattacharyya and Vyas, 2022) to alleviate the impact of well events. It consists in predicting via machine learning an initialization closer to the solution than the standard initialization, hence acceler-
ating the Newton's method convergence while still guaranteeing the correct solver solution. We use a fixed reservoir geometry and configuration where a well event occurs at a specific location. We aim to alleviate the impact of the well opening for a wide range of scenarios.

## 3.1. Hybrid Newton algorithm

The hybrid Newton algorithm is a modification of the standard Newton algorithm with respect to the initialization. Indeed, the hybrid method proposes to initialize with a more accurate solution instead of using the standard initialization (i.e. the previous time step solution). By using an initialization closer to the solution, we aim to converge in fewer iterations.

In this article, our primary focus is on the hybrid guess construction $X_{ML}^{n+1}$ and its performance compared to a reference initial guess. We will compare them through their impact on the number of Newton iterations.

## 3.2. Initial guess construction

In the hybrid Newton algorithm, a prediction guess $X_{ML}^{n+1} := (P_{pred}, S_{pred})$ is required. We propose the following construction for the pressure and saturation:

1. $P_{pred}$: Implicit Pressure solver

2. $S_{pred}$: Prediction using supervised learning

### 3.2.1. Pressure

During a well event, the pressure discontinuity is global in the reservoir and instantaneous in time. We propose to use the solution of an implicit pressure solver $P_{imp}$ (IMP) (Sheldon and Cardwell, 1959) at the step $n+1$ as a prediction guess $P_{pred}$. The implicit pressure solver solves the linear elliptic equation (9) and catches the main global variations of pressure but does not catch the small local variations of pressure due to saturation variations. Moreover, the implicit pressure solver only requires to solve a linear system. Therefore, we consider that it is not worth it to train a machine learning model as we have a good approximation. The implicit pressure solver solves the following equation:

$$div(v(P_{imp}, S^n)) = q_T, \qquad (9)$$

with $v = v_g + v_{lw}$, $v, v_g$ and $v_{lw}$ are respectively the total velocity, the gas velocity and the water velocity inside the reservoir. $q_T = q_g + q_{lw}$, $q_T, q_g, q_{lw}$ are respectively the total, gas and water injection flow rates. In practice, only gas is injected, therefore $q_T = q_g$.
### 3.2.2. Saturation

During a well event, the saturation discontinuity inside the reservoir is local and near to the well. One could use an implicit saturation solver (IMS) which would require a Newton's method to solve a non linear system smaller than one in the fully implicit formulation but which would still present the same convergence issue. Therefore, we propose to predict the saturation using a neural network considering that the inference of a neural network is rather fast compared to multiples Newton iterations of an implicit saturation solver. The predicted saturation is denoted $S_{ML}$.

The standard/classic methodology uses the initial guess $X_{init}^{n+1} = X_{sol}^{n} = (P^n, S^n)$ and the hybrid methodology uses the initial guess $X_{ML}^{n+1} = (P_{imp}, S_{ML})$. Given that we aim to conduct a fair assessment of the impact of a saturation predictive model, we compare the hybrid initialization with a reference initialization. Therefore we use in this article the following initial guess for the reference methodology $X_{init}^{n+1} = (P_{imp}, S^n)$.

In the end, we will compare in terms of Newton iterations the two following initialization:

1. Reference/standard initialization: $X_{init}^{n+1} = (P_{imp}, S^n)$

2. Hybrid initialization: $X_{ML}^{n+1} = (P_{imp}, S_{ML})$

## 3.3. Neural Network architecture

The objective is to predict using reservoir and well information the global saturation state reached after a well event. A well event at a specific location and a specific well geometry can be described with few parameters: an injection well flow $q_g$ and a time-step dt. A neural network architecture with good predictive capability for different physics-based processes is required.

### Fourier Neural Operator

Neural networks are commonly used to learn relationships between finite-dimensional spaces, but they can struggle to adapt to changes in governing equations or conditions (Fuks and Tchelepi, 2020). The Fourier Neural Operator (FNO) (Li et al., 2020) addresses this issue by learning relationships between infinite-dimensional spaces using data-driven methods. This allows the FNO to understand the rules governing an entire family of partial differential equations. Additionally, the FNO improves computational efficiency by converting convolution operations in neural networks to multiplication through the use of discrete Fourier transforms.

### Selected Architecture

We use the architecture presented in figure 2 based on an uplifting dense layer, four Fourier Layers (see figure 3) and two dense layers. We use the Gaussian Error Linear Units (GELU) (Hendrycks and Gimpel, 2016) as activation function for each layer. We use two different versions of this general architecture for the training on the two test cases. Nc and Ni are parameters respectively representing the number of channels and the number of inputs.

The architecture of a Fourier layer (fig. 3) is composed of two parts, one that apply Fourier transform, a linear
<table><caption>**Table 1**
Petrophysical Parameters</caption><thead><tr><th></th><th>Barrier zone</th><th>Drain zone</th></tr></thead><tbody><tr><td>Porosity [-]</td><td>0.2</td><td>0.2</td></tr><tr><td>Permeability[m²]</td><td>10<sup>-15</sup></td><td>100 × 10<sup>-15</sup></td></tr></tbody></table>

transform on the lower Fourier modes and a filter on the higher mode, then it applies the inverse Fourier transform. The other part is composed of a local linear transformation applied to the original input. Finally, the output of the two parts are added together and an activation function is applied.

## Input Features

There are multiple features that can be considered as input features for the selected neural architecture. We select four possible input features, the time-step dt, the injection flow rate $q_g$, the initial saturation $S^n$ and the implicit pressure $P_{imp}$. As $S^n$ and $P_{imp}$ have the same shape as the output feature, they can be used straightforward. However, as $q_g$ and dt are scalars, they need to be reshaped. We propose to reshape dt into a constant map of value dt everywhere. For the injection flow rate $q_g$, we reshape it to a map of value zero everywhere except at the well location where it takes the value $q_g$.

# 4. Test case and Database Generation

## 4.1. SHPCO2 benchmark

### Context

As a practical use case, we use the SHPCO2 (Haeberlein, 2011) benchmark to test our methods. This benchmark was created for modelling reactive transport for CO₂ geological storage.

### Reservoir configuration

We adapt the original 2D SHPCO2 benchmark by removing the gas zone and replacing it by a well at its center as depicted in figure 5.

The domain after modification is separated in two zones, the first zone called "Barrier zone" is colored in green in the figure 5 and the second zone called "Drain zone" is formed of the rest of the domain. These two domains have different petrophysical properties (see table 1).

The boundary conditions of the adapted SHPCO2 (see figure 5) are presented in table 2 with Injector 1 and Injector 2 are used to impose pressure boundary conditions.

We consider the mesh geometry corresponding to the small (S) mesh size of the SHPCO2 benchmark. It consists of a Cartesian mesh discretized into 95 cells along the x-axis and 60 cells along the y-axis, resulting in a total of 5700 cells.
<table><caption>**Table 2**
Boundary condition parameters</caption><thead><tr><th></th><th>Pressure[MPa]</th><th>Composition</th></tr></thead><tbody><tr><td>Injector 1</td><td>11</td><td>Water</td></tr><tr><td>Injector 2</td><td>10.5</td><td>Water</td></tr><tr><td>Productor</td><td>10</td><td>-</td></tr></tbody></table>

## 4.2. Test cases

We propose two different test cases based on the SHPCO2 benchmark and its reservoir configuration presented in section 4.1. The test case 1 has constant initial saturation $S_0$ while the test case 2 has more realistic initial saturations.

### 4.2.1. Test case 1

#### Database Generation

We launch simulations with a constant reservoir configuration except for three parameters, $S_0$ the initial saturation, $q_g$ the well injection flow rate and $dt$ the time-step. We allow a maximum of 200 Newton iterations to converge and a stopping criterion of $10^{-6}$.

We generate 5004 parameter combinations through a Latin Hypercube Sampling strategy Stein (1987) within the following ranges: $S_0 \in [0, 0.6]$, $|q_g| \in [10^{-5}, 10^{-3}]$ $1/s$ which corresponds to a well pressure $\in$]10, 20] MPa and $dt \in$ [0.1, 10] years. Note that $P_{imp}$ is obtained using the Implicit Pressure solver (IMP). An example of the resulting pressure and saturation fields is illustrated in Figure 6.

### 4.2.2. Test case 2

In the test case 1, we use constant reservoir saturation maps as initial saturation $S_0$ and realise one well opening with a particular well injection flow rate $q_g$ and time-step $dt$. In the test case 2, we use more realistic initial saturation maps. To do so, we realise consecutive well opening and closure events (see figure 7). We use Latin Hypercube Sampling strategy to generate parameter combinations. The initial parameters are $q_g$, $dt$ and $S^n$. After the first simulation, we close the well (i.e $q_g = 0$) and launch another simulation using the previous reservoir state obtained and a new time-step $dt$ (sampled through Latin Hypercube Sampling). Finally we open the well and launch a simulation with $q_g$ and $dt$ as parameters. The opening and closure step are repeated as many times as needed (see figure 7). The reservoir state is saved at every well opening or closure.

We generate data with 5 well openings and 4 closures and 3600 parameter combinations for each step. The parameters are sampled using a Latin Hypercube Sampling strategy within the following ranges: $S_0 \in [0, 0.6]$, $q_g \in [10^{-3}, 10^{-5}]$ $1/s$ and $dt \in$ [1, 10] years in seconds.

We perform simulations for each parameter combination and we only consider data where the well injected flow rate $q_g$ is not null. We have 5 wells openings (i.e $q_g > 0$) and 3600 parameter combinations for each one. Therefore,
there is a total of 18000 samples. When splitting the samples in train and test sets, data coming from a same scenario are sent together in a set (i.e 5 per 5). An example of the resulting pressure and saturation fields is illustrated in Figure 8.

# 5. Results

## 5.1. Neural Network training

### 5.1.1. Test case 1

We use the Neural Network architecture presented in figure 2 with $N_c = 32$ channels in the Fourier layers and $\{q_g, dt, S^n\}$ as inputs (i.e $N_i = 3$). The input parameter $dt$ is a scalar, therefore, we reshape it in a constant map of shape (95, 60). Moreover, $q_g$ is also a scalar. We reshape it in a (95,60) map which values are zeroes everywhere except at the well location where it takes $q_g$ as value.

We split the data in a train and test sets with a 80/20 splitting ratio. We train the model on a NVIDIA V100 GPU during 27 hours using Adam optimizer, a batch size of 10, a learning rate of $5 \times 10^{-5}$, a momentum of 0.9, a weight decay of $10^{-4}$ and keep the model parameters corresponding to the lowest test loss value. We show the relative L2 loss evolution in the figure 9. The lowest test loss value is $1.9 \times 10^{-3}$ reached at epoch 17285. The corresponding train loss value is $2.0 \times 10^{-3}$.

### 5.1.2. Test case 2

We use the Neural Network architecture presented in figure 2 with $N_c = 64$ channels in the Fourier layers. As the case is more complex than the test case 1, the neural network is harder. To alleviate this complexity, the implicit pressure is added to the input features. We therefore use $\{P_{imp}, q_g, dt, S^n\}$ as input features (i.e $N_i = 4$). The input parameter $dt$ is a scalar, therefore, we reshape it in a constant map of shape (95, 60). Moreover, $q_g$ is also a scalar. We reshape it in a (95,60) map which values are zeroes everywhere except at the well location where it takes $q_g$ as value. $S^n$ and $P_{imp}$ can be used straightforward.

We split the data in a train and test set with a 80/20 splitting ratio. We train the model on a NVIDIA V100 GPU during 132 hours using Adam optimizer, a batch size of 128, a momentum of 0.9, a weight decay of $10^{-4}$ and keep the model parameters corresponding to the lowest test loss value. We start with a learning rate of $10^{-4}$, at iteration number 1000, we decrease it to $5 \times 10^{-5}$, then $10^{-5}$ at iteration number 1600 and finally we set the learning rate to $5 \times 10^{-7}$ at iteration number 3600 and until the end.

We show the L2 loss evolution in the figure 10. The lowest test loss value is $8.7 \times 10^{-4}$ reached at epoch 7295. The corresponding train loss value is $9.3 \times 10^{-4}$.
## 5.2. Results

In this section, we present the results obtained using the hybrid methodology, first on the scenario presented in the figure 1 and then on the two tests cases. We compare the performances in term of Newton iterations obtained between the standard and hybrid approaches. We chose the number of Newton iterations as a metric because measuring CPU time for each method using the under-optimized Python framework we developed (YADS) would not accurately reflect the performance of an industrial solver. However, the number of Newton iterations can be compared with that of industrial software, providing a more meaningful benchmark.

### 5.2.1. Single prediction example

We apply the hybrid methodology on the example scenario presented in figure 1 only during the well event, this implies that every step outside of the well event requires the same number of Newton iterations for the standard and hybrid approaches. The results are shown in figure 11.

We observe a significant reduction in the number of Newton iterations during the well event. Immediately after the well opening, the standard approach (depicted by blue crosses) requires 40 iterations, whereas the hybrid approach (represented by green crosses) only necessitates 3 iterations. Furthermore, throughout the entire well event, the hybrid approach consistently requires fewer Newton iterations compared to the standard approach, resulting in a 68% reduction in Newton iterations during the well event and a 38% reduction across the entire scenario.

### 5.2.2. Test case 1

We launch simulations with the same parameters combinations of test case 1, $X_{init}^{n+1} = (P_{imp}, S^n)$ and $X_{ML}^{n+1} = (P_{imp}, S_{ML})$ as initial guesses. $P_{imp}$ is calculated using the Implicit Pressure Solver and $S_{ML}$ using the model obtained in the previous section. The results are presented in figure 12.

We observe that the hybrid methodology facilitates Newton's method initialization directly within the domain of quadratic convergence, resulting in a maximum of 5 Newton iterations for the training set and 4 for the test set. In essence, the more challenging the problem using the standard methodology, the greater the potential benefit from employing the hybrid methodology. With the hybrid formulation, we achieve an average speedup of 54%, translating to 54% fewer Newton iterations compared to the standard methodology for the training set and a 53% reduction for the test set.

### 5.2.3. Test case 2

We launch simulations with the same parameters combinations of test case 2, $X_{init}^{n+1} = (P_{imp}, S^n)$ and $X_{ML}^{n+1} = (P_{imp}, S_{ML})$ as initial guesses. $P_{imp}$ is calculated using the Implicit Pressure Solver and $S_{ML}$ using the model obtained in the previous section. The results are presented in figure 13.
Considering every simulations and using the hybrid methodology, we speed up by 39%, i.e. 39% less Newton iterations than the standard methodology the computations for the training set and by 38.7% for the test set.

## 5.3. Discussion

Considering a wide range of injection scenarios, we show that it is possible to learn the impact of a well event on a reservoir. We speed up by 53% the handling of well events for the test case 1 and by 38% for the test case 2. Moreover the hybrid Newton methodology is quite general and can be applied to any problem that requires an important number of costly iterations. Finally, we observe that the more challenging the problem is using the standard methodology, the greater the potential benefit from employing the hybrid methodology. However, there are some limiting issues that needs to be considered.

### 5.3.1. Offline cost

The data generation and training phase of this study are resource-intensive. For instance, Test Case 2 necessitated 132 hours of training and a substantial volume of full reservoir simulations for data generation. While this configuration is suitable for proof of concept purposes, it lacks scalability for direct integration into industrial software. Future efforts should focus on reducing these offline costs, such as by shortening training time or lowering dataset generation costs.

### 5.3.2. Constant discretization

We use a specific discretization (SHPCO2 S mesh) for the data generation and we predict on the same discretization. This implies that the model would not work for different meshes. A new model has to be generated. An extension of this work could involve testing it on 3D meshes and/or non-cartesian meshes. Moreover, the idea of a model invariant to discretization is developed through Neural Operators (Kovachki et al., 2021). A model could be trained using data from different discretizations and predict the solution on multiple discretizations.

### 5.3.3. Constant well position

The methodology is applied on a constant grid with a constant well position. While the pressure variations are global during a well event, the saturation variations are local. Therefore, if we change the well position, the model prediction will not be accurate. To alleviate this issue, a local approach could be used, i.e. create a model that predicts the saturation only around a well.

## 6. Conclusion and perspectives

We proposed in this article a proof of concept methodology to alleviate the impact of well events during the numerical simulation of CO₂ storage in the subsurface. We complement the standard numerical algorithm by predicting an initialization of Newton's method directly in the domain of quadratic convergence using a supervised learning ap-
398

399

400

401

402

proach based on recently developed Fourier Neural Operators. Our results show a significant decrease in the number of Newton iterations required for convergence, while ensuring the convergence to the correct solution. Moreover the hybrid Newton methodology is quite general and can be applied as non-linear preconditioning to any problem that requires an important number of costly iterations. Finally, we observe that the more challenging the problem is using the standard methodology, the greater the potential benefit from utilizing the hybrid methodology.

403

404

405

406

407

408

Reservoir engineers seek the best location for a new well. Therefore, enhancing the computational efficiency of managing well events across a wide range of positions presents an intriguing opportunity. To achieve this, a local approach could be adopted, where a model is developed to predict saturation specifically within the near-well region. Additionally, this local approach offers the advantage of reducing training costs. Furthermore, reservoir engineers perform history matching and run simulations to forecast optimal well conditions. In the context of a fixed well location, the global hybrid Newton method can be applied and prove valuable.

409

# 7. Acknowledgments

410

411

The authors thank Sorbonne University, IFP Energies nouvelles and DIM Math Innov for supporting financially and physically this research work developed as part of Antoine Lechevallier Ph.D thesis work.
Code availability section

Name of the code/library: YADS: Yet Another Darcy Solver

Contact: e-mail and phone number anle@norceresearch.no

Program language: Python

Software required: Paraview for visualization (Optional)

The source codes are available for downloading at the link: https://github.com/PINN-Well-opening-and-closing-events/Yads

References

Aghili, J., Franck, E., Hild, R., Michel-Dansac, V., Vigon, V., 2024. Accelerating the convergence of newton's method for nonlinear elliptic pdes using fourier neural operators. arXiv:2403.03021.

Ahmed, E., Klemetsdal, Raynaud, X., Møyner, O., Nilsen, H.M., 2022. Adaptive Timestepping, Linearization, and A Posteriori Error Control for Multiphase Flow of Immiscible Fluids in Porous Media with Wells. SPE Journal, 1-21URL: https://doi.org/10.2118/203974-PA, doi:10.2118/203974-PA, arXiv:https://onepetro.org/SJ/article-pdf/doi/10.2118/203974-PA/3022338/spe-203974-pa.pdf.

Amini, S., Mohaghegh, S., 2019. Application of machine learning and artificial intelligence in proxy modeling for fluid flow in porous media. Fluids 4. URL: https://www.mdpi.com/2311-5521/4/3/126, doi:10.3390/fluids4030126.

Bhattacharyya, S., Vyas, A., 2022. A novel methodology for fast reservoir simulation of single-phase gas reservoirs using machine learning. Heliyon 8, e12067. URL: https://www.sciencedirect.com/science/article/pii/S2405844022033552, doi:https://doi.org/10.1016/j.heliyon.2022.e12067.

Cai, X.c., Keyes, D., 2001. Nonlinearly preconditioned inexact newton algorithms. SIAM Journal on Scientific Computing 24. doi:10.1137/S106482750037620X.

Cai, X.C., Keyes, D.E., Marcinkowski, L., 2002. Non-linear additive schwarz preconditioners and application in computational fluid dynamics. International Journal for Numerical Methods in Fluids 40, 1463-1470. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.404, doi:https://doi.org/10.1002/fld.404, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.404.

Chen, B., Harp, D.R., Lin, Y., Keating, E.H., Pawar, R.J., 2018. Geologic CO₂ sequestration monitoring design: A machine learning and uncertainty quantification based approach. Applied Energy 225, 332-345. URL: https://ideas.repec.org/a/eee/appene/v225y2018icp332-345.html, doi:10.1016/j.apenergy.2018.0.

Chen, Z., Huan, G., Li, B., 2004. An improved impes method for two-phase flow in porous media. Transport in Porous Media 54, 361-376. doi:10.1023/B:TIPM.0000003667.86625.15.

Chen, Z., Huan, G., Ma, Y., 2006. Computational Methods for Multiphase Flows in Porous Media. volume 2 of Computational science and engineering. Society for Industrial and Applied Mathematics. URL: https://www.worldcat.org/title/computational-methods-for-multiphase-flows-in-porous-media/oclc/7346429515&referrer=brief_results.

Chu, A.K., Benson, S.M., Wen, G., 2023. Deep-learning-based flow prediction for co2 storage in shale&ndash;sandstone formations. Energies 16. URL: https://www.mdpi.com/1996-1073/16/1/246, doi:10.3390/en16010246.

Dawson, C., Klie, H., Wheeler, M., Woodward, C., 1997. A parallel, implicit, cell-centered method for two-phase flow with a preconditioned newton-krylov solver. Computational Geosciences 1, 215-249. doi:10.1023/A:1011521413158.
Eymard, R., Thierry, G., Herbin, R., 2000. Finite volume methods 7. doi:10.1016/S1570-8659(00)07005-8.

Fuks, O., Tchelepi, H., 2020. Limitations of physics informed machine learning for nonlinear two-phase transport in porous media.

Haeberlein, F., 2011. Time Space Domain Decomposition Methods for Reactive Transport — Application to CO₂ Geological Storage. Theses.
Université Paris-Nord - Paris XIII. URL: https://tel.archives-ouvertes.fr/tel-00634507.

Hendrycks, D., Gimpel, K., 2016. Gaussian error linear units (gelus). URL: https://arxiv.org/abs/1606.08415, doi:10.48550/ARXIV.1606.08415.

Huang, J., Wang, H., Yang, H., 2020. Int-deep: A deep learning initialized iterative method for nonlinear problems. Journal of Computational Physics 419, 109675. URL: https://www.sciencedirect.com/science/article/pii/S0021999120304496, doi:https://doi.org/10.1016/j.jcp.2020.109675.

Hubbert, M.K., 1956. Darcy's Law and the Field Equations of the Flow of Underground Fluids. Transactions of the AIME 207, 222-239. URL: https://doi.org/10.2118/749-G, doi:10.2118/749-G, arXiv:https://onepetro.org/TRANS/article-pdf/207/01/222/2176441/spe-749-g.pdf.

Hwang, F.N., Cai, X.C., 2005. A parallel additive schwarz preconditioned inexact newton algorithm for incompressible navier-stokes equations. Journal of Computational Physics 204, 666-691. doi:10.1016/j.jcp.2004.10.025.

Jiang, Z., Zhu, M., Li, D., Li, Q., Yuan, Y.O., Lu, L., 2023. Fourier-mionet: Fourier-enhanced multiple-input neural operators for multiphase modeling of geological carbon sequestration. arXiv:2303.04778.

Klemetsdal, S., Møyner, O., Lie, K.A., 2019. Robust Nonlinear Newton Solver With Adaptive Interface-Localized Trust Regions. SPE Journal 24, 1576-1594. URL: https://doi.org/10.2118/195682-PA, doi:10.2118/195682-PA. _eprint: https://onepetro.org/SJ/article-pdf/24/04/1576/2123843/spe-195682-pa.pdf.

Kovachki, N., Li, Z., Liu, B., Azizzadenesheli, K., Bhattacharya, K., Stuart, A., Anandkumar, A., 2021. Neural operator: Learning maps between function spaces. URL: https://arxiv.org/abs/2108.08481, doi:10.48550/ARXIV.2108.08481.

Lagaris, I., Likas, A., Fotiadis, D., 1998. Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Neural Networks 9, 987-1000. URL: https://doi.org/10.1109/2F72.712178, doi:10.1109/72.712178.

LeVeque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge Texts in Applied Mathematics, Cambridge University Press. doi:10.1017/CBO9780511791253.

Li, Z., Kovachki, N.B., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A.M., Anandkumar, A., 2020. Fourier neural operator for parametric partial differential equations. CoRR abs/2010.08895. URL: https://arxiv.org/abs/2010.08895, arXiv:2010.08895.

Li, Z., Zheng, H., Kovachki, N.B., Jin, D., Chen, H., Liu, B., Azizzadenesheli, K., Anandkumar, A., 2021. Physics-informed neural operator for learning partial differential equations. CoRR abs/2111.03794. URL: https://arxiv.org/abs/2111.03794, arXiv:2111.03794.

Lie, K.A., 2019. An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User Guide for the MATLAB Reservoir Simulation Toolbox (MRST). Cambridge University Press.

Liu, L., Keyes, D., 2015. Field-split preconditioned inexact newton algorithms. SIAM Journal on Scientific Computing 37, 1388-1409. doi:10.1137/140970379.

Liu, L., Keyes, D., Sun, S., 2013. Fully implicit two-phase reservoir simulation with the additive schwarz preconditioned inexact newton method, in: SPE Reservoir Characterization and Simulation Conference and Exhibition, Society of Petroleum Engineers (SPE). doi:10.2118/166062-ms. kAUST Repository Item: Exported on 2020-10-01.

Lu, L., Jin, P., Karniadakis, G.E., 2019. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. CoRR abs/1910.03193. URL: http://arxiv.org/abs/1910.03193, arXiv:1910.03193.
486 Møyner, O., 2024. Jutuldarcy.jl - a fully differentiable high-performance reservoir simulator based on automatic differentiation 2024, 1-9.
doi:https://doi.org/10.3997/2214-4609.202437111.

487 Odot, A., Haferssas, R., Cotin, S., 2021. Deepphysics: a physics aware deep learning framework for real-time simulation. CoRR abs/2109.09491.
URL: https://arxiv.org/abs/2109.09491, arXiv:2109.09491.

489 Ouaki, F., 2013. Etude de schémas multi-échelles pour la simulation de réservoir. Ph.D. thesis. Ecole Polytechnique X.

491 Palomino Monteagudo, J., Firoozabadi, A., 2007. Comparison of fully implicit and impes formulations for simulation of water injection in fractured and unfractured media. International Journal for Numerical Methods in Engineering 69, 698-728. doi:10.1002/nme.1783.

493 Peaceman, D., 1978. Interpretation of Well-Block Pressures in Numerical Reservoir Simulation (includes associated paper 6988).
Society of Petroleum Engineers Journal 18, 183-194. URL: https://doi.org/10.2118/6893-PA, doi:10.2118/6893-PA,
arXiv:https://onepetro.org/spejournal/article-pdf/18/03/183/2158027/spe-6893-pa.pdf.

496 Pour, K.M., Voskov, D., Bruhn, D., 2023. Nonlinear solver based on trust region approximation for CO2 utilization and storage in subsurface reservoir. Geoenergy Science and Engineering 225, 211698. URL: https://www.sciencedirect.com/science/article/pii/S2949891023002853, doi:https://doi.org/10.1016/j.geoen.2023.211698.

499 Quandalle, P., Savary, D., 1989. An implicit in pressure and saturations approach to fully compositional simulation, in: SPE Reservoir Simulation Conference, SPE. pp. SPE-18423.

501 Raissi, M., Perdikaris, P., Karniadakis, G., 2019. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, 686-707. URL: https://www.sciencedirect.com/science/article/pii/S0021999118307125, doi:https://doi.org/10.1016/j.jcp.2018.10.045.

504 Raonić, B., Molinaro, R., Ryck, T.D., Rohner, T., Bartolucci, F., Alaifari, R., Mishra, S., de Bézenac, E., 2023. Convolutional neural operators for robust and accurate learning of pdes. arXiv:2302.01178.

506 Rasmussen, A.F., Sandve, T.H., Bao, K., Lauser, A., Hove, J., Skaflestad, B., Klöfkorn, R., Blatt, M., Rustad, A.B., Sævareid, O., Lie, K.A., Thune, A., 2021. The open porous media flow reservoir simulator. Computers & Mathematics with Applications 81, 159-185. URL: https://www.sciencedirect.com/science/article/pii/S0898122120302182, doi:https://doi.org/10.1016/j.camwa.2020.05.014. development and Application of Open-source Software for Problems with Numerical PDEs.

510 Schlumberger, 2013. ECLIPSE: Technical Description 2013.2. pp. 391-391.

511 Sheldon, J.W., Cardwell, W.T.J., 1959. One-dimensional, incompressible, noncapillary, two-phase fluid flow in a porous medium. Transactions of the AIME 216, 290-296.

513 Shokouhi, P., Kumar, V., Prathipati, S., Hosseini, S.A., Giles, C.L., Kifer, D., 2021. Physics-informed deep learning for prediction of co2 stor-
age site response. Journal of Contaminant Hydrology 241, 103835. URL: https://www.sciencedirect.com/science/article/pii/
S0169772221000747, doi:https://doi.org/10.1016/j.jconhyd.2021.103835.

516 Skogestad, J.O., Keilegavlen, E., Nordbotten, J., 2013. Domain decomposition strategies for nonlinear flow problems in porous media. Journal of Computational Physics 234, 439-451. doi:10.1016/j.jcp.2012.10.001.

518 Stein, M., 1987. Large sample properties of simulations using latin hypercube sampling. Technometrics 29, 143-151. URL:
https://www.tandfonline.com/doi/abs/10.1080/00401706.1987.10488205, doi:10.1080/00401706.1987.10488205,
arXiv:https://www.tandfonline.com/doi/pdf/10.1080/00401706.1987.10488205.

521 Takamoto, M., Praditia, T., Leiteritz, R., MacKinlay, D., Alesiani, F., Pflüger, D., Niepert, M., 2022. Pdebench: An extensive benchmark for scientific machine learning, in: Advances in Neural Information Processing Systems, Curran Associates, Inc., pp. 1596-1611. URL: https://proceedings.neurips.cc/paper_files/paper/2022/file/0a9747136d411fb83f0cf81820d44afb-Paper-Datasets_
and_Benchmarks.pdf.

Tang, H., Kong, Q., Morris, J.P., 2023. Multi-fidelity fourier neural operator for fast modeling of large-scale geological carbon storage. arXiv:2308.09113.

Tariq, Z., Yildirim, E.U., Gudala, M., Yan, B., Sun, S., Hoteit, H., 2023. Spatial-temporal prediction of minerals dissolution and precipitation using deep learning techniques: An implication to geological carbon sequestration. Fuel 341, 127677. URL: https://www.sciencedirect.com/science/article/pii/S0016236123002909, doi:https://doi.org/10.1016/j.fuel.2023.127677.

Wen, G., Hay, C., Benson, S.M., 2021. Ccsnet: A deep learning modeling suite for co2 storage. Advances in Water Resources 155, 104009. URL: https://www.sciencedirect.com/science/article/pii/S0309170821001640, doi:https://doi.org/10.1016/j.advwatres.2021.104009.

Wen, G., Li, Z., Azizzadenesheli, K., Anandkumar, A., Benson, S.M., 2022. U-fno—an enhanced fourier neural operator-based deep-learning model for multiphase flow. Advances in Water Resources 163, 104180. URL: https://www.sciencedirect.com/science/article/pii/S0309170822000562, doi:https://doi.org/10.1016/j.advwatres.2022.104180.

Wen, G., Li, Z., Long, Q., Azizzadenesheli, K., Anandkumar, A., Benson, S.M., 2023. Real-time high-resolution co2 geological storage prediction using nested fourier neural operators. Energy Environ. Sci. 16, 1732–1741. URL: http://dx.doi.org/10.1039/D2EE04204E, doi:10.1039/D2EE04204E.

Witte, P.A., Hewett, R., Chandra, R., 2022. Industry-scale co2 flow simulations with model-parallel fourier neural operators, in: NeurIPS 2022 Workshop on Tackling Climate Change with Machine Learning. URL: https://www.climatechange.ai/papers/neurips2022/78.

Witte, P.A., Konuk, T., Skjetne, E., Chandra, R., 2023. Fast co2 saturation simulations on large-scale geomodels with artificial intelligence-based wavelet neural operators. International Journal of Greenhouse Gas Control 126, 103880. URL: https://www.sciencedirect.com/science/article/pii/S1750583623000506, doi:https://doi.org/10.1016/j.ijggc.2023.103880.

Wu, H., Lubbers, N., Viswanathan, H.S., Pollyea, R.M., 2021. A multi-dimensional parametric study of variability in multi-phase flow dynamics during geologic co2 sequestration accelerated with machine learning. Applied Energy 287, 116580. URL: https://www.sciencedirect.com/science/article/pii/S0306261921001264, doi:https://doi.org/10.1016/j.apenergy.2021.116580.

Yahia, I.B., Merlet, J.P., Papegay, Y., 2021. Mixing neural networks and the Newton method for the kinematics of simple cable-driven parallel robots with sagging cables, in: ICAR 2021 - 20th International Conference on advanced robotics, Ljulbjana, Slovenia. URL: https://hal.inria.fr/hal-03385003, doi:10.1109/ICAR53236.2021.9659400.

Yan, B., Chen, B., Robert Harp, D., Jia, W., Pawar, R.J., 2022. A robust deep learning workflow to predict multiphase flow behavior during geological co2 sequestration injection and post-injection periods. Journal of Hydrology 607, 127542. URL: https://www.sciencedirect.com/science/article/pii/S0022169422001172, doi:https://doi.org/10.1016/j.jhydrol.2022.127542.

Yang, H., Yang, C., Sun, S., 2016. Active-set reduced-space methods with nonlinear elimination for two-phase flow problems in porous media. SIAM Journal on Scientific Computing 38. doi:10.1137/15M1041882.

You, J., Ampomah, W., Sun, Q., 2020. Development and application of a machine learning based multi-objective optimization workflow for co2-eor projects. Fuel 264, 116758. URL: https://www.sciencedirect.com/science/article/pii/S001623611932112X, doi:https://doi.org/10.1016/j.fuel.2019.116758.
# 8. Appendix

## 8.1. Mesh definition and Fully Implicit discretization

### Mesh definition and notations

Let $\mathcal{K}_h$ a finite set of non empty, open and disjoint polyhedra forming a partition of $\Omega$ such that $h = \max_{k \in \mathcal{K}_h} h_k$ where $h_k$ is the diameter of element $k \in \mathcal{K}_h$. The set $\mathcal{K}_h$ thus forms a mesh of $\Omega$ whose elements $k \in \mathcal{K}_h$ are the cells. We define $N$ as the number of cells in the mesh $\mathcal{K}_h$. Similarly, we denote $\mathcal{F}_h$ the set of faces in the mesh. Moreover, we define $\mathcal{F}_h^i$ the set of internal faces of the mesh

$$ \mathcal{F}_h^i = \{\sigma \in \mathcal{F}_h | \sigma \not\subset \partial\Omega.\} $$

We denote $\mathcal{F}_h^k$ the set of boundary faces

$$ \mathcal{F}_h^k = \{\sigma \in \mathcal{F}_h | \sigma \subset \partial\Omega\}. $$

The set of faces $\mathcal{F}_h$ of the mesh $\mathcal{K}_h$ then verifies $\mathcal{F}_h = \mathcal{F}_h^i \cup \mathcal{F}_h^k$. We denote $\mathcal{F}_{h,k}$ the set of faces of a cell $k \in \mathcal{K}_h$:

$$ \mathcal{F}_{h,k} = \{\sigma \in \mathcal{F}_h | \sigma \subset \partial k\}. $$

Moreover, we denote $\mathcal{F}_{h,k}^i$ the set of internal faces of a cell $k \in \mathcal{K}_h$:

$$ \mathcal{F}_{h,k}^i = \mathcal{F}_{h,k} \cap \mathcal{F}_h^i, $$

and $\mathcal{F}_{h,k}^b$, the set of faces of a cell $k$ that are located on the border of the domain $\Omega$

$$ \mathcal{F}_{h,k}^b = \mathcal{F}_{h,k} \cap \mathcal{F}_h^k. $$

$\mathcal{F}_{h,k}^b = \emptyset$ if the cell $k$ does not share faces with $\partial\Omega$.

Boundary conditions are imposed on a border $\Gamma_D \subset \partial\Omega$. We then denote $\mathcal{F}_h^D$ the set of border faces included in

$\Gamma_D$:
$$ \mathcal{F}_h^D = \{ \sigma \in \mathcal{F}_h^k | \sigma \subset \Gamma_D \}, $$

and $\mathcal{F}_{h,k}^D$ the set of faces of the cell $k \in \mathcal{K}_h$, that are located on $\Gamma_D$

$$ \mathcal{F}_{h,k}^D = \mathcal{F}_{h,k} \cap \mathcal{F}_h^D. $$

We denote by $\Pi$ a set of wells, $\mathcal{K}_i^\Pi$ the set of cells perforated by a well $i$ and $\mathcal{K}_{k,i}^\Pi$, the cell $k$ perforated by a well $i \in \Pi$

$$ \mathcal{K}_{k,i}^\Pi = \{ k \in \mathcal{K}_i^\Pi | i \in \Pi \}. $$

In a similar way to $\mathcal{F}_{h,k}$, we note $\mathcal{K}_{h,\sigma}$ the set of cells $k \in \mathcal{K}_h$ located on both side of a face $\sigma \in \mathcal{F}_h$:

$$ \mathcal{K}_{h,\sigma} = \{ k, l \in \mathcal{K}_h | \sigma = \partial k \cap \partial l \}. $$

For a cell $k \in \mathcal{K}_h$, we identify its center as $x_k$. Similarly, for any face $\sigma \in \mathcal{F}_h$, $x_\sigma$ is the center of the face. For any cell $k \in \mathcal{K}_h$, and for any face $\sigma \in \mathcal{F}_{h,k}$, we note:

$$ d_{k,\sigma} = dist(x_k, x_\sigma). $$

Similarly, for any $k_1, k_2 \in \mathcal{K}_h$:

$$ d_{k_1,k_2} = dist(x_{k_1}, x_{k_2}). $$

For a face $\sigma \in \mathcal{F}_h$, $|\sigma|$ is the face area. In the same way, for a cell $k \in \mathcal{K}_h$, $|k|$ is the cell volume.

## Fully Implicit discretization

The fully implicit discretization of (2), (3), (4) gives $\forall k \in \mathcal{K}_h$ and for $\alpha \in \{w, g\}$:
$$|k|\phi_k \frac{S_{\alpha,k}^{n+1} - S_{\alpha,k}^n}{dt} + \left( \sum_{\sigma \in \mathcal{F}_{h,k}^i} F_{\alpha,\sigma}^i(P^{n+1}, S_{\alpha}^{n+1}) + \sum_{\sigma \in \mathcal{F}_{h,k}^D} F_{\alpha,\sigma}^D(P^{n+1}, S_{\alpha}^{n+1}) \right) + \sum_{i \in \Pi | k \in \mathcal{K}_i^\Pi} F_{\alpha,k}^{\Pi,n+1}(P^{n+1}, S_{\alpha}^{n+1}) = 0. \quad (10)$$

582

The phase fluxes are approximated as follows:

583

For a face $\sigma \in \mathcal{F}_{h,k}^i$ such that $\sigma = \partial k \cap \partial l$:

$$F_{\alpha,\sigma}^i(P^{n+1}, S^{n+1}) = \lambda_{\alpha,\sigma}^{i,n+1} T_{\sigma}^i(P_k^{n+1} - P_l^{n+1}),$$

584

where $T_{\sigma}^i$ the face transmissivity is calculated using a harmonic average:

585

$$T_{\sigma}^{i} = |\sigma| \frac{K_{k}K_{l}}{K_{l}d_{k,\sigma} + K_{k}d_{l,\sigma}},$$

and

$$\lambda_{\alpha,\sigma}^{i,n+1} = \begin{cases} \lambda_{\alpha}(S_{\alpha,k}^{n+1}) & \text{if } T_{\sigma}^{i}(P_{k}^{n+1} - P_{l}^{n+1}) \ge 0, \\ \lambda_{\alpha}(S_{\alpha,l}^{n+1}) & \text{if } T_{\sigma}^{i}(P_{k}^{n+1} - P_{l}^{n+1}) < 0, \end{cases}$$

586

and with the mobility $\lambda_{\alpha}$ of the phase $\alpha \in \{w, g\}$ defined as:

$$\lambda_{\alpha}(S_{\alpha}) = \frac{kr_{\alpha}(S_{\alpha})}{\mu_{\alpha}},$$

587

for faces in $\mathcal{F}_{h,k}^b \cap \mathcal{F}_{h,\sigma}^D$

$$F_{\alpha,\sigma}^D(P^{n+1}, S^{n+1}) = \lambda_{\alpha,\sigma}^{D,n+1} T_{\sigma}^D(P_k^{n+1} - P_{\sigma}),$$

588

where $T_{\sigma}^D$ is given by:

$$T_{\sigma}^{D} = |\sigma| \frac{K_{k}}{d_{k,\sigma}},$$
589 and

$$ \lambda_{\alpha, \sigma}^{D,n+1} = \begin{cases} \lambda_{\alpha}(S_{\alpha,k}^{n+1}) & \text{if } T_{\sigma}^{D}(P_{k}^{n+1} - P^{\sigma}) \ge 0, \\ \lambda_{\alpha}(S_{b,\sigma}) = & \text{if } T_{\sigma}^{D}(P_{k}^{n+1} - P_{\sigma}) < 0. \end{cases} $$

590 For cells $k \in \mathcal{K}_i^{\Pi}$

$$ F_{\alpha,k}^{\Pi,n+1}(P^{n+1}, S^{n+1}) = \lambda_{\alpha,i}^{\Pi,n+1} IP_k(P_k^{n+1} - P_i^{\Pi}), $$

591 where the production index of the cell k $IP_k$ is given by Peaceman's formula Peaceman (1978):

$$ IP_k = 2\pi \frac{K_k}{\log(r_e/r_w)}, $$

592 with $r_e \approx 0.14\sqrt{dx^2 + dy^2}$, dx and dy are the respective perforated cell length and width of a cartesian mesh and

593 $r_w$ is the well radius (of 0.1 meter in this work). Finally:

$$ \lambda_{\alpha,i}^{\Pi,n+1} = \begin{cases} \lambda_{\alpha}(S_{\alpha,k}^{n+1}) & \text{if } IP_k(P_k^{n+1} - P_i^{\Pi}) \ge 0, \\ \lambda_{\alpha}(S_{\alpha,i}^{\Pi}) = & \text{if } IP_k(P_k^{n+1} - P_i^{\Pi}) < 0. \end{cases} $$

594 with $S_i^{\Pi}$ the gas saturation imposed at the well i.
List of Figures

<table><tr><td>1</td><td>Example of scenario and its impact on Newton's method. The scenario is composed of a well event opened after 2 time steps (red dashed line) and closed at 12 time steps (red dotted line).</td><td>25</td></tr><tr><td>2</td><td>Selected neural network architecture composed by an uplifting dense layer, four Fourier layers, a dense layer and a projection dense layer. Nc and Ni are parameters respectively representing the number of channels and the number of inputs. 95 and 60 correspond to the spatial discretization of the reservoir.</td><td>26</td></tr><tr><td>3</td><td>The Fourier layer starts with an input vector v, applies the Fourier transform F to it, then performs a linear transformation R on the lower Fourier modes while filtering out the higher modes. The inverse Fourier transform F<sup>-1</sup> is then applied. A local linear transformation W is applied to the original input vector v. The output of the top and bottom operations are then added together and an activation function is applied.</td><td>27</td></tr><tr><td>4</td><td>Qualitative view of Neural network input feature possible assembly.</td><td>28</td></tr><tr><td>5</td><td>Adapted 2D SHPCO2 case geometry</td><td>29</td></tr><tr><td>6</td><td>Test case 1 example of reservoir Pressure (left) and Saturation (right) obtained after a time-step with S0 = 3.81 × 10<sup>-4</sup>, q<sub>g</sub> = 7.61 × 10<sup>-4</sup> 1/s and dt = 2.4 × 10<sup>8</sup>s.</td><td>30</td></tr><tr><td>7</td><td>Test case 2 workflow with multiple well openings and closures. A step of time with a null well flow is realised between each closure and opening.</td><td>31</td></tr><tr><td>8</td><td>Example of reservoir Pressure (left) and Saturation (right) obtained after a time-step with q<sub>g</sub> = 9.8 × 10<sup>-4</sup>1/s and dt = 3.1 × 10<sup>8</sup>s. The simulation required 12 Newton iterations to converge.</td><td>32</td></tr><tr><td>9</td><td>L2 loss evolution through epochs for test case 1. The lowest test loss value is 1.9 × 10<sup>-3</sup> reached at epoch 17285. The corresponding train loss value is 2.0 × 10<sup>-3</sup>.</td><td>33</td></tr><tr><td>10</td><td>L2 loss evolution through epochs for the test case 2. The lowest test loss value is 8.7 × 10<sup>-4</sup> reached at epoch 7295. The corresponding train loss value is 9.3 × 10<sup>-4</sup>. We start with a learning rate of 10<sup>-4</sup>, at iteration number 1000, we change it to 5 × 10<sup>-5</sup>, then 10<sup>-5</sup> at iteration number 1600 and finally we set the learning rate to 5 × 10<sup>-7</sup> at iteration number 3600 and until the end.</td><td>34</td></tr><tr><td>11</td><td>Hybrid Newton approach applied to the example scenario and its impact on Newton's method.</td><td>35</td></tr><tr><td>12</td><td>Test case 1 scatter plot of the number of Newton iterations needed to converge using standard methodology versus using hybrid methodology on the train set (left) and on the test set (right). The color bar shows the distribution of Newton iterations using standard and Hybrid methodologies for the train and test set respectively.</td><td>36</td></tr><tr><td>13</td><td>Test case 2 scatter plot of the number of Newton iterations needed to converge using standard methodology versus using hybrid methodology on the train set (left) and on the test set (right). The color bar shows the distribution of Newton iterations using standard and Hybrid methodologies for the train and test set respectively.</td><td>37</td></tr></table>
Hybrid Newton Method for CO₂ Storage Simulation Acceleration

<figure><img src="image_1.png" /><figcaption>**Figure 1:** Example of scenario and its impact on Newton's method. The scenario is composed of a well event opened after 2 time steps (red dashed line) and closed at 12 time steps (red dotted line).</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>**Figure 2:** Selected neural network architecture composed by an uplifting dense layer, four Fourier layers, a dense layer and a projection dense layer. Nc and Ni are parameters respectively representing the number of channels and the number of inputs. 95 and 60 correspond to the spatial discretization of the reservoir.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 3:** The Fourier layer starts with an input vector **v**, applies the Fourier transform **F** to it, then performs a linear transformation **R** on the lower Fourier modes while filtering out the higher modes. The inverse Fourier transform **F**⁻¹ is then applied. A local linear transformation **W** is applied to the original input vector **v**. The output of the top and bottom operations are then added together and an activation function is applied.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 4:** Qualitative view of Neural network input feature possible assembly.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 5:** Adapted 2D SHPCO2 case geometry</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 6:** Test case 1 example of reservoir Pressure (left) and Saturation (right) obtained after a time-step with $S_0 = 3.81 \times 10^{-4}$, $q_g = 7.61 \times 10^{-4}$ 1/s and $dt = 2.4 \times 10^8$ s.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 7:** Test case 2 workflow with multiple well openings and closures. A step of time with a null well flow is realised between each closure and opening.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 8:** Example of reservoir Pressure (left) and Saturation (right) obtained after a time-step with $q_g = 9.8 \times 10^{-4}$1/s and $dt = 3.1 \times 10^8$ s. The simulation required 12 Newton iterations to converge.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 9:** L2 loss evolution through epochs for test case 1. The lowest test loss value is 1.9 × 10⁻³ reached at epoch 17285. The corresponding train loss value is 2.0 × 10⁻³.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 10:** L2 loss evolution through epochs for the test case 2. The lowest test loss value is 8.7 × 10⁻⁴ reached at epoch 7295. The corresponding train loss value is 9.3 × 10⁻⁴. We start with a learning rate of 10⁻⁴, at iteration number 1000, we change it to 5 × 10⁻⁵, then 10⁻⁵ at iteration number 1600 and finally we set the learning rate to 5 × 10⁻⁷ at iteration number 3600 and until the end.</figcaption></figure>
Hybrid Newton Method for CO₂ Storage Simulation Acceleration

<figure><img src="image_11.png" /><figcaption>**Figure 11:** Hybrid Newton approach applied to the example scenario and its impact on Newton's method.</figcaption></figure>
Hybrid Newton Method for CO₂ Storage Simulation Acceleration

<figure><img src="image_12.png" /><figcaption>**Figure 12:** Test case 1 scatter plot of the number of Newton iterations needed to converge using standard methodology versus using hybrid methodology on the train set (left) and on the test set (right). The color bar shows the distribution of Newton iterations using standard and Hybrid methodologies for the train and test set respectively.</figcaption></figure>
Hybrid Newton Method for CO₂ Storage Simulation Acceleration

<figure><img src="image_13.png" /><figcaption>**Figure 13:** Test case 2 scatter plot of the number of Newton iterations needed to converge using standard methodology versus using hybrid methodology on the train set (left) and on the test set (right). The color bar shows the distribution of Newton iterations using standard and Hybrid methodologies for the train and test set respectively.</figcaption></figure>