## 1 Introduction

The management of subsurface flow-related resources requires running computationally expensive simulations. This includes oil/gas production, Carbon Capture, Sequestration, and groundwater management. The cost gets even higher when there are optimization problems involved. For example, for the well placement optimization (WPO) problem, the goal is often to search the best well locations for each well to achieve the highest overall profit or net present value (NPV). Here, thousands of simulation runs need to be performed to find the optimal solution. For a practical optimization project, the whole process of running a simulation may take a few weeks or even months. With these concerns, researchers are currently studying the surrogate modeling to replace the full simulation so that those projects can be performed within a more reasonable time frame. One baseline model that has been previously explored is the surrogate using a low-fidelity simulation model with machine-learning error correction

[tang2021use]. This surrogate model achieves speedup for a factor of 10 to 100 with less than 1% difference in optimal NPV compared with running full simulation. This paper introduces a new surrogate model that can learn from a small set of realizations running from simulation and has the inductive capability to forecast the rest of the cases. This surrogate model is built based on Graph Neural Network (GNN) framework, so we term it as GSimGNN.)Unlike previous work [sanchez2020learning][sanchez2018graph] focusing on generating graph using a physical substance (particles or objects), GSimGNN directly encodes well locations and grids as nodes and build edges between their neighbors. Pfaff et al. [pfaff2020learning] had a similar approach. They build graphs using meshes, and they create the concept of world edges to connect meshes that are close in the distance. But the problem they solved was using the finite element method (FEM), which is different from ours (finite volume method). Our study shows that the GSimGNN can learn the grid-based multiphase flow simulations with acceptable error for well optimization-related applications.

## 2 Simulated Data Generation Method

This study generates the database for training, validation, and testing using a simulator because of the broad application of simulation in practice and industry, as we mentioned in the previous section. Currently, we generate a two-dimensional channelized reservoir model with oil and water. As shown in Fig. 1, we have the channelized permeability field, and the rest simulation parameters are described in Table. 1.

To begin with, we set up the model as a Buckley-Leverret problem that contains incompressible flow and rock and straight-line relative permeability relationship (mobility ratio is 1 for both water and oil phase). This results in the fixed pressure field map after the transient state of the simulation (usually reaches a steady state after the first time step). Therefore, the only changing state variable during the simulation is saturation. The Buckley-Leverret problem is a classical problem widely studied in subsurface flow theory, making it a solid starting point in this study.

No. | Property Name | Value & Unit |
---|---|---|

1. | Grid Number | |

2. | Grid Size | m m m |

3. | Porosity | |

4. | Initial Water Saturation | |

5. | Initial Pressure | bar |

6. | Injector BHP | bar |

7. | Producer BHP | bar |

8. | Total Time | days |

9. | Single Time Step | days |

10. | Time Step Number per Simulation | |

11. | Number of Injector | |

12. | Number of Producer |

One of the most critical applications of surrogate modeling in reservoir management is the well placement optimization (WPO), as it has significant impacts on the environment and economics. In the current study, we prepare the dataset by generating different well configurations for all wells ( injectors and producers as mentioned in Table. 1. This means that in each configuration, the placement of all wells is completely random and different from other configurations. By having this setup, our surrogate model can learn different scenarios from different well configurations and, therefore can be applied to WPO problems in practice.

## 3 GSimGNN Model Framework

### 3.1 Graph Generation

The simulation grids contain oil, water, and rock in a waterflooding reservoir simulation. In this study, the authors define a node as a grid shown in the Figure 2. It indicates the conversion from simulation grids to GSimGNN nodes. Due to the pressure gradient, the fluid (oil and water) flows between two neighbor grids. Thus the edge is built between them. The figure 3 shows three-by-three two dimensional grids with one injector and one producer setup.

In the initial state, The volume of oil and water and pressure field stay constant, but when pressure changes due to the kick of sink or source term, the volume will change with time. For the arbitrary grid i, the mass conserved in the following equation (Equation 1).

(1) |

where subscript indicates phases, superscript indicates components. is phase density, is the mass fraction of component in phase , is the Darcy velocity (given below), represents the source/sink term, is time, is porosity, and is the phase saturation. The Darcy velocity for phase , , is given by

(2) |

where

is the absolute permeability tensor,

, and denote relative permeability, viscosity and fluid pressure of phase , respectively. In this study, the system has two components and two phases: water and oil. To solve water and oil saturation for the arbitrary grid. The node needs to have the features shown in the Table 2. The author used one-hot encoding to represent different types of grids shown in the Figure

4.No. | Feature Name | Unit |
---|---|---|

1. | Pressure | Bar |

2. | Saturation | - |

3. | Permeability | Darcy |

### 3.2 Graph Neural Network Framework for Simulation

The current section explains the GNN workflow trained with the graph data generated from simulation data. The overall structure of the GNN workflow contains an encoder, a processor, and a decoder ([sanchez2020learning]), as is shown in Fig. 7. It takes the input of each node feature and outputs the updated state variable (saturation in this study) of every node for the next time step.

To begin with, an encoder network is performed on every single node in the graph. The object of this process is to embed the features and project them to a standard dimension. The feature embedding can then be used for further treatment. Next, a processor is applied to the embedded features. This process involves the interaction of different nodes in the graph, and the interactions are computed among nodes based on their connection relationship. The computation is achieved by doing message passing and aggregation between nodes. The state’s dynamic transition among nodes is then completed during the process. The processor outputs the updated embedded features for the next step. Finally, a decoder is used to generate the final output. Similar to the encoder, the decoder is applied to every single node in the graph. It projects the embedded output back to the desired dimension as the final output.

### 3.3 GSimGNN Implementation Details

This section discusses the implementation details of our GSimGNN surrogate model for reservoir simulation. Our implementation is done using Pytorch Geometric (PyG) library and referencing the previous related work

[sanchez2020learning].The input of the current study contains different features, i.e., current pressure, current saturation, node permeability, node well bottom hole pressure (BHP, default value for non-well cells), one-hot encoding to identify different types of cells (with 3 dimensions). The output currently only contains future saturation as the pressure field is steady during the Buckley-Leverett setup simulation. Both input and output are assigned to each node within the graph during the graph generation.

The encoder of this study is a multi-layer perceptron (MLP) of

linear layers with an activation layer after each. As mentioned in the previous section, the encoder is performed on every single node of the graph to a standard dimension of .Our processor is composed of single-step processor. Each of the single-step processors first collects information from neighbor nodes. The collected information is then stacked with current node features before being processed by another MLP of layers plus activation. The output of MLP is then aggregated with those from other neighbor nodes. The whole process is repeated times for each node so that each node can be updated with enough information from neighbors as well as their neighbors.

Like the encoder, our decoder is another MLP of linear ldayers with an activation layer after each. The decoder would be applied to each node to project them back to the output dimension which is in the current study (future saturation).

## 4 Results and Discussion

As mentioned in the previous section and Table. 1, we generated multiple simulation cases with different random well locations for 5 injectors and 5 producers. Each generated well configuration is simulated for 3000 total simulation days with 150-time steps. There are in total 100 well configurations generated in the current study which includes a total of 15000-time steps for GNN to learn. The dataset is split into 6000 samples for training (equivalent to 40 simulation cases), 3000 samples for validation and model tuning (equivalent to 20 simulation cases), and 6000 samples for testing (equivalent to 40 simulation cases).

We need to mention here that our GNN model is the original setup for predicting one time step during the training stage. An excellent one time step prediction can be achieved quite easily without much tuning during the testing stage. However, as the simulator performs multiple-step predictions, we also need to do rollout predictions with GNN (the next step prediction is based on the previous step prediction instead of given hard data). This could introduce accumulated errors that result in deviations from simulation results after multiple steps (in our case, more than 50 steps of rollout would usually introduce obvious errors). The problem can be resolved by using L1 loss (on top of the original L2 loss) plus Gaussian noise during the training stage.

Next, we demonstrate the GNN performance during the rollout with several figures. Here, we selected 5 out of 40 test case results for presentation purposes. We present both saturation maps and oil/water well rates. As is shown in Fig. 8, we plotted the saturation field map at 1500 days and 3000 days for simulation results, GNN results, and their difference. The difference between GNN results and simulation results at 1500 days ranges from 0.05 to -0.15 and 3000 days ranges from 0.1 to -0.2. As mentioned in the previous paragraph, the error accumulates during the rollout which explains the pattern of increasing error shown in the difference map. In visual, the saturation shape of both time steps generated from GNN matches well with the simulation results. This indicates that our GNN model learns to simulate multiphase flow in 2D. In Fig. 9, we pick three producers (with relatively higher rates) to show the ability of our GNN model to predict well rates which can be used for the well control and placement optimization. The red line (simulation results) and blackline (GNN results) align well with each other with a close trend near the end of the simulation. Similar patterns can be observed for the water rate as well. For other saturation maps (Fig. 10, 12, 14, 16) and well rates (Fig. 11, 13, 15, 17), we have similar observations.