Widely used in all phases of development of an oil field, reservoir simulations are applied to model flow dynamics and predict reservoir performance from preliminary exploration until the full development stage of petroleum production. A deep understanding of the model, discretization, and solution methods can be found in [20, 21].

A reservoir simulator is a combination of rock and geological structure properties, fluid properties, and a mathematical representation of the subsurface flow dynamics. A two-phase (oil and water) two-dimensional reservoir simulator is based on flow equations represented in Eq. 1. ∇ [ k → → k r , j ( S j ) μ j ρ j ∇ P j ] + q j − ∂ ∂ t ( ϕ S j ρ j ) = 0 , (1) where the subscript j represents the phase oil or water, k is the absolute permeability, k r is the relative permeability, µ is the viscosity, ρ is the fluid density, and ϕ is porosity. The solution we seek here is P j , S j .

It is possible to identify three terms on Eq. 1: the first one is the flux term, followed by the source/sink term, and the last one is the accumulation term. Although we only describe a two-phase flow system, these equations can easily generalize for multiphase flow simulations. As stated before, the reservoir simulation equations can be recast in a systems framework, where the nonlinear equations can be linearized into a time-varying state-space model as depicted in Figure 1. Note that while the controls u are imposed on the wells, the state x evolution is calculated based on the discretized partial differential equation (Eq. 1), and the output y is observed on the wells.

In general, the wells can be controlled by flow rates or bottom-hole pressure (BHP). Here we will assume that the producer wells are controlled by BHP and injector wells are controlled by injection rate ( q i n j ). So, we can define a control vector composed of the controls for all wells at a timestep t as: u t = [ BHP t q inj t ] , where BHP t = [ BHP 1 t … BHP p t ] T , q inj t = [ q inj , 1 t … q inj , i t ] T , p is the number of producers, and i the number of injectors.

The state dynamical evolution will be solved for each timestep t. The state is represented as: x t = [ P t S t ] , where P t = [ P 1 t … P n t ] T , S t = [ S 1 t … S n t ] T , and n is the number of gridblocks.

The output for the producers is the oil flow rate ( q o ) and the water flow rate ( q w ) and for the injectors is the bottom hole pressure (BHP). Thus, the output vector can be defined as: y t = [ q o t q w t BHP t ] , (2) where q o t = [ q o , 1 t … q o , p t ] T , q w t = [ q w , 1 t … q w , p t ] T , BHP t = [ BHP 1 t … BHP i t ] T . These values are calculated using the Peaceman equation [22]. For example, for a producer well, the flow rate of phase j can be calculated as: q j = k r , j B j μ j WI ( P j − B H P ) , (3) where WI is the well index and can be calculated as: WI = 2 π α k h ln ( r o r w ) + SKIN , where α is a unit conversion factor, h is the height of the reservoir gridblock, r o is the grid block dimension equivalent radius, r w is the wellbore radius, and SKIN is the skin factor.

Following the assumptions as in [3] and the application of fully implicit discretization it is possible to write the residual form of Eq. 1 as: g ( x t + 1 , u t + 1 ) = F ( x t + 1 ) + Acc ( x t + 1 , x t ) + Q ( x t + 1 , u t + 1 ) = 0 , (4) where we can identify the flux term F ( x t + 1 ) , the accumulation term Acc ( x t + 1 , x t ) , and the source and sink term Q ( x t + 1 , u t + 1 ) .

This equation can be solved by the Newton’s method by defining the Jacobian matrix J = ∂ g ∂ x . So, the state evolution can be calculated as: x t + 1 = x t − ( J t + 1 ) − 1 [ F ( x t + 1 ) + Acc ( x t + 1 , x t ) + Q ( x t + 1 , u t + 1 ) ] . (5)

To apply reduced-order modeling techniques to reservoir simulation efficiently, one can linearize the residual Eq. 4 (or Eq. 5). In Ref. [23], the authors presented the linearization for a single-phase 2D reservoir based on the simulator’s Jacobian matrix. Van Doren et al. [24] and Heijn et al. [25] developed similar approaches for a two-phase reservoir. These linearization methods were developed aligned to a control system approach, and the reader can find a comprehensive review of the control system applied to reservoir simulation in [3]. Committing with the control system approach, a state-space representation of a linear system can be written as Eq. 6. The authors cited in this paragraph manipulated equations similar to Eq. 4 or Eq. 5 to linearize them and write as a linear system. { x ˙ = A c x + B c u y = C c x + D c u , (6) where x represents the state, x ˙ the state time derivative ( d x d t ), and u the inputs (controls) and y the system output. A c is called state matrix, B c input matrix, C c output matrix, and D c is the feed-through matrix. The subscript c stands for continuous time. It is worth mentioning that on Eq. 6 we are presenting the state evolution equation (first one) and the output equation (second one), but the focus of linearization method is to be applied on the state evolution equation.

In general, it is necessary to choose a linearization point in time and use the Jacobian matrix at this time to calculate previously presented matrices. If the system dynamic strongly changes over time, the choice of a fixed linearization time can lead to large estimation errors. Techniques like TPWL (trajectory piecewise linearization) have been used to overcome these difficulties, where it is possible to approximate matrices A , B , C , and D at each interest time.

In the next section, we will present the linearization method TPWL that, combined with a model order reduction technique called proper orthogonal decomposition (POD), has been successfully applied to reservoir simulation [6].

Trajectory piecewise linearization is a method composed of two main steps. The first step is a training procedure (offline processing) that aims to capture state snapshots of a high-fidelity solution for the system to be linearized. These state snapshots are stored along with the Jacobian matrix of each timestep.

The purpose of the second step (online processing) is to predict the state on the next timestep x t + 1 based on the current state x t and the control state u t . The idea here is to use the stored states on the training step in order to predict the new state. First, in the group of the stored state snapshots, we identified the snapshot x i closest to x t . Then, we use the stored state and the Jacobian matrix for i and i + 1 to calculate the state x t + 1 using the following equation: x t + 1 = x i + 1 − ( J i + 1 ) − 1 [ F i + 1 + Ac c i + 1 + ∂ Ac c i + 1 ∂ x i ( x t − x i ) + Q ( x i + 1 , u t + 1 ) ] . (7) Previously, we wrote the state evolution equation of the reservoir simulation solution (Eq. 5) as a linearized system (Eq. 6). Using similar approach, we can rewrite the TPWL equation (Eq. 7) as a linearized system: { x t + 1 = A i x t + B i u t y t + 1 = C i x t + 1 + D i u t , (8) where A , B , C , and D are the time discretized version of the continuous matrices, and the superscript i represents the linearized matrices using information from the snapshot i.

To solve Eq. 7 and to compute matrices A i , B i , C i , and D i it is necessary to calculate the inverse of the Jacobian matrix , which is computationally costly. Furthermore, storing a full state and a full Jacobian matrix for each timestep on the training runs requires a huge memory amount. The next section will show how to use a model order reduced method to solve these problems.

The states stored on the training step of the TPWL can be arranged on: X = [ x 1 x 2 … x n ] . It is possible to take the singular value decomposition (SVD) of X . Based on some criteria, one can choose to keep only the l largest singular values. The POD basis Φ is defined as the l right singular vectors, and the reduced/latent state space z is then defined as a projection of the original space x based on the POD basis: z = Φ T x . (9) Using the relation x ≈ Φ z we can rewrite Eq. 7 as: J i + 1 Φ ( z t + 1 − z i + 1 ) = − [ F i + 1 + Ac c i + 1 + ∂ Ac c i + 1 ∂ x i Φ ( z t − z i ) + Q ( x i + 1 , u t + 1 ) ] , (10) and multiplying both sides of the above equation by Φ T : Φ T J i + 1 Φ ( z t + 1 − z i + 1 ) = − Φ T [ F i + 1 + Ac c i + 1 + ∂ Ac c i + 1 ∂ x i Φ ( z t − z i ) + Q ( x i + 1 , u t + 1 ) ] . (11) Now, defining the reduced Jacobian matrix as J r i + 1 = Φ T J i + 1 Φ and the other reduced terms as F r i + 1 = Φ T F i + 1 , Ac c r i + 1 = Φ T Ac c i + 1 , ( ∂ Ac c i + 1 / ∂ x i ) r = Φ T ( ∂ Ac c i + 1 / ∂ x i ) Φ , and Q r = Φ T Q , we can rewrite Eq. 11 as the POD-TPWL equation: z t + 1 = z i + 1 − ( J r i + 1 ) − 1 [ F r i + 1 + Ac c r i + 1 + ( ∂ Ac c i + 1 ∂ x i ) r ( z t − z i ) + Q r ( x i + 1 , u t + 1 ) ] . (12) By truncating the latent space, we deal with smaller dimension matrices when compared with the TPWL equation presented previously.

One can rewrite Eq. 12 as a linear time variant system: { z t + 1 = A r i z t + B r i u t y t + 1 = C r i z t + 1 + D r i u t , (13) where A r i , B r i , C r i , and D r i are the reduced version of the ones presented on Eq. 8.

The complexity of calculating matrices A r i , B r i , C r i , and D r i (or its equivalents on Eqs. 6 and 8) can increase in more realistic cases (black-oil or compositional simulators). Another drawback of these methods is the mandatory access to the reservoir simulator code to extract the Jacobian matrix and other data structures from the simulator, which may turn impractical the use of these methods with commercial reservoir simulators.

In [1], the authors develop a method to calculate matrices A r t and B r t based on deep learning, using state snapshots as training data. Ref. [2] extended this idea by adding a physical loss function related to reservoir simulation. Here we expand on this and demonstrate the feasibility of extracting not only states but outputs.

This method was proposed in [1] for model learning and control of a nonlinear dynamical system using raw pixel images as input. It uses a convolutional autoencoder coupled with a control linear system approach to predict the time evolution of the system state.

In Figure 2, we present a schematic diagram showing the main elements for training this model. The inputs used on the training procedure are the elements shaded in gray ( x t , x t + 1 , u t , and Δ t t ). The variables with a hat ( ⋅ ^ ) are the estimated values. The paths presented with continuous arrows are used in both training and prediction steps, and those with dashed arrows are only used on the training procedures.

The autoencoder is composed of an encoder that aims to transform the state in the original dimension space ( x t ) to a reduced (latent) space ( z t ), and a decoder that works as an inverted encoder to project back the reduced space state to the original space. One can relate the encoder with the POD projection (Eq. 9) since both have the ability to transform the original state x on the latent space state z.

The transition element uses the state in latent space to generate the matrices A r t and B r t . These matrices are used to handle the dynamical system evolution using a linear control system approach, as in Eq. 14. Since this equation operates on the latent space, we can relate A r t and B r t with their reduced version used on Eq. 13. z ^ t + 1 = A r t z t + B r t u t . (14) It should be pointed out that matrices A r t and B r t change for each different input state x t , like in other linearization procedures applied to reduced-order models techniques (e.g., TPWL—trajectory piecewise linearization).

Other important elements on this method are the 3 loss functions ( ℒ ) that will be used during the training procedure. The reconstruction loss function ℒ rec is introduced to guarantee the autoencoder (encoder and decoder) capability to reconstruct x ^ t as close as possible to x t . x ^ t is calculated using: x ^ t = Decoder ( Encoder ( x t ) ) . (15) The reconstruction loss function is defined as: ( ℒ rec ) i = { ‖ x t − x ^ t ‖ 2 2 } i , (16) where i represents the sample index.

Minimizing the prediction loss guarantees the accuracy of the time evolution of the dynamical system. In this case, we define the prediction loss function as: x ^ t + 1 = Decoder ( Transition ( Encoder ( x t ) ) ) , (17) where, x ^ t + 1 is the prediction of the state variables at t + 1 , and then is possible to calculate the prediction loss: ( ℒ pred ) i = { | | x t + 1 − x ^ t + 1 | | 2 2 } i . (18) The third loss function of the E2C model is the transition loss function, which tries to guarantee the model competence to evolve the reduced state z in time. This mimics the MOR idea of developing a projection framework in which the reduced state can evolve in time based on the reduced Jacobian matrix. First, we need to calculate: z ^ t + 1 = Transition ( Encoder ( x t ) ) , and then: z t + 1 = Encoder ( x t + 1 ) , and with these two equations it is possible to calculate the transition loss (Eq. 19): ( ℒ trans ) i = { ‖ z t + 1 − z ^ t + 1 ‖ 2 2 } i . (19) Here, we briefly described the main elements of the E2C model. In the next section, we will show an improvement of this model, adding a physical loss function specifically designed to reservoir simulation problems.

Inspired by the E2C model, Jin et al. [2] proposed an improvement on this model by adding a physical loss function to it. The model proposed by them focuses on the problem of building a proxy for a petroleum numerical reservoir simulator to be used on a well control optimization process. In Eq. 28 of [2], the authors introduced a physical loss function, reproduced below: ( ℒ p ) i = { ‖ k ⋅ [ ( ∇ p t − ∇ p ^ t ) recon + ( ∇ p t + 1 − ∇ p ^ t + 1 ) pred ] ‖ 2 2 } i + + γ { ‖ ( q w , t − q ^ w , t ) recon + ( q w , t + 1 − q ^ w , t + 1 ) recon ‖ } i . (20) This physical loss function is composed of two parts. The first one contains the mismatch between the predicted pressure gradient between adjacent grid blocks and the input pressure gradient multiplied by the permeability. This part of the proposed loss function intends to assure that the fluid flow between adjacent grid blocks in the predicted model is as close as possible to the input/true model. Simply put, this is basically a mass conservative requirement. Here, this first part of the physical loss function will be named the flux loss function.

The second part of the physical loss function deals with a mismatch in the producer wells’ flow rate. In fact, as described by the authors, it tries to guarantee that the producers’ grid block pressure in the predicted and reconstructed states are as close as possible to its equivalent on the input/true model. Instead of adding physics to the neural network, this second part adds more weight to the producers’ grid block pressure during the training procedure. In fact the authors in [2] did not use the flow rate on their physical loss function; instead, they used only the producers’ well grid block pressure. Since we will propose a new way to handle well data, we will omit this second part of the study’s proposed physical loss function. Even though we are omitting this part of the physical loss function, we will have it in our implementation of the method proposed by [2], for comparison purposes. Our implementation of their loss function (Eq. 28 of [2]) is presented here on Eq. 35.

We present in Figure 3 a diagram of the autoencoder with the flux loss function terms calculation. Comparing this figure with Figure 2, it is possible to notice the addition of the gradient calculation on both sides of the chart and also the calculation of the flux loss functions both in the reconstruction and prediction parts, which can be expressed as: ( ℒ flux , rec ) i = k { ‖ ∇ p t − ∇ p ^ t ‖ 2 2 } i , (21) and ( ℒ flux , pred ) i = k { ‖ ∇ p t + 1 − ∇ p ^ t + 1 ‖ 2 2 } i , (22) where k is the permeability.

The flux loss function can be defined as: ( ℒ flux ) i = ( ℒ flux , rec ) i + ( ℒ flux , pred ) i . (23) In addition to the flux loss function, on the top right of Figure 3, it is possible to see the calculation of the output y ^ t + 1 . The model outputs are well bottom hole pressure for the injector wells, and water and oil flow rates for the producer wells, in the case presented here. The calculation of these quantities from the state variables (pressures and saturations) is made by applying the Peaceman equation ([22]). This way to estimate the well data is here denoted by proposition 1. Note, that this is the actual proposed method in [2], implemented here for comparison purposes. In the following sections, we will propose and compare two other ways of estimating well data using E2C like models.

Since we are building a proxy model for a petroleum reservoir numerical simulator to be applied in a well control optimization framework, we would like to emphasize the well data prediction (BHP and flow rates), not only on the state (pressure and saturation) calculation; or, in other words, we want to reconstruct the output and not only the states. To achieve this, we propose a well data loss function named ℒ well _ data 2 as a reference to proposition 2. We present the diagram for proposition 2 in Figure 4, where it is possible to see the addition of the well data loss function on the right upper corner. This loss function will make the well data predicted by the model y ^ t + 1 as close as possible to the true value y t + 1 . Note that both y and y ^ are composed by oil and water flow rates and bottom hole pressures, as defined in Eq. 2. The input used to train the model is the well data calculated by a reservoir simulation. Proposition 2 well data loss function is defined as: ( ℒ well _ data 2 ) i = { ‖ y t + 1 − y ^ t + 1 ‖ 2 2 } i . (24) To get a more reliable state prediction and, as consequence, a better well data estimation, we introduce a two-phase flux loss function (Eq. 25) for a two-phase simulator. ( ℒ flux 2 Ph ) i = ( ℒ flux 2 Ph , rec ) i + ( ℒ flux 2 Ph , pred ) i , (25) where its terms are defined in Eqs 26, 27: ( ℒ flux 2 Ph , rec ) i = k { ‖ k r , o ( S w t ) ∇ p t − k r , o ( S ^ w t ) ∇ p ^ t ‖ 2 2 + ‖ k r , w ( S w t ) ∇ p t − k r , w ( S ^ w t ) ∇ p ^ t ‖ 2 2 } i , (26) and ( ℒ flux 2 Ph , pred ) i = k { ‖ k r , o ( S w t + 1 ) ∇ p t + 1 − k r , o ( S ^ w t + 1 ) ∇ p ^ t + 1 ‖ 2 2 + ‖ k r , w ( S w t + 1 ) ∇ p t + 1 − k r , w ( S ^ w t + 1 ) ∇ p ^ t + 1 ‖ 2 2 } i , (27) where k is the permeability, k r , j is the relative permeability to the phase j, and S w is the water saturation.

The addition of proposition 2 loss function ℒ well _ data 2 to the training process is a simple step but can improve the model’s ability to estimate well data. Although the model’s structure is very similar to proposition 1 (Figure 3), we expect to have more reliable well data estimation because we are adding this information to the training procedure. On proposition 1, the well data was not part of the training procedure. It was calculated on the validation/test steps based on the states.

Inspired by the E2C model, we extend its idea to the output data and named it as Embed to Control and Observe (E2CO). On the original E2C model, a transition network receiving as input the state in the reduced space was employed to generate matrix A r t and B r t , which was used to evolve the state in the latent space based on Eq. 14. We introduce another transition network, namely, the transition output. This network will receive information from the state on latent space z t to generate matrix C r t and D r t in order to be able to estimate well data y ^ t + 1 as in . y ^ t + 1 = C r t z ^ t + 1 + D r t u t . (28) The E2CO diagram is presented in Figure 5, where it is possible to see the new transition output network. We also create a loss functions ℒ well _ data 3 (as reference to proposition 3) that provides information to train the neural network transition output parameters. This loss function has the same structure as Eq. 24.

One of the main advantages of E2CO, when compared to propositions 1 and 2, is that here we do not rely on the state variables in the original space to predict well data. Instead of this, E2CO uses the state on the reduced space to create matrices C r t and D r t and uses a control system approach to estimate the output data. This resembles the case of system identification as in the control system community.

The diagrams presented in Figures 2–5 are used to train each proposed method’s networks. In these diagrams, the loss functions have an important role in the training procedure. They are defined to provide information about the system dynamics and honor the physics presented in the reservoir simulation process.

When the training procedure is finished, and the neural network parameters (weights, bias, and filters) are optimally found, it is time to use them to predict the system behavior. The prediction diagrams are significantly simpler than the training ones. The prediction diagram for propositions 1 and 2 is presented in Figure 6. In this figure, the time evolution happens in the latent space (z), potentially reducing the computational cost of the prediction step. However, to predict well data using propositions 1 and 2, it is necessary to have the state on the original space (x) at each desired time step to predict the well data at that time. That being said, one can note that it is necessary to pass z ^ t + 1 through the decoder to calculate x ^ t + 1 and further apply the Peaceman equation to calculate well data y ^ t + 1 on each time step. The application of the decoder to predict the well data can lead to a computational cost increase.

For proposition 3, there is no need to decode the latent space since the well data is generated on that latent space using a specific network for this purpose (Transition Output). Figure 7 shows that both the state time evolution and the well data estimation happen in latent space. On proposition 3, the decoder’s use is limited to the cases where one would like to estimate the state on the original space (x).

In the following sections, we describe the networks used in each part of the proposed models. The structure of the networks used here is very similar to the ones used in [2]. The three proposed ways of calculating well data use the same elements described in this section. Although the elements are the same, each proposed method’s parameters will differ due to the different model/network structure and different loss functions used in the training procedure.

We implemented all the methods proposed here using Python 3.7.4 [26] and the framework TensorFlow 2.3.1 [27] with Keras API [28]. The dimensions of layers, number and size of filters, and output size of each block on the encoder and decoder used are presented in Tables 1 and 2.

The transformation block has dimension of n Trans = 200 for the transition network, and the last dense layers before the matrices calculations have dimensions of l z 2 for A r t and l z ( n i + n p ) for B r t , where n i and n p are, respectively, the number of injector and producer wells. The output of these dense layers is then reshaped to match each matrix dimension. For the transition output network we are using n Trans _ WD = 20 for the transformation block dimension, and the last dense layers have dimensions of l z ( n i + 2 n p ) for C r t and ( n i + n p ) ( n i + 2 n p ) for D r t .

We used the Adam optimization algorithm with a learning rate of 1 × 10 − 4 during the training procedure with a batch size of 4. The sample is randomly selected to compose the training batch. We run the training procedure for 100 epochs.

We initially trained our models on CPUs. However, we change to GPUs since the speed up is expressive. We provide in Table 3 a non-rigorous estimation of the training time, where it is possible to see the differences between training on CPU and GPU. More details on the GPUs NVIDIA Tesla used can be found on www.hprc.tamu.edu.

We utilized a commercial numerical reservoir simulator (IMEX from [30]) to generate the data employed on the training and validation steps of our propositions. However, any commercial off-the-shelf reservoir simulator can be used to train our model. A Cartesian grid with 60 × 60 × 1 gridblocks and an oil–water fluid model were used. A fixed permeability field was generated and is presented in Figure 11 (left). The relative permeability curves used in the simulation model are shown in Figure 11 (right).

The model has five producers and four injector wells. The producer wells are controlled by bottom hole pressure (BHP), and the injectors are controlled by the water injection rate. This model runs through 2,000 days, changing controls every 100 days. Some properties used in the simulation model are presented in Table 4.

The objective of constructing a fast proxy model is to eventually connect it to control optimization frameworks whereby multiple calls of the surrogate model are used to determine a control set that optimizes an objective function (e.g., Net Present Value and sweep efficiency). In this case, we expect that our proxy might have the competence to predict well data under different sets of controls, which is usually unknown before the optimization is complete. To mimic the behavior of a varying input setting, we will use 300 simulations to train the proposed models, and each simulation will run through 20 time periods.

For the producers, the bottom-hole pressure was sampled over a uniform distribution U ( 260,275 ) kgf/cm². To avoid having the same volume of water injected in all control sets, we used a different approach to create the injectors controls. For each control set i, we sampled a base injection rate q i n j _ b a s e , i ∼ U ( 300,950 ) m³/day, and for each time period t, we sampled a perturbation q i n j _ p e r t , i , t ∼ U ( − 80,90 ) m³/day. The injection flow rate is then calculated as q i n j ( i , t ) = q i n j _ b a s e , i + q i n j _ p e r t , i t .

In Figure 12 we present an example of the injectors controls where readers can confirm the injection rate variability. The variability in the producer’s controls can be observed in Figure 13.

Using each one of the 300 control sets generated, we run the reservoir simulator and store the state (pressure and water saturation) at each one of the 20 time periods. The training set will be composed by 6,000 samples, each one containing:

• x t → state on timestep t

In the case presented here, t = 0,1 , … , 19 , where t = 0 represents the initial state that also needs to be stored as the first snapshot of the state ( x 0 ). The validation set contains 100 control sets that will create 2,000 samples. The validation set generation follows a similar approach but with a different random number generator seed. This is to enforce that training and validation are performed using different control inputs.

Each proposition (propositions 1 to 3, as constructed before) was trained as described in Section 2.7 using the training data set. We then used the models (networks) trained for each proposition on the prediction architecture, as shown in Figures 6, 7, and evaluated each model’s quality on the validation set. We also used this validation set to choose the better hyperparameters for each model. This is a tedious but necessary step in looking for a low prediction error model.