Saturation and Pressure Prediction for Multi-Layer Irregular Reservoirs with Variable Well Patterns

Abstract

The well pattern and boundary shape of reservoirs determine the distribution of the remaining oil distribution to a large extent, especially for small-scale reservoir blocks. However, it is difficult to replicate experiences from other reservoirs directly to predict the remaining oil distribution because of the variety of irregular boundary shapes and corresponding well patterns. Meanwhile, the regular well pattern can hardly suit irregular boundary shapes. In this paper, we propose a well placement method for undeveloped irregular reservoirs and a multi-step prediction framework to predict both oil saturation and pressure fields for any reservoir shape and well pattern. To boost the physical information of input characteristics, a feature amplification approach based on physical formulae is initially presented. Then, 3D convolution technology is employed for the first time in 3D reservoir prediction to increase the spatial information in the vertical direction of the reservoir in the input. Moreover, to complete the two-field prediction, the concept of multi-task learning is adopted for the first time, improving the rationality of the forecast. Through the loss-based ablation test, we found that the operation we adopt will increase the accuracy of prediction to some extent. By testing on both manually designed and real irregular-shape reservoirs, our method is proven to be an accurate and fast oil saturation prediction method with its prediction loss less than 0.01 and calculation time less than 10 s in the future one year.

1. Introduction

In the process of reservoir development, it is necessary to estimate the oil saturation distribution. According to the distribution of the remaining oil, the current development status of the reservoir can be known, and only after the oil saturation distribution field is obtained can the next development target be determined, which provides a basis for injection-production adjustment and new drilling in the later stage. However, for irregular reservoirs, the boundary configuration is very variable and well placement needs to be adjusted according to that, and different reservoirs have different permeability distributions and production regimes. Therefore, it is difficult to find a widely used residual oil prediction method for quantitative calculations in irregular reservoirs.

Methods for obtaining a saturation field are mainly divided into three kinds. One is to calculate using formulae. Leverett predicted the static saturation distribution for whole reservoirs [1] according to the oil saturation–height function. Another example is Archie, who makes use of the formation resistivity and porosity to calculate water saturation [2]. In recent years, the most commonly used method has been the flow unit [3] and EWU [4] model, which calculates the seepage rate and drainage efficiency between wells. Because their calculation formula is relatively simple, the calculation speed is very fast.

With the continuous enrichment of reservoir research content and the upgrading of research methods, researchers have put forward higher requirements for the accuracy and function of research methods. Therefore, the second class is based on numerical simulators, which can show the state of reservoirs to a large extent and is the most commonly used method at present [5]. At present, researchers have carried out research on the distribution of the remaining oil on both macroscopic and microscopic scales. Microscopically, by using digital rock physics technology [6] with the volume of fluid methods [7] and other methods, the exploration of the distribution law of remaining oil under different pore-throat structures is realized [8,9]. It helps technicians recognize the reason for the existence of remaining oil from the mechanism and find an effective method to reduce the remaining oil saturation under certain types of pore-throat structures [10]. Macroscopically, numerical simulation methods have largely simulated the flow and seepage phenomenon of conventional reservoirs [11]. There have been many detailed problems from different types of reservoir-solved simulation methods, such as water drainage management [12], water imbibition [13], corresponding effects [14], and even fracture effects of carbon reservoirs [15]. Furthermore, with the improvement of the understanding of the mechanism and the continuous upgrading of methods, more and more reservoir types can also be used [16]. However, the numerical simulator needs to make a geological model and fit history information [17,18], which is time-consuming and laborious. Meanwhile, irregular boundaries can distort the mesh shape and increase the difficulty of convergence.

The third class is based on neural networks, which are the most labor-saving, but not yet popularly used by people due to practice limitations. For example, according to previous articles [18,19,20], the prediction of remaining oil distribution by the neural network is mainly for single-layer integrated reservoirs and regular patterns of wells, and the working system is fixed. Therefore, it is difficult to carry out practical applications. In addition, another type of neural network method is based on the LSTM method [21], which requires thousands of numerical simulations to generate enough samples, resulting in a long sample generation process.

With the rapid development of computer vision technology, we have found that field data can also be regarded as a kind of high-dimensional picture. At present, relatively advanced work in the field of image generation mainly includes sample set generation [22], super-resolution [23], and style transfer [24], and a set of mature network frameworks has been formed to extract image features [25,26]. However, there are still some differences between the image and our saturation field. First, the image is two-dimensional data, while our field data is three-dimensional. Second, the image features only have three channels, but our field data have more.

Since we need to obtain the pressure field and saturation field at the same time, they are coupled in the calculation process. The concept of multi-task learning also needs to be employed for reference [27]. There are many branches of multi-task learning used in different domains, such as computer vision [28], natural language processing [29] and reinforcement learning [30]. Generally, tasks tend to accomplish prediction with the assistance of their own specific task indicators, which commonly are their instinct styles or arranged numbers. Unfortunately, our data do not have such clearly pre-defined indicators for the two fields.

In this paper, we propose a saturation field and pressure field prediction method for multilayer irregular reservoirs. This method can be used directly on new reservoirs with irregular boundary shapes after training, without the need to regenerate any sample based on the new reservoir for retraining. Technically, we first propose the most comprehensive input feature selection based on the numerical simulation formulae and propose some new physically meaningful feature combinations to make the prediction process more reasonable. Furthermore, for the first time, 3D convolution technology is employed to extract features in 3D field data from multiple scales in order to predict more correctly. Moreover, taking into account the advantages of implicit algorithms, we first make use of the shared trunk [31] and branch net [32] from multi-task network technology to complete the interaction between the oil saturation field and pressure field in the calculation process, which further improves the accuracy of prediction. Finally, we propose a new well placement method named the constrained-Kmeans well-placement method for irregular boundaries, leading to higher oil recovery than commonly used regular well patterns.

In the following sections, the theoretical background introduces the basic reservoir numerical simulation method and proposes the required feature types and feature combination patterns. Then, in the Section 3, we introduce the workflow and give the specific structure of the 3D convolutional module and shared trunk module in the network. In the Section 4, we describe the training process in detail by providing detailed training parameters and sample sources. In the Section 5, we analyze the application of our method in two cases and show the effectiveness of our method based on the results. Finally, the conclusions are summarized.

2. Theoretical Background

In the beginning, it is important to determine the rough mapping process to guide the construction of the network. Therefore, we resort to the formulae of both reservoir engineering and numerical simulation. Here, a simplification formula provides a clear goal for neural network mapping and a new idea for feature engineering.

For a two-phase grid reservoir model, we only pay attention to calculating the oil saturation and pressure of one phase named $\beta$ (in this paper, we use water). The oil saturation of another phase can be calculated based on this way.

According to the material balance equation and the continuity equation, we can obtain the following formula:

$$\frac{{\left({\rho }_{\beta }{S}_{\beta }\varphi V\right)}_i^{n+1}-{\left({\rho }_{\beta }{S}_{\beta }\varphi V\right)}_i^n}{\Delta t}={\displaystyle \sum }_{j\in {\eta }_i}{\left({\rho }_{\beta }{\lambda }_{\beta }\right)}_{ij+1/2}^{n+1}{\gamma }_{ij}\left({\varphi }_{\beta j}^{n+1}-{\varphi }_{\beta i}^{n+1}\right)+Q_{\beta i}^{n+1}$$

(1)

Further $\lambda$, $\gamma$, and $\Phi$ are expressed as follows:

$${\lambda }_{\beta }=\frac{K{k}_{r\beta }}{{\mu }_{\beta }}$$

(2)

$${\gamma }_{ij}=\frac{A_{ij}{k}_{ij+1/2}}{d_i+d_j}$$

(3)

$$Q_{\beta i}^{n+1}=\frac{\rho K{k}_{r\beta }A}{{\mu }_{\beta }}{\displaystyle \sum }_{j\in {\eta }_i}\left({\varphi }_{\beta j}^{n+1}-{\varphi }_{\beta i}^{n+1}\right)+W_i+I_i$$

(4)

$${\varphi }_{\beta i}^{n+1}=P_{\beta i}^{n+1}-{\left({\rho }_{\beta }\right)}_{ij+1/2}^{n+1}g{D}_i$$

(5)

For each grid of formation, ${{\displaystyle \sum }}_{j\in {\eta }_i}\left({\varphi }_{\beta j}^{n+1}-{\varphi }_{\beta i}^{n+1}\right)$ means the sum of the difference in the potential of grids between grid $i$ and its neighboring grids and the $\varphi$ could be expressed by $P$ according to Equation (5).

Then, ${{\displaystyle \sum }}_{j\in {\eta }_i}\left({\varphi }_{\beta j}^{n+1}-{\varphi }_{\beta i}^{n+1}\right)$ in Equation (1) can be transformed into $\Delta {{\displaystyle \sum }}_{j\in \left(x,y\right)}\frac{\left(P_{\beta j+1}^{n+1}-P_{\beta j}^{n+1}\right)-\left(P_{\beta j}^{n+1}-P_{\beta j-1}^{n+1}\right)}{2\ast \Delta }$, where $x$, $y$ is the direction of coordinates. Then, the right side of the plus sign can be regarded as a Laplace aggregation and can further turn into $\Delta L{P}^{n+1}$ when grids are in the same scale in both $x$ and $y$ directions.

Finally, the right side of the equal in Equation (7) becomes:

$$\frac{\rho K{k}_{r\beta }}{{\mu }_{\beta }}{\displaystyle \sum }_{j\in {\eta }_i}\frac{A_{ij}{k}_{ij+1/2}}{d_i+d_j}\left({\varphi }_{\beta j}^{n+1}-{\varphi }_{\beta i}^{n+1}\right)+\frac{\rho K{k}_{r\beta }A}{{\mu }_{\beta }}{\displaystyle \sum }_{j\in {\eta }_i}\left({\varphi }_{\beta j}^{n+1}-{\varphi }_{\beta i}^{n+1}\right)+W_i+I_i$$

(6)

and then further becomes:

$$\frac{\rho K{k}_r}{\mu }\Delta \ast L^l\ast \left(1+C_{\rho }gD\right)P^{n+1}-W_i-I_i$$

(7)

where $l$ means power $l$ of $L$, denotes the scale of Laplace operation in the cube, and will correlate to the length of time step.

Thus, Equation (1) becomes:

$${\left(\rho S\varphi V\right)}^{n+1}+\Delta t\Delta \ast L^l\ast \left(1+C_{\rho }gD\right)P^{n+1}\frac{\rho K{k}_r}{\mu }={\left(\rho S\varphi V\right)}^n-\Delta t{W}_i-\Delta t{I}_i$$

(8)

Here, we focus on water phase saturation and water relative permeability and thus omit the subscript symbol for simplicity.

In addition, we can find $k_r$ in Equation (6) can be expressed by the Corey function and varies with $S$ as shown in Equation (9).

$$k_{rw}={(k_{rw})}_{S_{orw}}{[\frac{S_w-S_{wc}}{1-S_{wc}-S_{orw}}]}^{C_w}$$

(9)

To simplify our equation, we assume that the liquid is incompressible. Then, Equation (7) will transform to the following form after merging Equations (8) and (9).

For each grid:

$$\begin{array}{l}\left({\varphi }_0-{\varphi }_0{C}_{\varphi }\right) \mu \rho V S_w^{n+1}+{\varphi }_0{C}_{\varphi }\mu \rho V{P}^{n+1}{S}_w^{n+1}+t\Delta L_l\left(1+gD\right)P^{n+1}\rho K {(k_{rw})}_{S_{orw}}{[\frac{S_w^{n+1}-S_{wc}}{1-S_{wc}-S_{orw}}]}^{C_w}\\ =\left({\varphi }_0-{\varphi }_0{C}_{\varphi }\right) \mu \rho V S_w^n+{\varphi }_0{C}_{\varphi }\mu \rho V{P}^n{S}_w^n-t\mu W^{n+1}+t\mu I^{n+1}\end{array}$$

(10)

Here we assume the exponential value $C_w$ in Corey function is 4 and $S_w^n-S_{WC}=\widetilde{S_w^n}$, Then Equation (10) becomes:

$$\begin{array}{l}\left({\varphi }_0-{\varphi }_0{C}_{\varphi }\right)\mu \rho V \widetilde{S_w^{n+1}}+\left({\varphi }_0-{\varphi }_0{C}_{\varphi }\right)\mu \rho V{S}_{wc}+{\varphi }_0{C}_{\varphi }\mu \rho V{P}^{n+1}\widetilde{S_w^{n+1}}+{\varphi }_0{C}_{\varphi }\mu \rho V{P}^{n+1}{S}_{wc}^{}\\ +t\Delta L^l\left(1+gD\right)\rho K {(k_{rw})}_{S_{orw}}\frac{1}{{\left(1-S_{wc}-S_{orw}\right)}^4}{P}^{n+1}\widetilde{{ \mathrm{S}}_w^{n+1}{}^4}\\ =\left({\varphi }_0-{\varphi }_0{C}_{\varphi }\right)\mu \rho V\widetilde{S_w^n}+\left({\varphi }_0-{\varphi }_0{C}_{\varphi }\right)\mu \rho V{S}_{wc}+{\varphi }_0{C}_{\varphi }\mu \rho V{P}^n\widetilde{S_w^n}+{\varphi }_0{C}_{\varphi }\mu \rho V{P}^n{S}_{wc}^{}-t\mathsf{\mu }{W}^{n+1}-t\mu I^{n+1}\end{array}$$

(11)

When expressing constant coefficients of every variable with single symbols, Equation (11) becomes:

$$a\widetilde{ s_w^{n+1}}+b+c{P}^{n+1}\widetilde{ s_w^{n+1}}+d{P}^{n+1}+eK{P}^{n+1}{\widetilde{{ \mathrm{S}}_w^{n+1}}}^4=f\widetilde{S_w^n}+g+h{P}^n\widetilde{ s_w^n}+i{P}^n+t\mu W^{n+1}-t\mu I^{n+1}$$

(12)

Thus:

$$P^{n+1}\sim \left(W, I, K, D, P^n,\widetilde{S^n},P^n\widetilde{S^n},\dots ,\widetilde{S^{n+1}}\right)$$

(13)

and:

$$S^{n+1}\sim \left(W, I,K,D, P^n,\widetilde{S^n},P^n\widetilde{S^n},\dots ,P^{n+1}\right)$$

(14)

It is worth noticing that, for compressible fluid, the form of $P^n$ is going to come in terms of squares and cubes.

Combining Equations (11)–(14), we can draw three conclusions. First, if we want to get the oil saturation $S^{n+1}$ and $P^{n+1}$ at a later time step, seven variables and their combinations from the previous time step are needed, including the permeability field $K$, saturation at the previous step $S^n$ and at the current step $S^{n+1}$, pressure $P^n$ at the previous step and $P^{n+1}$ at the current step, the production and injection volume of the well $W_i$ and $I_i$, and some further non-linear combinations. Second, the $l$ power of $L$ is difficult to determine. When time goes by, the range of network diffusion gradually increases, leading to the number of $l$ being non-stationary. Third, in these implicit equations, saturation and pressure integrate together in a very complex way, making it difficult to calculate by analytical methods.

3. Methodology

In response to these problems, our framework contains several operations to deal with them. First, in the feature engineering part, we collect all required features according to the formulae and combine them in a non-linear way by means of feature amplification in the feature engineering process. Since the various features in the physical calculation process are expressed in the form of weights, we adopt 0–1 normalization to match the effect of the weight. Then, we employ a self-learned method in real time for the hidden state of fields through a self-defined encoder network that contains multi-layer convolution of field data and collects multiple-scale features by itself, such as dealing with uncertain power $l$ of $L$. Finally, we use a decoder network to achieve the decoupling of part of features in two target fields with a distillation operation to accomplish interaction with each other during output.

The whole workflow of our method is shown in Figure 1. We will describe some key parts of our method in the following subsections.

3.1. Feature Engineering

There are two things going on in this part of the work: feature group amplification and feature normalization. First, we increase the types of features group like they are in Equation (14), such as $P^n$ ∗ $\widetilde{S^n}$, By doing so, we release the difficulties of non-linear fitting and reduce the risk of overfitting. In total, we finish the feature group of {D, $K$, $W^{n+1}$, $I^{n+1}$,$P^n$, $\widetilde{S^n}$, $P^n\widetilde{S^n}$} and use them as our input features.

After that, we cast the data normalization. Here, we find that it is better to keep these features acting like the weight because that makes the size of the value can express the importance more clearly. Therefore, we adopt the interval [0,1] as the final target range. The final canonical equation is Equation (15)

$${\tilde{x}}^i=\frac{x^i-x^i{}_{min}}{x^i{}_{max}-x^i{}_{min}}$$

(15)

where $i$ is the dimension index of data, $\tilde{x}$ is the normalized data, $x$ is the raw data, $x_{min}$ is the minimum value of $x$, and $x_{max}$ is the maximum value of $x$.

3.2. Network Structure

The neural network structure is mainly divided into two parts: the encoder part is used to convert input data into hidden features of the reservoir, and the decoder part is used to decode hidden features into observed states.

3.2.1. Encoder Network

The function of this part is to extract the features in many scales and provide information for the next saturation and permeability fields, as shown in Figure 2.

This structure is designed and refers to the encoder part of the VGG network [33]. But our field data has 3 dimensions in shape; thus, we upgrade the convolution layers from 2D to 3D [34], as shown in Extractor modules in Figure 3. Precisely, it extracts features from three directions, such as a cube. To be specific, this part consists of multiple Extractor modules, each of which is constructed by one pooling layer and two 3D-convolution layers followed by the Leaky Relu function [35] as an activation function. Meanwhile, two convolutional layers share different convolution kernel sizes for multi-scale abstraction. Precisely, the larger-size convolution kernel increases the number of features, while the smaller size convolution kernel improves the mapping ability. Plus, the middle pooling layer reduces the dimension and results in stacking features in multi-scales. After multiple extractors, the input data will be fully compressed and amplified. In practice, the field data is compressed to a scale of 3 × 3 × 3, and its number of channels increases to 256.

3.2.2. Decoder Network

The decoder network is constructed to generate a saturation field and pressure field according to the hidden features calculated by the encoder part above. The entire network contains 3 functional modules, namely the generator module, the model distillation module [36] and the GRU module [37], as shown in Figure 4. To be specific, the generator module restores the size to the original size. Then, the distillation module accomplishes the interaction between the oil saturation field and pressure field in Equation (7). Finally, the GRU module combines shared features and independent features into complete intrinsic features for each field.

In detail, the generator module consists of three parts and its structure is shown in Figure 5. The first part is up-sampling processing, which is to expand the size of the three-dimensional field. The second and third parts are deconvolution processing, which is the reverse of the convolution operation mentioned earlier and the size of kernels changes corresponding to the encoder. Another module is named model distillation, which merges features from two fields and generates the interactive features back for them. This module can be constructed in many ways and here we perform a naive concatenation of the feature followed by two 3D convolution layers with kernel size 1 by 1 by 1. The third module GRU contains many formulae (Equations (16)–(19)) or fuses the distilled interactive features into the original branch so that they will include both mutual effects and their unique properties. Finally, the two-layer convolution works for a better fit. In addition, the activation functions of the last two layers of networks are set to Tanh() [38] and Sigmoid() [39].

$$r_t=\sigma \left(W_{}^r\ast H_t+U_{}^r\ast H{L}_{t-1}\right)$$

(16)

Here, the Reset Gate determines the degree of change in hidden features, where $\sigma$ is an active function and matrixes $U^r$ and $W^r$ are weight parameters, and $H_t$ and $H{L}_{t-1}$ are features from the distillation unit and specific field generator.

$$i_t=\sigma \left(W\ast H_t+r_t\odot U\ast H{L}_{t-1}\right)$$

(17)

$i_t$ is all the hidden features of the new input control data, including self-specific features.

$$z_t=\sigma \left(W_{}^z{H}_t+U_{}^{z}H{L}_{t-1}\right)$$

(18)

$z_t$ is the Update Gate, which determines the change degree of origin features affected by fusion features, and it is also a weight parameter.

$$H{L}_t=z_t\odot H{L}_{t-1}+\left(1-z_t\right)\odot i_t$$

(19)

Then, $H{L}_t$ is the final updated specific features of each field.

3.3. Performance Metrics

We used the absolute error ($L_1$) as our loss function, which measures the absolute error between the target results and prediction values, as shown in Equation (20). We also found that L2 leads to blurred predictions. In addition, it should be noted that the sampling method of the weight parameters has a great influence on the result, and the experimental results show that it is easier to jump out of the range of the local optimal solution to realize the initialization of the weight by using Gaussian distribution.

$$L_1={\displaystyle \sum }_{i=1}^n\left|y_i-{\widehat{y}}_i\right|$$

(20)

where n is the total number of samples, $y_i$ and ${\widehat{y}}_i$ are the target results and the corresponding prediction value in the grid block, respectively.

3.4. Constrained-Kmeans Well-Placement Method

In response to the difficulty of deploying regular well patterns in irregular reservoirs, we propose a method named constrained-Kmeans for the deployment of irregular reservoir well patterns. We take advantage of the Kmeans [40] method to cluster spatial grid locations and deploy new Wells at central points. The shape of the model varies greatly, so the location of the new well is completely different. At the same time, considering the benefits problem, the number of wells needs to be constrained based on the total geological reserves, which are calculated by Equation (21). Therefore, single-well control reserves are promised. The deployment effect of this method is shown in Figure 6.

$$N=\frac{\varphi \ast \rho \ast \left(1-S_w\right)\ast X\ast Y\ast Z}{B_{oi}}$$

(21)

where X, Y, and Z represent the scale of reservoirs in the x, y and z directions, respectively. B_oi represents the oil volume factor of the initial underground state.

Yellow areas are active grids in the formation, purple areas are inactive, and blue scatters are positions of wells.

4. Complete Workflow

• Step 1. Data collection

In training mode, sample generation and sample pool building work are carried out with the help of the Gaussian process to generate N permeability fields and random shapes selected from the triangle, quadrilateral and pentagonal as boundaries. Then, we define random types of wells at random locations (active grids) and define random production strategies (injection and production rate) for these wells. We did it randomly because we hope to obtain the distribution of the entire sample. Then we run a simulator to obtain both the saturation field and pressure field over time steps. Next, we split the time series into individual time steps. Because the permeability field and pressure field are included, the whole sequence can be regarded as a Markov process, so the data for each time step is an independent sample. To be precise, that is, for each production time step, the oil saturation field and pressure field before production are used as the input data of the sample, and the data of the two fields after production are used as the predicted labels.

When it comes to testing or applications, target information collection is needed. Precisely, for an irregularly targeted reservoir, we first deploy a uniform well pattern according to the previously proposed constrained-Kmeans method to achieve even and complete control of reserves. After designing the willingness regime for all wells, we obtain basically all raw features, including the initial permeability field, initial pressure field, well location and well injection/production rate. Then, we go to step 2. The most initial saturation field and pressure field show a uniform value for each field plane.

• Step 2. Feature normalization and amplification

According to the feature normalization and amplification methods mentioned above, the numerical range of features has shrunk to between 0 and 1. At the same time, different types of features will be combined in different ways to increase the types of features in the input.

• Step 3. Network calculation

In the training model, 4000 iterations are used to complete the training. The validation set is used as the basis for model parameter selection to improve the generalization ability of the model.

When in the application, continuous prediction of saturation and pressure field is as follows: When an initial pressure field and saturation field are given, we recursively use our method to let it continuously predict both saturation and pressure field changes with time steps to obtain the remaining oil distribution and pressure changes at different times.

• Loss calculation and backpropagation

After obtaining the oil saturation field and pressure field predicted by the network, we apply the L1 loss function mentioned above to calculate the error between the predicted result and the labeled result. Specifically, in the training mode, we apply the Adam algorithm [41] to calculate the gradient information of parameters in the network and further update them given an assigned learning rate.

• Result analysis

We show the loss of both training data and validating data in the training mode and further analyze the contribution of inner modules embedded in the network based on the loss curve [42,43,44]. When it comes to application, we do a prediction of saturation and pressure field in future 9 timesteps from a random reservoir state and show the trend of the error between the predicted field and real ones.

5. Case Study

In this section, our network is trained with samples generated from reservoir blocks with random boundary shapes and corresponding irregular well patterns. The geological model used as a training set is built with a parameter sample from

Table 1. Geostatistical parameters used to generate the geological models.

Parameter	Value	Parameter	Value
dimensions	20 × 20 × 8	porosity	0.25
initial water saturation	0.17	residual oil saturation	0.23
oil viscosity	5 cp	water viscosity	1 cp
well injection pressure upper limit	600 bar	well production pressure bottom limit	100 bar
water injection rate	1–3000 m3/d	well production rate	1–3000 m3/d

. Then, we test it on three geological models [45,46] containing two manually designed models and one real model. Our goal is to recurrently predict saturation and pressure fields in the future 12 timesteps, where each step is a month.

5.1. Training Data Collection

This part describes the source of data and the way of data processing, which provides data support for the later network operation.

5.1.1. Sample Generation

We construct permeability fields by sampling from 200 permeability fields generated by the open-source package SGeMS. We assume that the discretized permeability field in a heterogeneous reservoir follows a log-normal distribution. The mean value of log-permeability is 7, and the variance is 0.5. Detailed parameters of the geological model are shown in Table 1. Then, we take a random one of a triangle, square and pentagon as the boundary shape for each reservoir, as shown in Figure 7. We found that most of the boundary shapes of unconventional reservoirs can be regarded as one of these three shapes. Vertically, we use boundary shape and property distribution as the same as each layer.

After that, we sample random positions to deploy wells and run them within the control of random production regimes. After collecting the results of the simulation and splitting them by timesteps, we obtain input data. To be precise, each sample includes the permeability field, the oil saturation field and pressure field at the current time, the position of wells, and the subsequent working regimes, including injection and production rate. What we are expecting is to output the oil saturation field and pressure field at the next moment.

5.1.2. Sample Pool Building

When running simulators 1000 times, we obtain a total of 1000 sequences, each of which contains 10 timesteps. Moreover, the interval of each timestep is 30 days. After collecting the results of the simulation and splitting them by timesteps, we obtain a total of 10,000 sample data points. To select the best model, we divide our data into a training dataset and a validating dataset by a ratio of 0.9 to 0.1. Precisely, the training set is used to update the parameters of the network according to the gradient calculated by the backpropagation process, while the validating set uses a trained network to calculate the loss on them and select the real better parameters. This operation overcomes the overfitting problem to some extent.

It can be seen from this point that our sample-generating method can obtain more samples than that of the general recursive network method because the recurrent network needs multiple complete sequences by the simulator to be samples, but our method splits each sequence into n samples where n equals the number of time steps [47].

5.2. Feature Engineering

First, we do the normalization process by scaling original input features into the range of [0,1]. Then, we transform sparse data, such as well features, into dense form (field data) to represent well-local information until all well-based information is turned into field type. For example, for the well position feature, if a grid has a well, it shows different values from others that do not have one. Then, the features are composed into many groups as amplified features. Detailed results of feature engineering are shown in Figure 8.

5.3. Network Training

This part mainly describes the training methods and training results. Through the training results, we can more clearly see the strength of our method.

5.3.1. Training Strategy

We split raw data according to the Cross Validation method [48]. Precisely, the data are split into many mini-batches, and 90% of data in each min-batch is used for training, while the rest is validating datasets to select the best parameters.

In backpropagation, we used a learning rate of 0.0005 and reduced it to 0.95 times per 100 epoch steps. In addition, we designed the final number of hidden layer features to be 512. It is worth mentioning that we find that a higher learning rate leads to worse convergence of the model to the local optimal solution. We also found that when the number of hidden layer features increases or decreases, the computing accuracy of the network will decline.

5.3.2. Training Loss

According to the loss after 3000 iterative steps, the sum loss of saturation and pressure for prediction in training and validating gradually converges to only 0.003 and 0.013, respectively, as shown in Figure 9.

5.4. Ablation Test

In this section, we show how every module is important when making predictions in our method.

By comparing the final training loss (sum of saturation loss and pressure loss) of the complete model and six incomplete models. We obtain a ranking of their losses, as shown in

Table 2. The training loss of the complete model and 5 incomplete model.

Model	2D Convolution	No Distillation	[−1 to 1] Normalization	Less Hidden Features	Less Hidden Layers	No Feature Amplification	Complete
sum loss (training)	0.0042	0.0034	0.1401	0.0150	0.0036	0.0032	0.0028
sum loss (validation)	0.0087	0.0225	0.1419	0.0198	0.0165	0.0150	0.0069

and Figure 10. Here, the most obvious difference is between the results from two different normalization methods. When input data are normalized into the range from −1 to 1, the network easily converges to a very unreasonable regional optimal point that is far from the real global optimal point. Besides, fewer hidden units severely ruin our prediction accuracy both in training and testing. This is caused by the insufficient expressing ability of hidden layers when equipped with not enough neurons.

Apart from the two modules mentioned above, the rest of the networks seem to work the same in the training stage. However, they are very different when focusing on the validating stage. First of all, when a network lacks distillation between two fields, the network may find it hard to concatenate more information from another branch in the Decoder and thus result in the omission of some important information. The less important factors are the number of hidden layers and feature amplification work, and they work with the same importance as each other. They work almost the same because they both work to provide information for the network. But the difference between them is hidden layers to ensure that important information in the forward propagation is transmitted, and the feature amplification module is to ensure that the input information is sufficient to predict the result. The closest to the complete model is the model using 2D convolution. The 2D convolution method is less useful than the 3D convolution method for reservoir data, as shown in the red line in Figure 10. The reason is that there is a lack of feature extraction in the vertical direction.

5.5. K-Means Placement Method

To show the advantage of this deployment method, we compare the remaining oil distribution of the reservoir within the boundary of the bottom right picture of Figure 11 after the same developing policy under our well placement and commonly used five-point placement. The results are shown in Figure 11 and Figure 12 where the blue points and black points in Figure 12 represent injection well and production well respectively. We can find that, under the same liquid production rate, the accumulative production of oil by our method is more than the five-point method and our method makes it easier to recover the remaining oil around the boundary.

5.6. Prediction Result

We tested three examples using trained network models, the first two of which are artificial models with randomly generated boundary and permeability fields, to test the generalization performance of our method. The third example is an actual small-scale reservoir used to test the practicality of our method.

5.6.1. Case 1: Two Manually Designed Reservoirs

We regenerate triangles and pentagon-shaped reservoirs with the same grid size and grid dimensions as the training dataset and test them from an undeveloped state for future 12 timesteps where each timestep is one month. By using the constrained-Kmeans method, the well placement of the reservoirs is shown in Figure 13. Their production and injection quantities are sampled randomly from a uniform distribution of 1–3000, while limited by a constant bottom hole pressure value.

We make use of 0–1 normalization and feature amplification to increase the number of meaningful features, as described above. Then, the features are imported into the network to predict both saturation and pressure at each step recurrently. In practice, except for the first timestep when we need to input the reservoir state field as a restart state, at each subsequent timestep, the network goes through based on the result of the previous step.

The losses calculated by the four algorithms are shown in

Table 3. Error values from the four types of measurements for the last predicted steps.

Loss Type	Case 1		Case 2
	Saturation	Pressure	Saturation	Pressure
MAE(L1)	0.0094	0.0095	0.0062	0.0123
MSE(L2)	0.0060	0.0072	0.0047	0.0083
R2	0.9820	0.9423	0.9925	0.9755
adjusted R2	0.9763	0.9393	0.9893	0.9701

. The error given by L1 loss is the highest because L1 pays more attention to the singular values with a large deviation between the predicted results and the real results, while L2 pays more attention to the approximation of the mean. In addition, R2 is similar to adjusted R2 and is still relatively close to 1, indicating that our prediction has a high degree of accuracy. Additionally, we can clearly see that for both case 1 and case 2, the predictive loss for pressure is always higher than that of the oil saturation field. This is because the data for the saturation field is relatively more steady compared to the pressure field, with only a large change in gradient at the front edge of displacement, while most of the region shows approximate numerical values, especially in high water saturation and non-water saturation areas. In comparison, the distribution of the pressure field is much more complex, not only showing the pressure difference between injection and production wells but also reflecting the superposition of the potential function generated by multiple injection and production wells at any position. Furthermore, case 2 shows a greater loss than case 1, which is due to the fact that each layer in case 2 is not completely horizontal, as in case 1. Therefore, in addition to predicting the impact of injection and production, it is also necessary to identify the gravitational effect caused by the difference in depth within the same layer, increasing the difficulty of prediction.

From the L1 loss of every step for both fields, as shown in Figure 14, we can find the trend that the error increases with the steps. As a result, when we make predictions, we need to trade the length of the prediction and the prediction accuracy. In addition, when the water drainage process steps into the final stage, the error becomes stable for the prediction of the oil saturation field. Figure 15 and Figure 16 show how the reservoir saturation and pressure field change from the initial state throughout the 12 timesteps. In brief, we exhibit layers 1, 3, 5 and 7 to represent the properties of the top, above-middle, bellow-middle and bottom layers of the reservoir. Literally, our network accurately predicts the oil saturation distribution in both horizontal and vertical directions. Meanwhile, the total time-consuming and final loss of our method shown in

Table 4. The time-consuming nature of both the network and simulator.

Method	Case 1	Case 2
Our method	1.134 s	1.128 s
simulator	15 s	13 s

illustrates that our method saves 90% time for prediction compared with the numerical simulator. As our method uses inner matrix linear calculations combined with activation functions to perform calculations, there are no iterations or convergence issues involved in the entire calculation process. Therefore, all calculations can be obtained in an analytical way, so the computing speed is much faster than numerical simulators.

5.6.2. Case 2: A Real Reservoir

To further verify the practicability of our method, we test it on an actual reservoir with irregular lithological boundaries. This model is an independent turbidite sand body with medium porosity, medium–low permeability percolation and anomalously high pressure. The entire reservoir is buried at a depth of 3200 m with a thickness of 9–16 m. The specific parameters are shown in

Table 5. Parameters of the real reservoir block.

Items	Value	Item	Value
area (km2)	0.84	porosity (%)	19.4
reserves (104 t)	32	permeability (mD)	3.1
thickness (m)	4.5	viscosity (mPa.s)	23.1
Depth (m)	3200	density (g/cm3)	0.8793
formation pressure (MPa)	44.1	freezing point (°C)	37
pressure coefficient	1.38	salinity (mg/L)	26,661
Temperature (°C)	123	formation water	CaCl2

. We used a grid size of 50 m × 50 m × 2 m to coarsen the original model to eventually generate 20 × 20 × 8 grids to match our network. The historical development was divided into two development stages: the elastic development stage and the water injection development stage. Specifically, the rapid decline in formation energy during elastic development resulted in low production from some of the producing wells. Then, water injection was started at well X3 to replenish the formation energy, as shown in Figure 17. However, due to the limited injection capacity, the energy supplement was not sufficient to maintain high production rates. In this case, considering that well X6 has started to produce water and the lack of formation energy is severe, the developer plans to turn well X6 into a water injection well in the future for more energy replenishment. We desire to forecast the feature states of reservoirs when design developing policy is used to make sure that it is profitable. The data in this article are calculated using Schlumberger’s numerical simulator. Specifically, the model uses the static properties of the actual reservoir, including permeability, saturation and net-to-gross ratio. Then, fluid modeling and relative permeability modeling are established using on-site PVT data and relative permeability data. Finally, the actual development regime is set as a control parameter to run the simulator and the saturation field and pressure field during the historical stage are obtained after a process of historical fitting work. In addition, we run more than 10 time-steps to the future as labels of upcoming predictions to test the accuracy of our method.

As a result, we input the production and injection regimes to be used in the coming year (12 steps and 1 month each) into our program and use the calculation results of the numerical simulator as label data. In order to use the available data as much as possible, we finetune our network with historical data. To prevent overfitting, we train them only in less than 30 epochs before forecasting the future. From the comparison shown in Figure 18, the prediction results start to deviate from labels gradually with the steps increasing the same as the manually designed geo-model. Meanwhile, as the water drive area in the oil saturation field is gradually fixed, the oil saturation variation in individual time steps becomes smaller and smaller, and the cumulative effect of errors diminishes so that the losses gradually smooth out. Furthermore, from Figure 19, it can be seen that the prediction error mainly originates from the saturation difference at the displacement front. Therefore, as the number of grids covered by the displacement front increases, the difficulty of predicting saturation increases, leading to an increase in loss. Gradually, the drainage area tends to be constant, leading to the position of the displacement front becoming fixed, and the variation in water saturation in the reservoir also stabilizes. Therefore, the prediction loss starts to decrease, and the cumulative loss becomes smoother, as shown in Figure 18. Additionally, in Case 2, the type of injection-production wells remains unchanged in the whole developing process, and their injection-production rates vary little, reflecting a very steady pressure field with little difference in the pressure field between each step. This explains why the prediction results for the pressure field in Case 2 are better than those of the oil saturation field.

Moreover, the prediction speed of our method is still much faster than that of the numerical simulator, as shown in Figure 20, and even faster for the larger geological model. On the other hand, the history fitting process can be done automatically by our method, which further saves the time needed for prediction preparation.

6. Summary and Conclusions

This paper proposes a neural network method for predicting the oil saturation field and permeability field for a 3D irregular reservoir. The method embeds physical formulas, adopts 3D convolution technology, and integrates the concept of multi-task learning to achieve accurate predictions of saturation and pressure fields. The implementation of cases shows that this method can achieve good prediction results within a certain time scale. Some detailed conclusions are shown below.

• A sequential saturation and pressure prediction method for a reservoir with an irregular boundary shape and an irregular well pattern is proposed. This method is faster than the numerical simulation method and can provide an accurate saturation and pressure field. After training, the method can predict the oil saturation and pressure field of the irregular boundary-shaped reservoir several times under any reasonable production regime conditions according to the state of the initial reservoir.
• A number of modules and technical processes in the framework have improved prediction performance to some extent. The constrained-Kmeans well placement method is helpful for uniform well layout in the irregular reservoir because it achieves a 3% increase in oil recovery compared with regular well patterns such as five-point methods.
• The modules we employed in our framework have their unique contributions. Data normalization plays the most important role in performance, as it determines the success of subsequent operations. Feature selection and feature combination according to theoretical formulas can provide more information for network calculation, thus reducing the loss to a lower degree. Compared to 2D convolution technology, 3D convolution technology is more favorable for the extraction and transformation of vertical multi-scale spatial features. As a result, it reduces to a fourth of the L1 loss value. In addition, the use of the multi-task concept by the distilling module further reduces loss by around 10% of loss value.
• This method is proved accurate in the prediction of the future 10–12 months and particularly within 5 months. The L1 loss of prediction is only 0.0025 and it only increases to 0.01 at timestep 12 in the real reservoir geo-model. Meanwhile, its time consumption is one-twelfth that of the simulator.

In future work, it is possible to use the Physics-informed Neural Network method to replace the neurons in our network with partial differential terms in differential equations to increase the physical information in the calculation process and better guide the prediction process.