Deep Neural Networks Predicting Oil Movement in a Development Unit

TL;DR

This paper introduces a deep learning-based workflow that significantly accelerates 3D reservoir simulation data analysis, enabling efficient forecasting of oil flow, pressure, and saturation dynamics for improved oilfield development planning.

Contribution

The novel Metamodel combines Variational Autoencoders and Recurrent Neural Networks to efficiently predict reservoir dynamics without compromising accuracy, vastly speeding up traditional simulation methods.

Findings

The Metamodel accurately reconstructs reservoir dynamics.

It forecasts flow rates, pressure, and saturation effectively.

The workflow is orders of magnitude faster than conventional methods.

Abstract

We present a novel technique for assessing the dynamics of multiphase fluid flow in the oil reservoir. We demonstrate an efficient workflow for handling the 3D reservoir simulation data in a way which is orders of magnitude faster than the conventional routine. The workflow (we call it "Metamodel") is based on a projection of the system dynamics into a latent variable space, using Variational Autoencoder model, where Recurrent Neural Network predicts the dynamics. We show that being trained on multiple results of the conventional reservoir modelling, the Metamodel does not compromise the accuracy of the reservoir dynamics reconstruction in a significant way. It allows forecasting not only the flow rates from the wells, but also the dynamics of pressure and fluid saturations within the reservoir. The results open a new perspective in the optimization of oilfield development as the scenario screening could be accelerated sufficiently.

Figures

Tables

Table 1. Table 1: Prediction error between forecasted and ground truth reservoir states f t subscript 𝑓 𝑡 f_{t} for different metamodels. All the values are averaged throughout the validation dataset.

Mean relative error for $f_t$ in % and its standard deviation
Model	$d_z$	Error $\pm$ std, %			RMSE
		$p$	$s_{oil}$	$s_{wat}$	$p$ (bar)	$s_{oil}$	$s_{wat}$
Linear+PCA	$150$	$16.5\pm 10.8$	$59.5\pm 15.4$	$38.5\pm 7.7$	32.88	0.26	0.27
Linear+VAE	$150$	$28.7\pm 20.2$	$35.4\pm 8.5$	$22.4\pm 8.8$	63.42	0.16	0.16
Linear+VAE	$250$	$11.1\pm 8.7$	$22.7\pm 12.4$	$12.8\pm 6.5$	21.42	0.10	0.09
GRU+PCA	$150$	$43.9\pm 30.1$	$67.3\pm 5.6$	$42.4\pm 9.4$	90.33	0.30	0.30
GRU+VAE	$150$	$8.0\pm 5.1$	$11.4\pm 6.2$	$6.7\pm 3.2$	19.13	0.05	0.05
GRU+VAE	$250$	$9.5\pm 6.5$	$14.7\pm 8.5$	$8.6\pm 4.2$	21.37	0.06	0.06
LSTM+VAE	$150$	$8.9\pm 6.9$	$12.8\pm 8.7$	$7.5\pm 4.6$	18.72	0.06	0.05

Equations

$$\mathbf{E}_{m,u}\Big{[}\frac{1}{T}\sum_{t=0}^{T}||\hat{F}(\hat{f}_{0:t},u_{0:t},m)-F(f_{0:t},u_{0:t},m)||^{2}_{2}\Big{]}\rightarrow\min_{\hat{F}},$$

$$\mathbf{E}_{m,u}\Big{[}\frac{1}{T}\sum_{t=0}^{T}||\hat{R}(\hat{f}_{0:t},u_{0:t},m)-R(f_{0:t},u_{0:t},m)||^{2}_{2}\Big{]}\rightarrow\min_{\hat{R}},$$

$$f_0=F^{init}(m),{\hat{f}}_0={\hat{F}}^{init}(m).$$

$$z_{t+1}=F^{latent}(z_{0:t},u_{0:t},m).$$

$$\mathbf{E}_{f_{0:T}\sim\mathbf{D}}\Big{[}\sum_{t=0}^{T}\delta(f_{t},G^{-1}_{\theta}(G_{\theta}(f_{t}))\Big{]}\rightarrow\min_{\theta}.$$

$$X=U\Sigma V^T.$$

$$G_{\theta }(f_t)={\mathbf{E}}_{z_t\sim q_{\theta }(z_t|f_t)}[z_t],$$

$$G_{\theta }^{-1}(z_t)={\mathbf{E}}_{f_t\sim p_{\theta }(f_t|z_t)}[f_t],$$

$${\hat{z}}_{t+1}=W_z{\hat{z}}_t+W_u{u}_t+W_{m}m+b,$$

$$a^l=\sigma (W^l{a}^{l-1}+b^l),$$

$$a^0=concatenation({\hat{z}}_t,u_t,m),$$

$${\hat{z}}_{t+1}=a^L,$$

$$h_{t+1}={\sigma }_r(W_r{h}_t+W_z{\hat{z}}_t+W_u{u}_t+W_{m}m+b),$$

$${\hat{z}}_{t+1}={\sigma }_o(W_o{h}_{t+1}+b_o),$$

$${\hat{z}}_0=G_{\theta }(F^{init}(m)),$$

$${\hat{z}}_0\approx W_{i}m+b_i.$$

$$V_t^{\alpha }(p_{ref})=V^{bulk}\varphi (p_{ref})\sum _{i=1}^{n_x{n}_y{n}_z}{s}_{ti}^{\alpha }{e}^{(c_r-c_{\alpha })(p_{ti}-p_{ref})},$$

$$r_t^{\alpha }(p_{ref})=\frac{V_t^{\alpha }(p_{ref})-V_{t-1}^{\alpha }(p_{ref})}{\Delta t}+q_t^{\alpha }(p_{ref}),$$

$$r_t^{\alpha }(p_{surface})=\frac{r_t^{\alpha }(p_{ref})}{b_{\alpha }(p_{ref})},$$

$$r_t^g(p_{surface})=\frac{r_t^g(p_{ref})}{b_g(p_{ref})}+\gamma r_t^o(p_{ref}),$$

Full text

a Skolkovo Institute of Science and Technology,

Nobelya Ulitsa 3, Moscow, Russia 121205

b Gazpromneft Science and Technology Centre,

75-79 liter D Moika River emb., Saint Petersburg, Russia 190000

K****eywords Metamodel  $\cdot$ Reservoir Model  $\cdot$ Reduced Order Model  $\cdot$ Neural Network  $\cdot$ Machine Learning

1 List of Abbreviations

[TABLE]

2 Introduction

Computer 3D modelling of oil flows through a porous medium is the most frequently used tool for an oil field development optimization and prediction of unknown reservoir properties by history-matching procedure history-matching-1990 . It is essential to have an accurate and fast enough model of oil flows to solve both problems.

The common modelling approach is to use a finite-difference hydrodynamical simulation (FDHS) ertekin-bars-2001 . FDHS uses computations on a discrete computational grid, and its accuracy is strongly dependent on a resolution of the grid. Being accurate, computations on fine grids may consume weeks to be performed for a big oil field. It is almost impossible to run iterative optimisation and history-matching methods using this kind of models in sufficient time.

There are many techniques aimed to speed up the simulation process: starting from simple downscaling and local grid refinement and ending with data-driven Reduced Order Modelling (ROM) approaches. Capacitance-Resistance Models (CRM) connectivitycrm-2005 ; kovalcrm-2015 ; fullycoupledcrm-2014 , for example, tries to fit nonlinear autoregressive model able to forecast fluid production rates. While CRM retains the physicality of the process, other models train autoregressive production model in a fully data-driven manner insim-2015 ; rpm-2016 ; shm-2layer-2015 . On the other hand, Galerkin projection methods are aimed to find the linearly reduced approximation of Partial Differential Equations (PDE) solved by FDHS drrnn-2018 ; mor-in-fd ; srom ; pwl-mor . The overall idea of the approach is to construct projection bases using Proper Orthogonal Decomposition (POD). For nonlinear PDEs straight Galerkin projection is computationally inefficient due to the need of computing fully-dimensional Jacobian matrix at each timestep. Some approaches tries to overcome this issue via additional approximation layer, e.g. Empirical Interpolation eim and Discrete Empirical Interpolation deim methods.

Petroleum society requires a fast and accurate model for three-phase flow simulation on a 3-dimensional computational grid. And all the discussed approaches cannot solve the task properly: either they are not able to produce solutions for spatial fields (pore pressure, saturation, etc.) or they can work only under nonrealistic greatly simplified settings.

In this work, we provided a methodology for a fast data-driven 3D reservoir modelling, based on machine learning techniques and aimed to approximate FDHS solutions on a subset of possible simulation options. We use “Metamodel” terminology GTApprox2016 to describe a purely data-driven 3D reservoir model approximating the solution given by a base model (e.g. FDHS), and we provide a detailed definition in Section 3.

We concentrated our attention on modelling so-called development units - symmetric subelements of a whole oilfield. We discussed scenarios with water injection and without injection at all. The approach can be straightforwardly used to train a metamodel for different injection fluids though this theme is out of the scope of this paper.

The idea of forecasting for just a development unit can be used to construct a metamodel of a whole oilfield via division it into several independent units. On the other hand, fully-convolutional neural networks can represent a way to scale the proposed approach to a whole oilfield. Experiments with whole and real cases are retained for future work.

Proposed metamodelling approach is based on the idea of forecasting in latent variable space. The same idea was used in e2c-2015 ; rce-2017 to approximate the dynamics of an environment with high-dimensional state space for model-based Reinforcement Learning. Authors propose to use Variational Autoencoders vae-2013 for the model order reduction and we try it in our framework.

Our major contributions are:

• •

The metamodelling framework based on forecasting in latent variable space.

• •

Metamodel of a field development unit composed of Convolutional Variational Autoencoder for dimensionality reduction and GRU Recurrent Neural Network gru-2014 for the latent dynamics prediction. The metamodel is able to predict both pressure-saturation distributions (Section 5) and production rates from wells (Section 6).

• •

Comparative analysis of different metamodels in the sense of FDHS approximation accuracy, see Section 8.

3 Metamodelling of Reservoir Flows

We mentioned that our approach is developed to work in a subset of possible simulation options. The overall operational envelope of the presented model has two major restrictions:

Geometry of the computational domain is restricted to a field development unit, 2. 2.

Reservoirs are homogeneous in horizontal dimensions.

It may be noticed that if wells are allocated on a regular grid and the reservoir is homogeneous in horizontal dimensions, then the dynamics will be the same for all elements of the grid. In petroleum engineering, such an element is usually referred to as a development unit. We investigate three common well allocation schemes (fig. 1):

• •

Staggered line drive (fig. 1a),

• •

Direct line drive (fig. 1b),

• •

Without injection (fig. 1c).

A model that forecasts dynamics for a development unit (actually for $1/4$-th of it) may be directly applied for symmetric reservoirs. Though real reservoirs frequently are neither strictly symmetric nor horizontally homogeneous, they may be discretised into development units of required form. The ensembling of units produced by such a discretisation was left for the future work. Throughout the rest of the paper, we consider forecasting dynamics and production rates for 1/4-th of a development unit with one production well and, possibly, one injection well.

We will think of a reservoir dynamics as a non-Markovian Decision Process (nMDP) with a time discretization step $\Delta t$. The nMPD is defined using 5-tuple $(\mathcal{F},\mathcal{U},\mathcal{M},F,R)$, where:

• •

$\mathcal{F}$ is a set of reservoir states $f_{t}\in\mathcal{F}$ representing distributions of pore volume pressure $p_{t}$ and oil, water and gas saturations $s_{t}^{o},s_{t}^{w},s_{t}^{g}$ within 3D coordinate grid at a fixed moment of time $t$. Since $s_{t}^{o}(x,y,z)+s_{t}^{w}(x,y,z)+s_{t}^{g}(x,y,z)=1$ we may drop one of the saturations from the reservoir state and compute it explicitly when needed (e.g. gas saturation).

• •

$\mathcal{U}$ is a set of control variables $u_{t}\in\mathcal{U}$ representing production bottomhole pressure and water injection rate (if injection is presented).

• •

$\mathcal{M}$ is a set of variables $m\in\mathcal{M}$ defining reservoir initial conditions, rock and fluids properties. This variable is invariant of time and is fixed for a concrete development unit. We refer to these variables as a metadata of a development unit.

• •

$F$ is a transition function that maps applied control sequence $u_{0:t}$, all previous states $f_{0:t}$ and metadata of development unit $m$ into a following reservoir state $f_{t+1}$, such that $f_{t+1}=F(f_{0:t},u_{0:t},m)$

• •

$R$ is a production rates function, such that $r_{t}=R(f_{0:t},u_{0:t},m)$, where $r_{t}$ is a vector of oil, water and gas production rates at time $t$. Production rates may be thought as a kind of analogy to reward in classical MPDs.

There are models able to restore the transition function $F$, and production rates function $R$ by solving corresponding partial differential equations (PDEs). We call them Base Models. Commonly, Base Models solve PDEs using finite-difference or finite-volume schemes on a discretised grid and differ from each other in used PDE’s modifications and numerical integration schemes. Further, we will concentrate on so-called finite-difference hydrodynamical simulators (FDHS), which solve three-phase fluid flow equations with gas liberation/dissolving from/into the oil options. Though, our approach can be straightforwardly generalised for other types of Base Models.

Thinking of FDHS as of an ideal flow model (in the sense of accuracy) we may want to speed up it by substituting it with a metamodel. To proceed, we need a formal definition of a reservoir metamodel (see also GTApprox2016 ):

Definition 1

Reservoir metamodel is a data-driven reservoir model aimed to approximate the transition function $F$ and production rates function $R$ represented by a Base Model.

Here and further we will assume the approximation quality to be measured in the sense of L2-norm. If we denote approximated transition and production functions by $\hat{F}$ and $\hat{R}$ respectively then the approximation problem may be formulated as follows:

[TABLE]

[TABLE]

where the sequence of approximated reservoir states $\hat{f}_{0:t}$ is obtained from the approximated transition function $\hat{F}$, and the initial states are computed as follows:

[TABLE]

In the following sections, we will discuss the training procedure of a data-driven reservoir metamodel starting from a description of the training set.

4 Training Data

In order to construct a metamodel we use supervised machine learning techniques, requiring a training set of examples of the Base Model dynamics. Training set $\mathbf{D}$ is composed of scenarios (episodes), where each scenario represents one run of FDHS and contains:

Input variables:

• •

Scenario metadata in compressed vector form: $m\in\mathbf{R}^{61}$ (for more details see A). Metadata is fixed throughout the scenario.

• •

Sequence of control variables $u_{0:T}$, $u_{t}\in\mathbf{R}^{2}$ consists of production bottomhole pressure and water injection rate at time $t$. 2. 2.

Output variables (FDHS solutions):

• •

Sequence of reservoir states $f_{0:T}$. Each reservoir state $f_{t}\in\mathbf{R}^{3\times n_{x}\times n_{y}\times n_{z}}$ is a tensor containing values of pore pressure, oil and water saturation for all cells of computational grid at time $t$, where $n_{x}$, $n_{y}$, $n_{z}$ are numbers of cells in corresponding dimensions.

• •

Sequence of oil, water and gas production rates $r_{0:T}$, $r_{t}\in\mathbf{R}^{3}$.

The output variables (2) were obtained from multiple FDHS runs on randomly generated input variables (1). The generative model of input variables is based on the analysis of real laboratory data and aimed to produce realistic samples (see B for details).

Each scenario has its own time horizon $T$. FDHS turns off if the oil production rate falls below the limit of 3 stb/day, with an upper bound: $T<30$ years.

5 Latent Variable Reservoir Dynamics

Direct application of basic ML techniques to the approximation of functions $F$ and $R$ is complicated by the so-called curse of dimensionality. Even for simplified Markovian version of the task (when $f_{t+1}=F(f_{t},u_{t},m)$), the basic linear regression model will require learning of more then $(3n_{x}n_{y}n_{z})^{2}$ parameters, which will lead to the overfitting of the metamodel for any realistic number of available training examples — the metamodel will not be able to reproduce training accuracy on scenarios not included in the training set. Moreover, such a model will not be efficient both in terms of computational speed and memory consumption.

One way to reduce the problem complexity is to find a bijective mapping (dimensionality reduction) $G:\mathcal{F}\rightarrow\mathcal{Z}$, where $\mathcal{Z}$ is a low-dimensional latent variable space, and to define the dynamics in this latent variable space:

[TABLE]

The graphical representation of this latent variable reservoir dynamics model is shown in Figure 2.

Full-order reservoir states $f_{t}$ might be recovered as $f_{t}=G^{-1}(z_{t})$. It is not required to recover $f_{t}$ at each timestep – all the dynamics can be computed in the latent variable space with the need to use inverse mapping $G^{-1}$ only for the timesteps of interest.

In the following subsections, we show some practical models for both dimensionality reduction function $G$ and latent dynamics function $F^{latent}$. Different combinations of these models will result in different reservoir metamodels.

5.1 Search for a mapping $G$ for the dimensionality reduction

From an ML point of view, dimensionality reduction problem may be solved by fitting the so-called autoencoding model. We will use autoencoders to represent the true bijective mapping $G$ and its inverse $G^{-1}$ with parametrized approximations $G_{\theta}$ (encoder) and $G^{-1}_{\theta}$ (decoder).

Autoencoder passes the input $f_{t}$ through a low-dimensional bottleneck via subsequent application of encoding and decoding functions, and is trained to minimise some deviation $\delta$ between autoencoder’s input and output:

[TABLE]

Throughout this work, we analysed two distinct autoencoding models:

• •

Principal Component Analysis,

• •

Convolutional Variational Autoencoder.

Principal Components Analysis (PCA) is a linear dimensionality reduction model and is based on Singular Value Decomposition (SVD) of the matrix $X$ composed of all tensors $f_{t}\in\mathbf{D}$ flattened into vectors:

[TABLE]

$d_{z}$ most significant singular values are retained to form truncated matrix of right singular vectors $\tilde{V}\in\mathbf{R}^{(3n_{x}n_{y}n_{z})\times d_{z}}$. $G_{\theta}(f_{t})=\tilde{V}^{T}f_{t}$ is now can be used as an encoder and $G_{\theta}^{-1}(z_{t})=\tilde{V}z_{t}$ – as a decoder.

PCA provides the best approximation accuracy throughout all linear models with fixed latent space dimensionality $d_{z}$ in the sense of Frobenius norm $\delta(f,f^{\prime})=||f-f^{\prime}||^{2}_{2}$.

However, PCA has two significant disadvantages: it is linear, and it ignores the spatial structure of the tensor $f_{t}$ due to the flattening procedure and fully-connected structure.

Convolutional Variational Autoencoder combines the application of Variational Autoencoding framework (VAE) vae-2013 with convolutional neural network modules.

Variational Autoencoder is a probabilistic model that tries to approximate not just a point estimate of $z_{t}$, but rather the posterior distribution over it with probability density function $q_{\theta}(z_{t}|f_{t})$. The task is solved by maximization of Evidence Lower Bound Objective (ELBO) proposed in the original paper.

The parametric form of the conditional distribution $q_{\theta}(z_{t}|f_{t})$ can be chosen to be Gaussian, and in this case, the ELBO maximization is closely related to the minimization of L2-norm of deviation in (5), so we stick with the Gaussian approach.

Due to the deterministic structure of latent dynamics models described below, we must eliminate the stochasticity of VAE. We use the following expressions for encoder and decoder respectively:

[TABLE]

[TABLE]

where $p_{\theta}(f_{t}|z_{t})$ is a likelihood function.

Convolutional neural networks colah-conv-2015 ; alexnet-2012 ; vgg-2014 are used to represent both likelihood function $p_{\theta}(f_{t}|z_{t})$ and posterior distribution $q_{\theta}(z_{t}|f_{t})$. The use of convolutional layers allows us to overcome both linearity and indifference to the spatial structure of tensors $f_{t}$. We use 3D convolutional layers with elementwise nonlinearity between them (Leaky Rectified Linear Unit – LReLU leakyrelu-2015 ). Dimensionality reduction comes from the stride of the convolution window.

Convolutions exploit spatial invariance property that is also presented in FDHS engine: time derivatives of a property in a computational cell are dependent only on the parameters in the neighboring cells. Spatial invariance allows using less number of trainable parameters $\theta$ which helps us to overcome overfitting problem.

Graphical representation of convolutional VAE is provided at Figure 3.

5.2 Approximation of the latent dynamics

Given a set of latent variables $z_{t}=G_{\theta}(f_{t})\;\forall f_{t}\in\mathbf{D}$ we can now approximate the dynamics $F^{latent}$ with $\hat{F}^{latent}$ using standard ML techniques.

Most of ML models assume the input variable of a fixed dimensionality, which restricts the direct applicability of these models to the approximation of non-Markovian dynamics $F^{latent}$ since its inputs $z_{0:t}$ and $u_{0:t}$ are of variable length. For this kind of models we may allow the approximation to be Markovian instead: $\hat{z}_{t+1}=\hat{F}^{latent}(\hat{z}_{t},u_{t},m)$.

Different possible choices of the latent dynamics model are provided below. We denote the set of model’s parameters by $\phi$ and will fit them to minimize the average L2 loss in the latent space.

Linear Regression (LR) model requires Markovian simplification and follows the equation:

[TABLE]

where parameters $\phi=\{W_{z},W_{u},W_{m},b\}$ are obtained analytically.

Fully-Connected Neural Network (FCNN) is another Markovian model able to approximate a much broader class of functions than LR. FCNN is composed of $L$ consecutive fully-connected layers with outputs $a^{l},\;\forall l\in[1,\;\dots,L]$:

[TABLE]

where $\sigma$ is an elementwise nonlinear function. Parameters $\phi=\bigcup_{l=1}^{L}\{W^{l},b^{l}\}$ might be trained via one of gradient-based optimization procedures.

Recurrent Neural Network with Gated Recurrent Units (GRU RNN) gru-2014 compared to FCNN and LR models is a non-Markovian one. GRU RNN is a modification of a simple Recurrent Neural Network, and its internal structure is out of the scope of this work. For the sake of brevity, we provide here the equations for a simple RNN model. For $\forall t\in[0,\;\dots,T-1]$:

[TABLE]

where $h$ is a hidden layer vector, $\sigma_{o}$ and $\sigma_{r}$ are nonlinear functions; parameters $\phi=\{W_{o},W_{r},W_{z},W_{u},W_{m},b_{o},b\}$ are obtained from gradient-based optimization procedure.

We also tried to compare GRU RNN with closely related LSTM lstm-1997 model. In the experimental section, we showed almost the same accuracy of GRU and LSTM given the fact that LSTM has a bigger amount of trainable parameters (it is slower).

5.3 Initial latent state approximation

Each of described above latent dynamics models requires an initial latent state $\hat{z}_{0}$ to start the computations for a scenario.

From the equation (3) we know, that initial reservoir state $f_{0}$ is a function of just metadata vector $m$: $f_{0}=F^{init}(m)$. Hence, $\hat{z}_{0}$ may be computed as follows:

[TABLE]

where $F^{init}$ is a known physics-based procedure aimed to calculate equilibristic initial reservoir state.

However, this approach requires explicit calculation of the high-dimensional tensor $f_{0}=F^{init}(m)$. Instead we may want to find the direct linear approximation for $\hat{z}_{0}$:

[TABLE]

6 Calculation of Well Production Rates

Besides the calculations of the dynamics, the metamodel is required to approximate the production rates function $R$.

To calculate the production rates from the well we use the physics-based model without trainable parameters. We exploit the fact that there are no additional sources or sinks in the development unit except for production and injection wells. Hence, the amount of oil, water, and gas produced might be calculated from the mass balance as the difference between subsequent reservoir states.

The volumetric amount of a fluid $\alpha$ at time $t$ can be calculated from the corresponding saturation field $s_{t}^{\alpha}$. Amounts of fluid from each computational cell should be brought to the same pressure $p_{ref}$ and then summed:

[TABLE]

where pressure and saturation distributions $p_{t}$, $s^{\alpha}_{t}$ are taken from the approximated reservoir state $\hat{f}_{t}$, $V_{bulk}$ is a bulk volume of computational cell, $\phi$ is porosity, $c_{r}$ and $c_{\alpha}$ are rock and fluid compressibilities, and the index $i$ denotes the number of a computational cell.

The reference pressure $p_{ref}$ should be chosen from the interval where there is no mass transfer involved between the oil and the gas phases. For example, it might be selected to be equal to the initial pressure at the upper part of the reservoir.

Thus we can calculate the liquid’s production rate at the reference pressure:

[TABLE]

where $q_{t}^{\alpha}$ is the injected amount of fluid $\alpha$, which is equal to zero for all fluids except for water.

We can now obtain the production rates at surface conditions as follows:

[TABLE]

where $b_{\alpha}$ is a formation volume factor of the fluid $\alpha$.

To take into account the influence of gas, liberated from the oil, we use the following modification of (17):

[TABLE]

where $\gamma$ is the gas content.

7 Summary

The overall structure of the metamodel consists of two distinct parts: dynamics and production rates models. The dynamics model itself decomposes into the dimensionality reduction model and the latent dynamics model.

In the experimental section (Section 8) we show the superiority of the dynamics model composed of VAE for dimensionality reduction and GRU RNN for latent dynamics forecasting against other investigated choices for these components of the dynamics model.

Algorithm 1 shows the overall training procedure for the metamodel and is suitable for different choices of metamodel’s components.

8 Results

To analyse the performance of proposed metamodelling approach, we generated a dataset (using the generative model described in B) and divided it into training $\mathbf{D}^{train}$ and validation $\mathbf{D}^{val}$ subsets, containing more than 5000 and 100 scenarios respectively. The dimensionality of the computational grid is chosen to be $(n_{x},n_{y},n_{z})=(41,60,10)$ and time discretization step: $\Delta t=14$ days.

The dataset contains scenarios of three types of possible development units, with either horizontal or vertical wells, different rock, and fluid properties, smooth or sharp applied control sequences $\mathbf{u}$.

3D ground truth solutions $\mathbf{f}$ and production rates $\mathbf{r}$ were obtained from Rock Flow Dynamics tNavigator simulator tnav used as the Base Model. We turned on the options of the simulator responsible for gas liberation and oil evaporation. “No flow” boundary conditions and equilibristic initial conditions were used throughout the simulations.

Following the Algorithm 1, we performed the training for two variations of the dimensionality reduction model: VAE and PCA, and four variations of latent dynamics model: LR, FCNN, GRU and LSTM RNNs. Two dimensionalities of the latent space were tested: $d_{z}=150$ and $d_{z}=250$.

We used ADAM adam-2014 optimizer to perform gradient-based optimization when needed. All the Neural Networks were implemented using standard NN modules from the PyTorch library pytorch for Python programming language. We used LReLU nonlinearity for all of the NNs, except the RNN modules where we use nonlinearities proposed in the original works gru-2014 ; lstm-1997 . Batch normalization batchnorm-2015 was used in all non-recurrent parts of NNs and layer normalization layrnorm-2016 in recurrent ones.

It should be mentioned that the trained metamodel needs no additional training to be applied to a new development unit. The ideal case for us is to construct a dataset broad enough, that the metamodel can be trained just once and then applied to arbitrary reservoirs.

Dynamics part of the metamodel was evaluated in the sence of relative deviation and Rooted Mean Squared Error (RMSE) between metamodel’s forecast and the ground truth for all scenarios from the validation dataset $\mathbf{D}^{val}$. The results of the analysis are provided in the Table 1 for different choices of metamodel’s components. Bold values denotes the best obtained approximation accuracy.

The FCNN model is not provided in the table due to its inability to output long predictions. The training scheme for FCNN is aimed to minimize the error of one-step prediction. This leads FCNN to the extreme instability when NN gets its outputs as inputs. Interestingly, that this effect does not appear in the LR model, which is structured similarly to FCNN model, due to the simplicity of linear regression.

Visual analysis of the dynamics forecasting performance for the 3D flow problem was made using the plots of mean pressure and fluid saturations versus time (Figure 4) and with the horizontal slices of pressure and saturation tensors at a fixed moment of time (Figure 5). Figures 4 and 5 show the results for the randomly chosen validation scenario with sharp applied control, staggered line drive allocation scheme, and vertical wells.

The accuracy of the procedure for computing production rates was tested on the best dynamics metamodel (VAE + GRU RNN, $d_{z}=150$) since it is strongly dependent on the dynamics accuracy. Mean relative deviation between predicted and ground truth rates for this model is provided on the Figure 6 where different groups of scenarios were analyzed separately. Visual check of the performance can be made on the plots of predicted daily production rates for randomly chosen validation scenario (Figure 7) – here we used the same scenario that was used in the dynamics visualization.

On average, the total on-wall time consumed by tNavigator to produce solutions for 30 years of simulation is 120 seconds. Computations were divided between 30 CPU cores. On the other hand, the metamodel consumes 10 seconds for the same simulation, if the GPU is used. Even in the case when the CPU is used, it consumes only 70 seconds. Though, convolutional layers are not adapted to be used with CPU. So, the metamodel is about x12 faster than the Base Model.

9 Discussion

We have developed the novel approach for prediction of reservoir performance enabling a new way of assessing the possible reservoir development scenario. Being much faster than a conventional reservoir modeling tool, the metamodel allows overlooking many more scenarios over a given time interval.

The approach is a pioneering one regarding its ability to output not only the flow rates at the wellheads, but also the actual dynamics of 3D fields of pressure, and fluid saturations over the computational cells of a conventional hydrodynamic model of a reservoir.

The approach itself is also somewhat unique because of its high-level structure containing methods of reduced order modelling (ROM) and the data-driven forecasting (DDF) tools. This high-level structure presents a very general framework for further enhancements as both ROM and the DDF methods are being developed these days very rapidly. We have tried just several ROM and DDF approaches and defined the ones demonstrating the best results.

There is no doubt that both building blocks of our workflow will progress in terms of accuracy and computational efficiency in training and especially prediction modes.

10 Further Work

We plan to continue our research in three principal directions. First is further experimenting with DDF and ROM architectures to increase the accuracy of the forecast and ensure the higher computational efficiency, as well as to try to combine data, generated by FDHS with different fidelities, when training a metamodel (see also MFGP2015 ; MFGP2017 ). Second is about going a scale up and try to adapt the metamodels for the full-scale reservoirs and oilfields containing many development units. Later will lead to more useful tools for applied reservoir engineering, including field development optimization and history matching. The third is the expansion of the metamodel’s operational envelope to the more complex fluid systems and their PVT representations and the low permeable reservoirs. The low permeable direction is likely to be supported by a training set generated with an approach presented in bezyan2019novel .

Appendix A Scenario metadata in compressed vector form

Variable $m$, defined in Section 3 and mentioned as metadata of the development unit, represents physical properties of the rock and fluids, and geometrical characteristics of the development unit and corresponding wells.

We use the metadata as an input of ML models, and most of them require inputs in a vectorised form with fixed dimensionality. Moreover, this dimensionality should be small to prevent the overfitting.

Some of the physical properties can be represented as scalars and used as is, being stacked into a vector. These are rock and water compressibilities, fluids densities, water viscosity, development unit size in each of three dimensions, initial water-oil and gas-oil contact depths, the depth of the top of the reservoir, initial pore pressure, saturation pressure and gas content.

Categorical properties, namely well allocation scheme and spatial orientation of wells, can be presented via one-hot encoding. The exact position of either vertical or horizontal well can be described by two values: the depth of the beginning of the perforation interval, the length of the perforation interval. We assume wells to have just one perforation interval.

The complicated part is to represent properties defined for each computational cell (porosity, absolute permeability) and table-defined functions (PVT tables, tables of relative permeabilities).

We exploit the restrictions imposed on the development unit to compress large fields such as porosity and permeability fields. We assume homogeneous rock in the sense of porosity, and layerwise-nonhomogeneous rock in the sense of permeability. So the porosity can be represented as a scalar, and permeability field forms the vector of length $n_{z}$ with values for each layer of the reservoir.

The problem with table-defined functions is that they have a non-fixed number of points in their tables. The problem can be resolved by interpolating these functions onto a predefined domain.

Such representation fixes the dimensionality but is rather high-dimensional. We compress the dimensionality using the PCA method.

Actually, the restrictions on the homogeneity of the development unit might be overcome using the same dimensionality reduction idea (possible with the neural network instead of PCA). Though, this approach requires construction of the generative model for the whole heterogeneous permeability field, which is a complicated task by itself.

Appendix B Generative model of realistic metadata and control variables

Generation of synthetic but realistic oil-fields is a very complex task by itself. We restricted ourselves on the generation of just pieces of an oil-field (development units) with the layered structure.

But still, the detailed description of the process is too large to be shown in the paper. The generative model should provide us with the randomised samples of the reservoir geometry, porosity and absolute permeability fields, initial conditions, PVT and relative permeability tables, exact positions of the wells, and so on.

Here we present the overall ideas used to obtain the generative model:

• •

Collect a large set of values of the interest from a variety of real oil-fields.

• •

Divide all the properties into interindependent groups.

• •

Visualize joint empirical distributions for each group.

• •

For each group, fit parameters of chosen parametric distribution to match the empirical data using the max-likelihood method. We used following parametric families of distributions: Gaussian, log-Gaussian, Uniform.

For example, we grouped values of porosity and absolute permeability and modelled them jointly. The form of the empirical distribution is similar to two dimensional Gaussian if the permeabilities are logarithmized. Samples from the resulting distribution are shown in Figure 8 together with the points from the dataset.

Table-defined functions are generated from the point estimates paired with standard correlation functions, such as Corey, Al-Marhoun, Khan et al., Vasquez and Beggs correlations [33].

For the time-series (bottomhole pressures and injection rates) we used three possible forms: constant, exponential, sharp (see Fig. 9). The form is sampled randomly, and actual values are obtained from the Uniform distribution.