[orcid=0000-0003-3065-8244] [1]

# Combined mechanistic and machine learning method for construction of oil
reservoir permeability map consistent with well test measurements

E.A. Kanin<EMAIL_ADDRESS>A.A. Garipova  S.A. Boronin  V.V.
Vanovsky  A.L. Vainshtein  A.A. Afanasyev  A.A. Osiptsov  E.V. Burnaev
Skolkovo Institute of Science and Technology (Skoltech), Bolshoy Boulevard 30,
bld. 1, Moscow, Russia 121205  Institute of Mechanics, Moscow State
University, Michurinsky Pr., 1, Moscow, Russia 119192

###### Abstract

We propose a new method for construction of the absolute permeability map
consistent with the interpreted results of well logging and well test
measurements in oil reservoirs. Nadaraya-Watson kernel regression is used to
approximate two-dimensional spatial distribution of the rock permeability.
Parameters of the kernel regression are tuned by solving the optimization
problem in which, for each well placed in an oil reservoir, we minimize the
difference between the actual and predicted values of (i) absolute
permeability at the well location (results of interpretation of well logging);
(ii) absolute integral permeability of the domain around the well and (iii)
skin factor (results of interpretation of well tests). Optimization task
(inverse problem) is solved via multiple solutions to forward problems, in
which we estimate the integral permeability of reservoir surrounding a well
and the skin factor by the surrogate model. The last one is developed using an
artificial neural network trained on the physics-based synthetic dataset
generated using the procedure comprising the numerical simulation of
bottomhole pressure decline curve in reservoir simulator followed by its
interpretation using a semi-analytical reservoir model. The developed method
for reservoir permeability map construction is applied to the available
reservoir model (Egg Model) with highly heterogeneous permeability
distribution due to the presence of highly-permeable channels. We showed that
the constructed permeability map is hydrodynamically similar to the original
one. Numerical simulations of production in the reservoir with constructed and
original permeability maps are quantitatively similar in terms of the pore
pressure and fluid saturations distribution at the end of the simulation
period. Moreover, we obtained an good match between the obtained results of
numerical simulations in terms of the flow rates and total volumes of produced
oil, water and injected water.

###### keywords:

absolute permeability well logging well test hydrodynamic modeling machine
learning artificial neural network optimization algorithms

##  1 Introduction

Building high-quality reservoir simulation models remains a complex task that
requires the synergy of several branches of geoscience and reservoir
engineering. Traditional approaches to geological model construction, in
particular, stochastic modeling, are very time-consuming with the main
disadvantage being uncertainty in resulting reservoir properties. For reliable
production simulation results, petroleum engineers have to solve the inverse
problem, namely, the history matching of the hydrodynamic model. It is an
iterative calibration process involving the alteration of the parameters of
the original geological model to match the production data. During geological
modeling, 2D maps of absolute or effective reservoir permeability are built
using static datasets at well locations, which consist of information obtained
usually from well logging and core studies (routine and specific core
analysis). Accounting for dynamic well test interpretation data is carried out
by various techniques but still poses a challenge for engineers.

The most common way to bring well test data in line with the static data is to
adjust the petrophysical relationship between rock porosity and permeability
or to tune the variogram used for spatial correlation. These procedures
involve manual adjustments according to experience and expertise of a
particular engineer. Therefore, the entire process can take a long time and
the result can be not optimal in view of hydrodynamic similarity in between
the real reservoir property distribution and the constructed permeability and
porosity map.

Zakirov et al. (2014, 2016, 2018) presented a highly-promising approach to
solve the problem of well logging and well test data fusion in constructing
the reservoir permeability map. The authors performed the geostatistically
driven history matching by adjoint methods with the set of control parameters
including properties of variogram, porosity-permeability relationship for each
rock facies, and data at control points. The algorithm is automated and has
been successfully validated at several synthetic cases and applied for
realistic oilfield models.

Kolesnikov et al. (2010) suggested calculating productivity index avoiding
reservoir simulations by finite difference approach, which resembles tensor
permeability upscaling methods for the rapid calibration of the geological
models with well test measurements and production profiles. Modifications of
permeability used to calculate productivity indexes allow matching
productivity indexes in geological models to the results of well test
interpretation. The authors tested the approach on the geological model of an
oilfield located in West Siberia. Comparison of actual productivity indexes
obtained using well test data with those calculated by the reservoir model
after its execution showed an acceptable agreement.

In the study by He et al. (2000), authors described a multi-step procedure for
the efficient generation of reservoir properties accounting for the dynamic
data obtained using stochastic models. Authors claimed that the set of
realizations obtained using this algorithm typically provides an acceptable
approximation to the probability density function for reservoir models so that
the proposed method can be used in Monte Carlo modeling. Authors reported that
their approach allowed one to comply with the well test data and retain a
heterogeneity typical of the original geological model. It consists of two
parts, namely, automatic and manual. The former part is a minimization problem
and the latter one is a data preparation task, which requires: (i) cutting a
sector from the reservoir model to apply a history matching process; (ii)
upscaling the sector model; (iii) determining a variogram using the coarse
model, and (iv) choosing the sets of variable and fixed parameters. Even
though the automatic part was quite efficient, the procedure could not be
utilized as a “production line” at that stage.

The method based on the Ensemble Kalman Filter (EnKF) is presented in Coutinho
et al. (2010) with the aim to update reservoir permeability distribution using
available bottomhole pressure profiles and well logging data, estimate skin
factors of reservoir layers and compute “effective” skin factor of wells in
the framework of multilayered reservoir models. The algorithm worked well for
synthetic cases, but it was not successful when applied to the real field
case. The authors considered two different approaches: updating layer
permeability multipliers with EnKF and double stochastic EnKF, but neither of
the two approaches provided a reasonable data match.

Ensemble Kalman Filter was successfully used to perform history matching for
real field cases. Evensen et al. (2007) applied this technique to a reservoir
located in the North Sea area to estimate permeability, porosity, initial
fluid contacts as well as vertical and fault transmissivity multipliers. A
similar approach is applied to oil saturated reservoir in the paper (Bianco et
al., 2007). The EnKF was also used to conduct history matching for a deepwater
formation with multiscale parametrization (Zhang and Oliver, 2011).

The evolution of EnKF based methods lead to the development of ensemble
smoothers (ES) and their use in history matching. In the paper by Evensen and
Eikrem (2018), authors discussed the algorithms for adjusting reservoir models
to comply with production rate data by ensemble smoothers. The advantage is
that ES allows considering so-called hyperparameters, which represent
geological model inputs and less computational expenses compared to EnKF. The
authors discussed approaches to reduce redundant information in production
data series and how to deal with errors correlated in time. Real field
application results have shown that the Iterative ensemble Smoother
formulation performed the best in comparison with other formulations.

As we can see, automated conditioning of dynamic well data is a complex task,
and the search for reliable algorithms is still ongoing. Machine learning
methods can be applied to significantly speed-up different stages of modeling
workflow, for example, seismic and well logging interpretation, facial and
petrophysical analysis. Developments in artificial intelligence technologies,
especially neural networks, allowed considering the data fusion problems from
a different perspective. Several studies were published in the last couple
years and dealt with the approaches based on neural networks as described
below.

In the study by Thanh and Sugai (2021), the authors proposed an enhanced
framework for modeling the distribution of lithofacies and petrophysical
properties of a sandstone reservoir containing fluvial channels. The
integrated method is proposed, which is based on (i) artificial neural network
(ANN); (ii) Sequential Gaussian Simulation, and (iii) object-based modeling.
ANN is used to predict the petrophysical properties of the reservoir by
combining seismic attributes and well logging data. Object-based modeling was
applied to distribute facies of channels in the 3D model, which allows for
resolving realistic depositional environments. The proposed modeling workflow
facilitates the reduction of the typical time required to conduct the history
matching procedure applied to a reservoir containing fluvial channels.

Bai and Tahmasebi (2020) suggested the algorithm to construct geologically
realistic subsurface models conditioned to well location data with the help of
a surrogate algorithm. While cross-correlation-based simulation (CCSIM) allows
effective reconstruction of reservoir models, it suffers from the requirement
to balance between the quality of realizations and the degree of point data
reproduction. The authors combined a pattern-based method with a convolutional
neural network (CNN) to overcome this challenge. Inside the surrogate
algorithm combining CCSIM and CNN, the former was used for grids in the
absence of real data as it allows the generation of high-quality geological
features. For the grids, where real data is available, CCSIM is utilized to
obtain an initial guess, while CNN is applied to improve initial realizations
and increase the accuracy of the real data reproduction. Using the spatial
distribution of real data, the proposed method allows one to determine missing
domains to cover the mismatched real data. Further, a refill of initial model
realizations containing missing zones is carried out using the trained model.

Titus et al. (2022) used ANNs to condition a surface-based geological model
(SBGM), constructed with a parametric non-uniform rational B-spline (NURBS)
approach to well data. ANNs were applied in the following way: (i) to map
input parameters of SGBM to types of facies in the vicinity of well locations
in the framework of the forward modeling step and (ii) to obtain the optimized
set of input parameters of SBGM using a back-propagation method, so that the
constructed SBGM complies with the types of facies obtained in well
measurements. The approach was tested on a synthetic 2D case, and it
demonstrated the ability to generate a set of realizations that matches the
well data. Moving from 2D facies distribution towards petrophysical
properties, the authors observed the potential of CNNs and RNNs (recurrent
neural networks): the former allows one to learn important features of
spatially correlated data, while the latter can be used to condition
individual surfaces to model temporal sequences. Moreover, the authors
mentioned that the proposed methodology can be reapplied to object-based
models, in which a complex reservoir geometry is described by parameterized
objects.

As we can see from the literature review, machine learning algorithms has been
successfully applied to the problem of production data conditioning and
history matching even though it is still a developing field. Currently there
is a gap in between existing methods of permeability map construction based on
machine learning and well test and well log data fusion. We believe that this
important component of geological model construction can be successfully
solved provided the corresponding artificial intelligence tool is developed.

We would like to highlight the importance to construct the absolute
permeability field approximation that is hydrodynamically similar to the
actual distribution around each well. In the opposite case, when the
similarity is absent, one can obtain the significant difference between
modeling results and production history when the approximate absolute
permeability cube is embedded into a hydrodynamic simulator. Permeability
distribution accounting for well test interpretation results can be an optimal
initial approximation for the history matching procedure. The main aim of the
present work is to develop a computationally efficient algorithm for
constructing the absolute permeability cube based on well logging and well
test data fusion. We apply numerical and semi-analytical hydrodynamic
modeling, and optimization algorithms. Computationally heavy components of
this chain of numerical algorithms, namely, reservoir production simulations
to obtain pressure decline curve and its consequent interpretation using the
analytical reservoir model, are replaced by a fast surrogate model developed
using machine learning (ML) algorithm. We demonstrate the capabilities of the
proposed combined mechanistic and ML approach using the synthetic case, while
the developed method can be used for field data with no modifications. The
proposed surrogate algorithm of well test and well log data fusion can be
implemented into a wide variety of existing algorithms of permeability map
construction to improve the initial guess to overall history matching process
by preserving hydrodynamic similarity in between the original and constructed
maps.

We organize the paper in the following way. Section 2 outlines the problem
formulation. Section 3 describes the methodology for building the absolute
permeability field approximation. In Section 4, we explain the procedure of
synthetic dataset generation using the numerical hydrodynamic simulator and
semi-analytical reservoir model and describe the surrogate model predicting
the integral permeability around the well and skin factor. Section 5
demonstrates the obtained results and their analyses. Section 6 provides the
main features and limitations of the proposed models for the construction of
the absolute permeability map as well as the directions for future
development. Finally, we summarize the main findings of our study and provide
potential directions for further research in the area in Section 7.

##  2 Problem formulation

The proposed methodology for building the absolute permeability field
approximation is outlined in Section 3 based on the example of the synthetic
reservoir model which we describe in the current Section. Figure 1 shows the
schematic drawing of the reservoir model.

Figure 1:  A synthetic reservoir model; $N$ vertical wells (blue cylinders)
cross the formation by perforating it along the entire thickness (red colored
zones).

Reservoir is approximated by a box with dimensions $\Delta X\times\Delta
Y\times\Delta Z$. It is fully penetrated by $N$ vertical wells which are
parallel to $z$-axis and located at the points
$\left\\{x_{i},y_{i}\right\\}\big{|}_{i=1}^{N}$. We study the formation with
homogeneous filtration-storage properties along the vertical direction. As a
result, the discussed model is two-dimensional, and our target reservoir
property is a function of the lateral coordinates only, $k_{0}=k_{0}(x,y)$,
where the subscript “0” denotes the real absolute permeability distribution,
which is needed to be approximated.

Assume that the absolute permeabilities ${k^{\mathrm{WL}}_{i}}$,
${k^{\mathrm{WT}}_{i}}$ and skin factor $S_{i}$ are known for each well.
Parameter ${k^{\mathrm{WL}}_{i}}$ corresponds to the absolute permeability at
the well location, so that ${k^{\mathrm{WL}}_{i}}=k_{0}(x_{i},y_{i})$ and
superscript “WL” stands for ’well logging’. Using the well logging
measurements, one can determine the rock porosity $\phi$ in the vicinity of
the well, which can be converted into absolute permeability via the dependence
$k^{\mathrm{WL}}(\phi)$ obtained using lab experiments on rock core samples.
Parameter ${k^{\mathrm{WT}}_{i}}$ is an integral permeability of the rock
surrounding the the well, so that

${k^{\mathrm{WT}}_{i}}=\mathcal{F}(k_{0}(x,y)),~{}(x,y)\in\mathcal{A}^{\mathcal{R}}_{i},$
(1)

where $\mathcal{A}^{\mathcal{R}}_{i}$ is a circle round the well $i$ of radius
$\mathcal{R}$, and superscript “WT” denotes ‘well test’. By function
$\mathcal{F}$ we denote the physics-based averaging of the permeability field
around the well inside the circle of radius $\mathcal{R}$, which can be
interpreted as a distance from the well up to which the pore pressure
disturbance propagate during the well test. Skin factor $S_{i}$ quantifying
the contrast between the rock permeability in the very vicinity of the well
and at a certain distance from it, is also a function of the absolute
permeability field around the well

$S_{i}=\mathcal{G}(k_{0}(x,y)),~{}(x,y)\in\mathcal{A}^{\mathcal{R}}_{i}.$ (2)

The absolute integral permeability ${k^{\mathrm{WT}}_{i}}$ and skin factor
$S_{i}$ can be found from the well test analysis. The standard method is
interpretation of the bottomhole pressure build-up curve during the well shut-
in (build-up test) (Horne, 1995). Using the approach proposed by Perrine
(1956) and Martin (1959) including the concepts of total mobility and total
compressibility, one can estimate the total mobility $\lambda_{t}$ and skin
factor values from the interpretation. Further, absolute permeability is
computed using the equation

$k^{\mathrm{WT}}=\lambda_{t}\left(\frac{k_{ro}}{\mu_{o}}+\frac{k_{rw}}{\mu_{w}}+\frac{k_{rg}}{\mu_{g}}\right)^{-1},$

where $k_{ro}$, $\mu_{o}$, $k_{rw}$, $\mu_{w}$, $k_{rg}$, $\mu_{g}$, are
relative permeabilities and viscosities of oil, water, and gas phase,
respectively. Saturations of each phase in the reservoir are required to
compute $k_{ro},k_{rw},k_{rg}$. Since saturations depend on time and distance
from the well during test, and there are certain difficulties in their
measurement, $k_{\mathrm{WT}}$ value is usually determined with significant
error.

In the current work, we construct the approximation $k=k(x,y)$ of the real
absolute permeability field $k_{0}=k_{0}(x,y)$ assuming that the well
locations $\left\\{x_{i},y_{i}\right\\}\big{|}_{i=1}^{N}$, permeability values
obtained from the well logging
$\left\\{k^{\mathrm{WL}}_{i}\right\\}\big{|}_{i=1}^{N}$, integral
permeabilities $\left\\{k^{\mathrm{WT}}_{i}\right\\}\big{|}_{i=1}^{N}$, and
skin factor values $\left\\{S_{i}\right\\}\big{|}_{i=1}^{N}$ are available.
Since we account for the integral permeability in the approximation, the
absolute permeability map around each well is hydrodynamically similar to the
real distribution, which allows to obtain acceptable match between the results
of reservoir simulations and observed production data already at the start of
history matching process.

##  3 Global model for reservoir permeability map construction

In the framework of current study we approximate the absolute permeability map
of reservoir $k_{0}=k_{0}(x,y)$ using Nadaraya-Watson kernel regression
(Nadaraya, 1964; Watson, 1964). According to this approach, functional
dependence of absolute permeability on the spatial coordinates has the
following form:

$k(x,y)=\frac{\sum_{i=1}^{N}\left[k^{\mathrm{near}}_{i}\mathcal{K}^{\mathrm{near}}(r-r_{i})+k^{\mathrm{far}}_{i}\mathcal{K}^{\mathrm{far}}(r-r_{i})\right]}{\sum_{i=1}^{N}\left[\mathcal{K}^{\mathrm{near}}(r-r_{i})+\mathcal{K}^{\mathrm{far}}(r-r_{i})\right]},$
(3)

where $r=\sqrt{x^{2}+y^{2}}$, $r_{i}=\sqrt{x_{i}^{2}+y_{i}^{2}}$ and kernel
functions:

$\mathcal{K}^{\mathrm{far}}(r)=\left(\frac{r}{r_{d}}\right)^{\alpha}\mathrm{exp}\left[-\left(\frac{r}{r_{d}}\right)^{\beta}\right],~{}~{}~{}~{}\mathcal{K}^{\mathrm{near}}(r)=\gamma~{}\mathrm{exp}\left[-\left(\frac{r}{r_{g}}\right)^{\delta}\right]$
(4)

Equations (3) and (4) contain $2N+6$ unknown parameters:

$\Phi=\left\\{\Phi_{i}\right\\}\big{|}_{i=1}^{N}=\left\\{k^{\mathrm{near}}_{i},k^{\mathrm{far}}_{i}\right\\}\big{|}_{i=1}^{N},\quad\Psi=\left\\{r_{d},r_{g},\alpha,\beta,\gamma,\delta\right\\},$

which we denote by $\Omega=\\{\Phi,\Psi\\}$. The former ones, $\Phi$, have the
dimension of permeability (in millidarcy range in the framework of current
study), and describe the contribution of permeabilities in the near and far
zones of each well to the reservoir absolute permeability map $k(x,y)$
calculated using Eq. (3). The contributions are weighted by the kernel
functions (4) parameterized by characteristic lengthscales $r_{d}$ and
$r_{g}$, which have the dimension of length, as well as dimensionless
variables $\alpha,\beta,\gamma,\delta$, which are the fixed and assumed to be
similar of all wells. A typical spacial distribution of kernel functions (10)
is shown in Fig. 1. One can observe that function
$\mathcal{K}^{\mathrm{near}}(r)$ takes maximum value at $r=0$, so that its the
main contribution to overall permeability map $k(x,\,y)$ is at the point of
the well location $(x_{i},y_{i})$. At the same time,
$\mathcal{K}^{\mathrm{far}}(r)$ reaches the maximum at the distance
$r=(\alpha/\beta)^{1/\beta}r_{d}$, which describes the far-field contribution
to absolute permeability map according to well test data at the distance
$r_{d}$ scaled by the combination of parameters $\alpha$ and $\beta$.

Figure 2:  Kernel functions $\mathcal{K}^{\mathrm{far}}(r)$ and
$\mathcal{K}^{\mathrm{near}}(r)$ (see Eq. (10)) at $r_{d}=150~{}\mathrm{m}$,
$r_{g}=30~{}\mathrm{m}$, $\alpha=1,~{}\beta=1,~{}\gamma=0.5,~{}\delta=0.5$.

Parameters $\Omega$ are calculated using the available formation properties
described in Section 2:
$\left\\{k^{\mathrm{WL}}_{i},k^{\mathrm{WT}}_{i},S_{i}\right\\}\big{|}_{i=1}^{N}$.
They are obtained via solution of the minimization problem:

$\min_{\Omega}\frac{1}{N}\sum_{i=1}^{N}\big{|}\big{|}\widetilde{\mathcal{X}}^{*}_{i}-\mathcal{X}^{*}_{i}\big{|}\big{|}_{1},$
(5)

where $||\cdot||_{1}$ denotes the $L_{1}$-norm. Vector $\mathcal{X}_{i}$ is
composed of the “true” values of absolute permeabilities obtained via well
logging and well tests as well as skin factors:
$\mathcal{X}_{i}=\left\\{k^{\mathrm{WL}}_{i},k^{\mathrm{WT}}_{i},S_{i}\right\\}$;
vector $\widetilde{\mathcal{X}}$ is composed of similar parameters estimated
using the approximate absolute permeability field (3),
$\widetilde{\mathcal{X}}_{i}=\left\\{k(x_{i},y_{i}),\widetilde{k}_{i},\widetilde{S}_{i}\right\\}$;
the superscript “*” denotes that the parameters are scaled before substituting
into target function (5). We apply standard scaling, so that vector has zero
mean and unit variance after the transformation with the parameters determined
from the normalization of matrix
$\mathcal{X}=\left\\{\mathcal{X}_{i}\right\\}\big{|}_{i=1}^{N}$. We solve the
minimization problem (5) using differential evolution optimization algorithm
implemented in SciPy library (Virtanen et al., 2020).

We predict the values of integral permeability
$\\{\widetilde{k}_{i}\\}_{i=1}^{N}$ and skin factor
$\\{\widetilde{S}_{i}\\}_{i=1}^{N}$ for each well using the surrogate model.
In fact, this model approximates functions $\mathcal{F}$ and $\mathcal{G}$
introduced in Section 2 (see Eqs. (1), (2)). Details of the development of the
surrogate model are formulated below in Section 4, while in the rest of this
section we discuss application of the developed model to the construction of
absolute permeability map.

Let us consider the well denoted by index $j$. We clip a square domain of size
$\mathcal{D}$ with sides parallel to the axes $x$ and $y$, which is formally
described as

$\mathcal{B}^{\mathcal{D}}_{j}=\left[x_{j}-\frac{\mathcal{D}}{2},x_{j}+\frac{\mathcal{D}}{2}\right]\times\left[y_{j}-\frac{\mathcal{D}}{2},y_{j}+\frac{\mathcal{D}}{2}\right].$

The absolute permeability distribution is given by equation (3), which we
rewrite in the following way

$\displaystyle k(x,y)=\frac{\mathcal{M}_{j}+\sum_{i=1,i\neq
j}^{N}\mathcal{M}_{i}}{\mathcal{N}_{j}+\sum_{i=1,i\neq
j}^{N}\mathcal{N}_{i}},$
$\displaystyle\mathcal{M}_{i}=k^{\mathrm{near}}_{i}\mathcal{K}^{\mathrm{near}}(r-r_{i})+k^{\mathrm{far}}_{i}\mathcal{K}^{\mathrm{far}}(r-r_{i}),$
$\displaystyle\mathcal{N}_{i}=\mathcal{K}^{\mathrm{near}}(r-r_{i})+\mathcal{K}^{\mathrm{far}}(r-r_{i}).$
(6)

Assuming that the distance between the wells is large enough, the contribution
from well $j$ into the absolute permeability field $k(x,y)$ is given by the
ratio $k^{\mathrm{well}}_{j}=\mathcal{M}_{j}/\mathcal{N}_{j}$, which is
described by the set of parameters $\\{\Phi_{j},\Psi\\}$. In turn, the impact
of the neighbouring wells on the permeability distribution $k(x,y)$ is
governed by the function $k^{\mathrm{neigh}}_{j}=\sum_{i=1,i\neq
j}^{N}\mathcal{M}_{i}/\sum_{i=1,i\neq j}^{N}\mathcal{N}_{i}$. One can
approximate the function $k^{\mathrm{neigh}}_{j}$ inside the zone
$\mathcal{B}^{\mathcal{D}}_{j}$ by a quadratic polynomial with two variables:

$\displaystyle\widetilde{k}^{\mathrm{neigh}}_{j}=a_{2,0}(x-x_{j})^{2}+a_{1,1}(x-x_{j})(y-y_{j})+a_{0,2}(y-y_{j})^{2}+$
$\displaystyle+a_{1,0}(x-x_{j})+a_{0,1}(y-y_{j})+a_{0,0}.$ (7)

Consequently, there is a one-to-one correspondence between the permeability
distribution around the well $j$ described as
$k(x,y),~{}(x,y)\in\mathcal{B}^{\mathcal{D}}_{j}$, and the set of 14
parameters:

$\chi=\left\\{\underbrace{k^{\mathrm{near}}_{j},k^{\mathrm{far}}_{j},}_{\Phi_{j}}\underbrace{r_{d},r_{g},\alpha,\beta,\gamma,\delta,}_{\Psi}\underbrace{a_{2,0},a_{1,1},a_{0,2},a_{1,0},a_{0,1},a_{0,0}}_{\Xi_{j}}\right\\},$
(8)

where we introduce notation $\Xi_{j}$ for the quadratic polynomial
coefficients. Based on the values of the parameters (8), our surrogate model
(Section 4) estimates the integral permeability $\widetilde{k}_{j}$ and the
skin factor $\widetilde{S}_{j}$.

##  4 Surrogate model for evaluation of integral permeability

In this section, we describe the surrogate model predicting the integral
permeability around the well (mimicking the results of well test
interpretation) and corresponding skin factor. The model is based on machine
learning algorithm, namely, artificial neural network (ANN) (Rosenblatt,
1958). ANN is trained on the physics-based synthetic dataset generated with
the help of the numerical hydrodynamic simulator MUFITS (Afanasyev, 2020) and
in-house semi-analytical reservoir model (Ozkan and Raghavan, 1991a).

We begin with the explanation of the synthetic dataset preparation procedure.
We create multiple instances of synthetic reservoir model similar to that
described in Section 2 (Figure 1), namely, a rectangular box with dimensions
$\Delta X=\Delta Y=4~{}\mathrm{km}$, $\Delta Z=10~{}\mathrm{m}$ (the thickness
value is for reference only since the model is effectively two-dimensional).

Let us consider the generation of a single realization of the synthetic
reservoir model. Vertical wells are placed randomly subject to the condition
that the distance between any two wells is greater then the predefined value
$\Delta d$, so that $\sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}}\geq\Delta d$
for any $i\neq j$. For a particular reservoir model, distance $\Delta d$ is
chosen randomly in the range $\Delta d\in[300,600]~{}\mathrm{m}$. The
placement of vertical wells continues until the algorithm can not find any
possible location for the next well. Figure 3 illustrates the results of
application of the described procedure at $\Delta d$ = 500 m.

Figure 3:  An example of the synthetic reservoir model generated according to
the algorithm described in the main text. Here, the minimal distance between
any two wells $\Delta d$ equals 500 m. Black crosses denote the positions of
the vertical wells, while the squares with dotted blue boundary mark the
domains $\mathcal{B}^{\mathcal{D}}$ used for the creation of the samples for
the synthetic dataset.

At the next step, the algorithm assigns the set of parameters for each well
$\left\\{\Phi_{i},\Psi_{i}\right\\}_{i=1}^{N}$, where $\Psi_{i}$ is the set of
kernel parameters $\\{r_{d},r_{g},\alpha,\beta,\gamma,\delta\\}$ specified for
each of the wells. This feature of the absolute permeability map
parametrization differs from the one utilized in the general modelling
approach (Section 3) allowing us to use diversified samples from a single
synthetic reservoir model into the overall dataset. The values of parameters
$\left\\{\Phi_{i},\Psi_{i}\right\\}_{i=1}^{N}$ are selected randomly in the
following ranges:

$\displaystyle
k^{\mathrm{near}},k^{\mathrm{far}}\in[1,15]~{}\mathrm{mD},~{}r_{d}\in[100,300]~{}\mathrm{m},~{}r_{g}\in[5,50]~{}\mathrm{m},$
$\displaystyle\alpha\in[0.5,2],~{}\beta\in[1,2],~{}\gamma\in[0.01,2],~{}\delta\in[0.05,1].$
(9)

The range for variables $k^{\mathrm{near}},k^{\mathrm{far}}$ are chosen in
accordance with the typical permeability values of an oilfield in Western
Siberia, while the intervals for the geometrical parameters $r_{d},r_{g}$ are
chosen according to the typical areas around wells covered by well logging and
well tests. We select experimentally the ranges for the remaining parameters
of the kernel regression, namely, $\alpha,\beta,\gamma,\delta$.

Next, we apply the modified version of kernel regression (3) to compute the
absolute permeability field:

$\displaystyle
k(x,y)=\frac{\sum_{i=1}^{N}\left[k^{\mathrm{near}}_{i}\mathcal{K}^{\mathrm{near}}_{i}(r-r_{i})+k^{\mathrm{far}}_{i}\mathcal{K}^{\mathrm{far}}_{i}(r-r_{i})\right]}{\sum_{i=1}^{N}\left[\mathcal{K}^{\mathrm{near}}_{i}(r-r_{i})+\mathcal{K}^{\mathrm{far}}_{i}(r-r_{i})\right]},$
$\displaystyle\mathcal{K}^{\mathrm{far}}_{i}(r)=\left(\frac{r}{r_{d,i}}\right)^{\alpha_{i}}\mathrm{exp}\left[-\left(\frac{r}{r_{d,i}}\right)^{\beta_{i}}\right],$
$\displaystyle\mathcal{K}^{\mathrm{near}}_{i}(r)=\gamma_{i}~{}\mathrm{exp}\left[-\left(\frac{r}{r_{g,i}}\right)^{\delta_{i}}\right]$
(10)

After that, the algorithm is applied to each well and performs the following
operations. It cuts a square of size $\mathcal{D}=500~{}\mathrm{m}$ with sides
parallel to the coordinate axes (domain $\mathcal{B}_{j}^{\mathcal{D}}$). The
algorithm identifies the contributions from well $j$ ($k_{j}^{\mathrm{well}}$,
see Figure 4b) and neighbor wells ($k_{j}^{\mathrm{neigh}}$, see Figure 4c)
into the absolute permeability field
$k(x,y),(x,y)\in\mathcal{B}_{j}^{\mathcal{D}}$ (see Figure 4a). The former
distribution is characterized by parameters $\\{\Phi_{j},\Psi_{j}\\}$, while
the latter one is governed by the quadratic polynomial coefficients (7); the
approximation of $k_{j}^{\mathrm{neigh}}$ by the polynomial is shown in Figure
4d.

Figure 4:  Plot (a) shows an example of the synthetic permeability
distribution $k(x,y)$ in the square domain $\mathcal{B}^{\mathcal{D}}$;
contributions from well $k^{\mathrm{well}}$ located in the center of the
domain and neighbor wells $k^{\mathrm{neigh}}$ are shown in plots (b) and (c);
distribution $k^{\mathrm{neigh}}$ is fitted by the quadratic polynomial (7),
and plot (d) shows the result of approximation.

Finally, we save the absolute permeability distribution
$k(x,y),(x,y)\in\mathcal{B}_{j}^{\mathcal{D}}$ into the text file for
consequent numerical simulations in MUFITS, while the vector
$\chi=\left\\{\Phi_{j},\Phi_{j},\Xi_{j}\right\\}$ is stored into the table of
input parameters. This completes the input features preparation procedure for
the surrogate model.

Now, we discuss the estimation of the output parameters for neural network,
namely, integral permeability and skin factor. The procedure can be called as
the synthetic well test, and it consists of two stages:

  1. 1.

the numerical simulations of a drawdown test: bottomhole pressure dynamics
$p^{\mathrm{num}}_{w}$ is simulated using MUFITS provided the fixed well flow
rate;

  2. 2.

interpretation of the bottomhole pressure dynamics using the semi-analytical
reservoir model; we minimize the $L_{2}$-norm of the vector composed of the
deviations between the numerical and semi-analytical $p^{\mathrm{s-a}}_{w}(t)$
pressure values in different time instants within the specified interval:

$\min_{\widetilde{k},\widetilde{S}}\big{|}\big{|}\mathcal{P}\big{|}\big{|}_{2},$
(11)

where
$\mathcal{P}=\\{p^{\mathrm{num}}_{w}(t_{i})-p^{\mathrm{s-a}}_{w}(t_{i})\\}_{i=1}^{T}$.
The minimization problem is solved using the gradient method Nelder–Mead
implemented in SciPy library (Virtanen et al., 2020).

We specify the following input parameters in the numerical and semi-analytical
hydrodynamic models:

  * •

formation dimensions are $\Delta x=\Delta y=500$ m, $\Delta z=10$ m;

  * •

well radius is $r_{w}=0.1$ m;

  * •

porosity is $\phi=15~{}\%$;

  * •

fluid parameters: viscosity is $\mu=2.5$ cP, total compressibility is
$c_{t}=2\cdot 10^{-4}$ bar-1, formation volume factor is $B=1.2$ m3/sm3;

  * •

initial condition is uniform pore pressure $p_{i}=250$ bar;

  * •

boundary conditions are constant pressure $p=p_{i}$ at lateral borders, top
and bottom boundaries are closed;

  * •

flow rate is $q=20$ m3/d;

  * •

production period is $t\in[0,180]$ d.

Each permeability distribution $k(x,y),(x,y)\in\mathcal{B}^{\mathcal{D}}$
obtained using Eq. (10) and stored in a text file is passed to numerical
hydrodynamic simulator. In the semi-analytical hydrodynamic model,
permeability $\widetilde{k}$ is uniform and its value, as well as the skin
factor $\widetilde{S}$, are determined by the solution of the minimization
problem (11). Consequently, $\widetilde{k}$ parameter is assumed to be the
integral permeability corresponding to the distribution $k(x,y)$ inside the
domain $\mathcal{B}^{\mathcal{D}}$, and the adjusted homogeneous permeability
distribution is hydrodynamically similar to the heterogeneous absolute
permeability field $k(x,y)$. Note that the described physics-based averaging
procedure for the absolute permeability field does not depend on fluid and
rock properties, boundary conditions and operation mode since the absolute
permeability is a geometrical parameter of the rock. Integral permeability is
determined according to the bottomhole pressure dynamics during the transient
production period at the interpretation stage.

In the numerical reservoir model, approximation is carried out using the mesh
with cell size of 15 m $\times$ 15 m. Local grid refinement is applied in the
vicinity of the wellbore with the cell size of 5 m. Note the the choice of
mesh resolution is based on numerical convergence tests. Filtration in the
reservoir is simulated using BLACKOIL module, and we enable the single-phase
fluid option, namely, a dead oil.

Semi-analytical reservoir model is represented by a fully-penetrating vertical
line-source well in a reservoir, the analytical solution is carried out in the
Laplace space. The inverse Laplace transformation is performed using Stehfest
numerical algorithm (Stehfest, 1970). To give some details on the analytical
solution, it is derived using the principle of superposition according to
which, the solution for a uniform-flux point source along the well trajectory
is integrated (Ozkan, 1988; Ozkan and Raghavan, 1991a). The latter function is
an analytical solution to 3D filtration equation for slightly compressible
fluid in the reservoir approximated by a rectangular box with homogeneous
properties, in which the uniform-flux point source is located. The
mathematical formulation of the point-source problem is as follows (Van
Everdingen and Hurst, 1949; Hovanessian, 1961):

$\displaystyle\frac{\widetilde{k}}{\mu}\left[\frac{\partial^{2}p}{\partial
x^{2}}+\frac{\partial^{2}p}{\partial y^{2}}+\frac{\partial^{2}p}{\partial
z^{2}}\right]=\phi c_{t}\frac{\partial p}{\partial t}+\frac{qB}{\Delta
z}\delta(x-x_{w})\delta(y-y_{w})\delta(z-z_{w}),$
$\displaystyle(x,y,z)\in[0,\Delta x]\times[0,\Delta y]\times[0,\Delta z];$
$\displaystyle p(t=0)=p_{i};$ $\displaystyle p(x=0)=p(x=\Delta
x)=p(y=0)=p(y=\Delta y)=p(z=0)=p(z=\Delta z)=p_{i};$

where coordinates $(x_{w},y_{w},z_{w})$ is the point source location and
$\delta(\cdot)$ is the Dirac delta function. Bottomhole pressure dynamics
$p^{\mathrm{s-a}}_{w}(t)$ corresponds to the pore pressure evolution at the
point $x=r_{w}$, $y=0$. The analytical solutions for a uniform-flux point
source and uniform-flux fully-penetrating vertical line-source well can be
found in papers by Ozkan (1988); Ozkan and Raghavan (1991a). We apply several
techniques improving the convergence of the series and shortening the
computation time as described in Ozkan (1988); Ozkan and Raghavan (1991b);
Ozkan (1994).

Verification of the developed semi-analytical reservoir filtration model is
carried out (see Fig. 5a). For this purpose we consider the formation with
uniform permeability of 10 mD and compare the bottomhole pressure behavior as
a function of time computed via the results obtained using numerical
$p^{\mathrm{num}}_{w}(t)$ and semi-analytical $p^{\mathrm{s-a}}_{w}(t)$ models
as implemented into MUFITS simulator and commercial software Kappa Saphir,
respectively. We obtain a good match between Kappa Saphir and in-house semi-
analytical model, while there is a small discrepancy between the results of
simulations conducted in MUFITS simulator and benchmark solution in Kappa
Saphir. The latter one can be attributed to implementation of the equation of
state describing slightly compressible fluid embedded into BLACKOIL module of
MUFITS simulator, where the density dependence on pressure includes both
linear the quadratic terms. In Figure 5b, we show the results of
interpretation of welltest in the reservoir with the permeability distribution
$k(x,y)$ presented in Figure 4a. Here, we compare the time dependencies
$p^{\mathrm{num}}_{w}(t)$ and $p^{\mathrm{s-a}}_{w}(t)$ computed via the
numerical simulator MUFITS and semi-analytical reservoir model, respectively.
The latter curve corresponds to the integral permeability 9.22 mD and skin
factor -0.49 obtained by solving the minimization problem (11).

Figure 5:  Plot (a) shows comparison of bottomhole pressure dynamics obtained
in in-house semi-analytical reservoir model (dashed black curve), numerical
model in MUFITS simulator (solid gray curve) and semi-analytical model
implemented into Kappa Saphir (dashed red curve); in plot (b) we show the
results of solution to minimization problem (11) (welltest interpretation
using semi-analytical model) in a reservoir with the non-uniform absolute
permeability distribution shown in Fig. 4a.

Now as we described the mechanistic modelling workflow to evaluate the
integral absolute rock permeability in the area surrounding a vertical well
and the skin factor, we consider a machine learning algorithm (namely,
artificial neural network or ANN) to develop the surrogate model. In Figure 6
we show its schematic representation.

Figure 6:  Schematic representation of an artificial neural network.

In the framework of current study, ANN solves the regression problem. As an
input, it takes the vector of 14 components, Eq. (8), describing the absolute
permeability field $k(x,y)$ around the well inside the square of 500 m size
and predicts the integral permeability $\widetilde{k}$ and skin factor
$\widetilde{S}$. ANN includes input and output layers as well as several
hidden layers. Each layer consists of nodes, and each node contains a number.
In the first layer, the quantity of nodes is equal to the number of the input
features (14 nodes in our case). In the last layer, the quantity of nodes
corresponds to the number of output features (two parameters in our case as
described above). The quantity of hidden layers and corresponding number of
nodes are hyperparameters, which are found by conducting a series of numerical
experiments. In the current machine learning model, we take 3 hidden layers
with 64, 128, and 64 nodes, respectively. All ANN layers are fully-connected,
so that the node $k$ in the layer $n$ is linked to all nodes in the subsequent
layer $n+1$.

ANN estimates the output parameters via the forward propagation. In this
procedure, the machine learning algorithm computes the values at each node,
e.g., $y_{n+1,k}$, where $n+1$ is the layer number and $k$ in the node number,
according to the following steps: (i) calculation of a linear combination of
values at the nodes in the layer $n$ ($\mathbf{y}_{n}$) with weights
$W_{n\rightarrow n+1,k}$; (ii) application of a non-linear activation function
$g(\cdot)$ to this linear combination. Steps (i) and (ii) can be summarized as
$\mathbf{y}_{n+1}=g(W_{n\rightarrow n+1}\mathbf{y}_{n}+\mathbf{b})$, where
$\mathbf{b}$ denotes the bias vector. In the present model, we utilize ReLu
activation function, $g(x)=\max(0,x)$. Non-linearity is used to propagate the
parameters in between all layers except for the connection between the last
hidden and output layers.

The weights of ANN, which can be represented as components of matrices
$W_{n\rightarrow n+1}$, are tuned via the gradient descent optimization
algorithm minimizing the loss function value and, in the regression task, it
is a mean square error (MSE):

$\mathcal{L}(y,\hat{y})=\frac{1}{n_{\text{s}}\cdot
n_{\text{out}}}\sum_{i=1}^{n_{\text{s}}}\sum_{j=1}^{n_{\text{out}}}(y_{i,j}-\hat{y}_{i,j})^{2},$

where $y_{i}$ and $\hat{y}_{i}$ denote true and predicted output vectors for
sample $i$, respectively; $n_{\text{out}}=2$ is the number of output features;
$n_{\text{s}}$ is the number of samples in a batch (mini-batch gradient
descent optimization is considered) and is determined experimentally, it
equals 32 in the present model. Computation of loss function gradients with
respect to weights is the backward propagation. We utilize Adam optimizer
(Kingma and Ba, 2014) as the gradient descent realization.

We apply ANN implemented in Scikit-learn library (Pedregosa et al., 2011),
namely, function MLPRegressor(). Synthetic dataset consisting of 40000 data
points is divided into training and test sets in the ratio of 4 to 1,
respectively. We tune hyperparameters of ANN including batch size, initial
learning rate and strength of the L2 regularization term based on the training
dataset using the random search algorithm and cross-validation technique.
Random search selects randomly hyperparameter values within the predefined
intervals and estimates the performance of machine learning model at a
specified set of hyperparameters using the cross-validation approach. This
procedure repeats 1000 times leading to the optimal combination of
hyperparameters corresponding to the best model performance. We evaluate the
ANN performance in terms of the mean absolute error (MAE). In the cross-
validation procedure, the training dataset is divided into $M$ equal parts,
$M-1$ partitions are used for training the machine learning algorithm, while
the remaining part is utilized for its validation. The procedure is repeated
$M$ times, and different validation part is taken at each iteration. As a
result, we obtain $M$ validation (MAE) scores and average them, and this
averaged metric is utilized by the random search algorithm to identify the
optimal values of hyperparameters.

The random search algorithm combined with the cross-validation approach
provides the following values of hyperparameters: batch size is 32, initial
learning rate is $2\cdot 10^{-3}$, prefactor before regularization term is
0.6. It is recommended to scale input features before applying ANN, which is
carried out using StandardScaler() function with parameters determined using
the training data. Additionally, early stopping technique is utilized to
prevent ANN from overfitting (MLPRegressor() includes this option and leaves
10% of the training dataset as a validation part to control the predictive
capability of ANN during training phase). When the hyperparameters are tuned,
ANN is trained at the training dataset, and its overall performance in terms
of MAE, MSE and coefficient of determination (R2) is estimated using the test
set.

Figure 7 shows the cross-plot with predictions of the integral permeability
(a) and skin factor (b) using the surrogate model based on ANN machine
learning algorithm. In this chart, the true values of output parameter are
plotted at the $x$-axis, while the predicted ones are at the $y$-axis. Note
that the cloud of points is oriented along the line of ideal prediction $y=x$.
Only few data points are poorly estimated (with an error larger than 15%),
e.g., the samples with skin factor close to the upper and lower limits and in
the vicinity of zero. We compute MAE, MSE and R2 scores using the predicted
and true values of the integral permeability and skin factor for the training
and test datasets separately. The results are summarized in Table 1. We
obtained the acceptable accuracy of the surrogate model based on ANN, so that
it can be used for calculations of the absolute permeability map according to
the methodology outlined in Section 3.

Figure 7:  Cross-plot with the predictions of the integral permeability (a) and skin-factor (b) by ANN; blue points are related to the train set, while the orange ones correspond to the test part.  Table 1: Surrogate model performance at train and test sets for the integral permeability and skin factor in terms of MAE, MSE and R2 metrics. Output parameter | Dataset | MAE | MSE | R2
---|---|---|---|---
integral permeability | train | 0.25 | 0.153 | 0.931
test | 0.256 | 0.153 | 0.925
skin factor | train | 0.08 | 0.016 | 0.909
test | 0.081 | 0.016 | 0.901

##  5 Results and discussion

In this section we present the results of simulations using the proposed
combined mechanistic and machine learning workflow. We apply the technique
outlined in Section 3 to approximate the absolute permeability map of
synthetic reservoir model, namely, “Egg Model” as described by Jansen et al.
(2014). In Fig. 8 we show the permeability map at the top view of the model.

Figure 8: The top view permeability map of reservoir in “Egg Model” (Jansen et
al., 2014); locations of producers (red) and injectors (blue) are shown by
crosses; colour map corresponds to the original permeability distribution;
dashed boxes denote reservoir domains considered in physics-based averaging of
the actual $k_{0}(x,y)$ and approximate $k(x,y)$ absolute permeability
distribution around each well.

The Egg Model consists of an ensemble of 101 three-dimensional absolute
permeability field realizations of a channelized oil formation. Permeability
field is given in the discrete form modelled with $60\times 60\times 7=25.000$
grid cells. The number of active blocks is 18553, and non-active cells are
located around the reservoir so that its shape resembles the form of an egg.
The channels have high permeability values, while the domains between them are
low-permeable. Such structure of the absolute permeability field resembles the
winding river patterns observed in the fluvial systems. In the current
analysis, we utilize a single realization of the absolute permeability field
and take its first layer. Therefore, 3600 original grid cells are considered
and passed to MUFITS simulator so that our hydrodynamic model is effectively
two-dimensional. The maximum and minimum values of the chosen permeability map
are scaled to the interval of $[1,15]$ mD applied for the generation of the
synthetic dataset (Section 4, equation (9)), and the obtained permeability
field is shown in Figure 8.

Let us discuss the geometrical parameters of the formation, fluid and rock
properties, initial and boundary (internal and external) conditions
incorporated into the Egg Model. Majority of the model parameters are similar
to that described in (Jansen et al., 2014). Numerical simulations of the two-
phase (oil-water) filtration is performed using MUFITS reservoir simulator. We
increase the original lateral reservoir size up to 2 km, and its thickness is
set to 10 m (a single layer of mesh cells). As a result, the reservoir
dimensions are 2 km $\times$ 2 km $\times$ 10 m, so that the grid block length
and width are 33.3 m, while the cell height is 10 m. Porosity distribution is
uniform. Oil and water are slightly compressible liquids with fixed
viscosities (BLACKOIL module of MUFITS simulator is applied), while the rock
is incompressible. Relative permeabilities of oil and water are governed by
Corey model, and the capillary pressure is absent. At the initial state, pore
pressure is uniform and equals 400 bar (the formation is located at 4 km
depth), and water saturation is non-zero. The formation is intersected by 12
vertical wells, 4 producers and 8 injectors, and their locations are marked by
red and blue crosses in Figure 8 and shown in Table 3. Water flooding is the
major production mechanism. Vertical wells operate under constant bottomhole
pressure (internal boundary condition), namely, 350 bar for producers and 450
bar for injectors, respectively. The simulation period is 30 years. All
external boundaries of the reservoir are closed. Values of model parameters
are summarized in Table 2.

Table 2: Parameters of the synthetic reservoir Egg Model. Parameter | Value
---|---
Grid-block size | 33.3 m $\times$ 33.3 m $\times$ 10 m
Porosity | 0.2
Oil compressibility | 10-5 bar-1
Water compressibility | 10-5 bar-1
Oil dynamic viscosity | 5 cP
Water dynamic viscosity | 1 cP
End-point relative permeability, oil | 1
End-point relative permeability, water | 1
Corey exponent, oil | 1
Corey exponent, water | 1
Residual-oil saturation | 0.1
Connate-water saturation | 0.1
Initial reservoir pressure | 400 bar
Initial water saturation | 0.1
Production well bottom-hole pressures | 350 bar
Injection well bottom-hole pressures | 700 bar
Well-bore radius | 0.1 m
Simulation time | 30 years

Further, we discuss the generation of actual absolute permeability
distributions $k_{0}(x,y)$ required for the construction the approximate
fields $k(x,y)$, i.e., $k^{\mathrm{WL}},k^{\mathrm{WT}},S$, using the proposed
method described in Section 3. The values
$\left\\{k^{\mathrm{WL}}_{i}\right\\}_{i=1}^{12}$ correspond to the absolute
permeability at the well locations $\left\\{x_{i},y_{i}\right\\}_{i=1}^{12}$,
$k^{\mathrm{WL}}_{i}=k(x_{i},y_{i})$, and are listed in Table 3. The integral
permeability and skin factor are computed using the synthetic well test
procedure described in Section 4. To a well with index $j$, the following
procedures are applied: (i) we cut a square of size
$\mathcal{D}=500~{}\mathrm{m}$ with sides parallel to the coordinate axes
(region $\mathcal{B}_{j}^{\mathcal{D}}$); (ii) pass the absolute permeability
field $k_{0}(x,y),(x,y)\in\mathcal{B}_{j}^{\mathcal{D}}$ into MUFITS simulator
and carry out the calculations of drawdown test (single-phase fluid problem,
production with constant flow rate); (iii) conduct interpretation of the
obtained bottomhole pressure behaviour using the semi-analytical reservoir
model via the solution of the minimization task (11). We apply the same
hydrodynamic model parameters as described in Section 4 during the synthetic
well test procedure except for the flow rate and production period that, which
are set to $q=10$ m3/d and $t\in[0,360]$ d, respectively.

Results of synthetic well test procedure according to steps (i) – (iii)
described above and applied to the production “PROD1” and injection “INJECT1”
wells are shown in Fig. 9. When the absolute permeability field is essentially
heterogeneous (e.g., the domain $\mathcal{B}^{\mathcal{D}}$ around the
producer “PROD2”, which is intersected by the highly permeable channel in
vertical direction), the solution to the minimization problem (11) providing
the excellent match between the numerical bottomhole pressure dynamics and
semi-analytical one does not exist and we found approximate values of
$k^{\mathrm{WT}}$ and $S$. Note that it is impossible to clip a square around
injectors located along the lateral border of the model, and, in this case, we
perform the synthetic well test working with a symmetry element (quarter or
half of the formation as shown in Fig. 8) in MUFITS simulator, example is
shown in Figure 9b. Estimated values of the integral permeability
$k^{\mathrm{WT}}$ and skin factor $S$ are summarized in Table 3.

Figure 9: Results of synthetic well test procedure applied to producer “PROD 1” (a) and injector “INJECT1” (b); in the left column of plot plots the actual permeability map $k_{0}(x,y),(x,y)\in\mathcal{B}^{\mathcal{D}}$ is shown, while the results of the numerical drawdown test and its interpretation are given in the right column of plots; since the majority of injectors are located along the reservoir boundary, the synthetic well test is performed using the symmetry element as shown in plot (b). Table 3: Wells locations $\left\\{x_{i},y_{i}\right\\}_{i=1}^{12}$, permeability values obtained from well logging $\left\\{k^{\mathrm{WL}}_{i}\right\\}_{i=1}^{12}$ and well test $\left\\{k^{\mathrm{WT}}_{i}\right\\}_{i=1}^{12}$ (or integral permeability) measurements as well as skin factor values $\left\\{S_{i}\right\\}_{i=1}^{12}$ in the Egg Model. Index | Well ID | x, m | y, m |  $k^{\mathrm{WL}}$, mD |  $k^{\mathrm{WT}}$, mD | $S$
---|---|---|---|---|---|---
1 | PROD1 | 517 | 1417 | 3.58 | 3.38 | -0.33
2 | PROD2 | 1150 | 1317 | 11.44 | 7.92 | -1.35
3 | PROD3 | 750 | 517 | 2.10 | 2.39 | 0.27
4 | PROD4 | 1383 | 583 | 3.10 | 3.08 | -0.19
5 | INJECT1 | 150 | 1883 | 3.33 | 3.38 | 0.10
6 | INJECT2 | 983 | 1750 | 5.19 | 3.72 | -1.31
7 | INJECT3 | 50 | 1150 | 4.22 | 5.00 | 0.45
8 | INJECT4 | 883 | 950 | 2.00 | 2.22 | 0.17
9 | INJECT5 | 1650 | 1150 | 6.51 | 5.14 | -1.01
10 | INJECT6 | 250 | 283 | 3.41 | 3.22 | -0.24
11 | INJECT7 | 1050 | 50 | 1.83 | 2.22 | 0.95
12 | INJECT8 | 1883 | 183 | 3.68 | 3.43 | -0.30

Using well locations $\left\\{x_{i},y_{i}\right\\}_{i=1}^{12}$, absolute
permeabilities $\left\\{k^{\mathrm{WL}}_{i}\right\\}_{i=1}^{12}$,
$\left\\{k^{\mathrm{WT}}_{i}\right\\}_{i=1}^{12}$, skin factor values
$\left\\{S\right\\}_{i=1}^{12}$ and ANN-based surrogate model described in
Section 4, we solve the optimization problem (5) and estimate the values of
kernel regression parameters
$\Omega=\left\\{\left\\{k^{\mathrm{near}}_{i},k^{\mathrm{far}}_{i}\right\\}\big{|}_{i=1}^{N},r_{d},r_{g},\alpha,\beta,\gamma,\delta\right\\}$.
Note that the solution of the optimization task (5) is not unique, so that
below we show a permeability map $k(x,y)$ approximating the actual
distribution $k_{0}(x,y)$ and providing an acceptable match in terms of the
oil production and water injection profiles obtained during the numerical
simulations using these maps.

High-permeability channels contain producer “PROD2“ and injectors “INJECT3”,
“INJECT7” so that corresponding reservoir domains in the constructed
approximate map have larger permeability (about 8 mD) as compared to the
remaining zones, where the absolute permeability varies in between 3 and 5 mD
(see Fig. 10a). Using the approximate permeability distribution $k(x,y)$, we
can determine the permeability values at the well locations
$\left\\{k(x_{i},y_{i})\right\\}_{i=1}^{12}$, and, using the surrogate model,
we estimate the integral permeability
$\left\\{\widetilde{k}_{i}\right\\}_{i=1}^{12}$ and skin factor
$\left\\{\widetilde{S}_{i}\right\\}_{i=1}^{12}$. The predicted and true values
of the properties $k^{\mathrm{WL}}$, $k^{\mathrm{WT}}$, $S$ are shown in the
cross-plots in Figure 10, plots (b) and (c). Note that the integral
permeability corresponding to reservoir area surrounding injector “INJECT2”
and skin factor for the producer “PROD2” are estimated with large error, while
the remaining values are fitted with acceptable accuracy.

Figure 10:  Plot (a) shows the comparison of actual absolute permeability map
$k_{0}(x,y)$ and its approximation $k(x,y)$ evaluated via equation (3) with
substituted parameters from Table 3 and variables obtained via the solution of
minimization problem (5); cross-plots (b) and (c) show actual values of
$k^{\mathrm{WL}}$, $k^{\mathrm{WT}}$, $S$ and the ones obtained using
approximate absolute permeability map $k(x,y)$.

Next, we carry out numerical reservoir modeling using MUFITS simulator with
the parameters listed in Table 2 combined with actual and approximate absolute
permeability distributions, and we denote these cases as “ACTUAL” and “APPROX”
for brevity. Figures 11(a) and (b) show pore pressure and water saturation
maps at the end of computation period (30 years), while the relative
difference between the compared cases are given in plots (c) and (d). The
approximate absolute permeability map allows to obtain qualitatively similar
distributions of pressure and water saturation upon a fairly long operation
period. Quantitative differences are small in terms of pressure (less than
4%). As expected, relative difference in saturations $\Delta
S_{\mathrm{water}}$ reaches large values in the high-permeable channels (e.g.,
for the domain around the producer “PROD2” it is 80%), since the approximate
permeability distribution can not reproduce them due to simple parametrization
of permeability maps using kernel functions in the form of exponents depending
on the distance to the wells. However, in low-permeability zones, the
differences in the water saturation fields between analyzed solutions are
small.

Figure 11: Plots (a) and (b) show distributions of pressure and water
saturation at the end of simulation period (30 years) computed via the
hydrodynamic simulator MUFITS in which the actual (ACTUAL) and approximate
(APPROX) absolute permeability distributions together with the parameters
outlined in Table 2 are incorporated. The relative differences between the
compared distributions are shown in panels (c) and (d) for pressure $\Delta
p=\left|p_{\mathrm{ACTUAL}}-p_{\mathrm{APPROX}}\right|/p_{\mathrm{ACTUAL}}$
and water saturation $\Delta
S_{\mathrm{water}}=\left|(S_{\mathrm{water}})_{\mathrm{ACTUAL}}-(S_{\mathrm{water}})_{\mathrm{APPROX}}\right|/(S_{\mathrm{water}})_{\mathrm{ACTUAL}}$
fields, respectively.

Now we analyze the temporal dependence of well flow rates. Fig. 12(a) shows
the dynamics of cumulative production of oil and water, as well as total
injected water volume. We obtained a good match between actual permeability
map (solid lines) and approximate permeability map (dashed) cases. We also
demonstrate oil production rate and water injection rate for all producers and
injectors in plots (b) and (c) of Fig. 12. A notable discrepancy between the
results of reservoir simulations using original and approximated permeability
maps is obtained only for injector “INJECT2”.

Figure 12: Plot (a) shows the dynamics of the cumulative oil (red lines) and
water (green lines) production, as well as total injected water volume
(magenta lines); dynamics of the oil production rate of producers “PROD1” –
“PROD4” are shown in plot (b), while plot (c) shows water injection rate
histories of injectors “INJECT1” – “INJECT8”; solid lines refer to the results
of numerical simulations using original permeability distribution (ACTUAL),
while dashed lines corresponds to that carried out using the approximate
permeability map (APPROX).

Based on the obtained results shown in Figs. 11 and 12, we conclude that the
approximate permeability field $k(x,y)$ (see Fig. 10a) is hydrodynamically
similar to the actual distribution $k_{0}(x,y)$ since it provides the
production/injection profiles close to that obtained using the original
permeability map.

##  6 Discussion

In the framework of two-dimensional reservoir model and Nadaraya-Watson (NW)
kernel regression we developed computationally efficient surrogate model for
evaluation of integral permeability of reservoir around vertical wells as well
as skin factor. As shown in the previous section, the developed model after
incorporation into algorithm for construction of global absolute reservoir
permeability map, allows us to obtain permeability distribution, which is
hydrodynamically similar to the original one. We stress that the value of the
obtained result is justified by essentially non-homogeneous permeability
distribution and two-phase filtration in the synthetic reservoir, so that the
considered synthetic case is close to real oilfield conditions.

Nevertheless, there are several limitations of the proposed general workflow
and surrogate model of integral permeability evaluation which we would like to
discuss.

  1. 1.

Relatively simple parameterization of permeability maps in the form of NW
kernel regression based on exponent functions depending on the distance to
wells. This form allows developing computationally efficient permeability cube
construction algorithm, while, as a result, we obtain relatively smooth
distributions, which does not allow resolving abrupt permeability variations
typical of real geological conditions. In particular, failure to reproduce
high-conductivity channels in fluvial geological conditions can result in
significant errors in predicting water breakthrough in oilfields with water
injectors;

  2. 2.

Two-dimensional problem formulation, which does not allow one to resolve
heterogeneity of rock properties in vertical direction and layered reservoir
structure;

  3. 3.

Missing hydraulic fractures in the problem formulation, in particular, in
semi-analytical model of reservoir, as well as horizontal well trajectory
(i.e. multifractured wells); as the semi-analytical reservoir model does not
take into account these completions, the resulting integral permeability can
be determined with a significant error, which in turn will spoil the
hydrodynamic similarity and lead to poor prediction of production/injection
rates in reservoir simulations using the constructed permeability maps.

The issue with simplified spacial parametrization of permeability map can be
overcome in several ways. Most obvious one is to use more sophisticated
functions in kernel regression algorithm, in particular, Fourier series (or
any other series of suitable basis functions) with a sufficient number of
modes to resolve typical scale of heterogeneity in current geological
conditions (e.g., the width of the high-permeability channels in Egg Model
described above). The drawback of this approach is potentially a very large
number of parameters describing permeability field (input features), which can
result in very poor performance of the developed surrogate model of integral
permeability evaluation. More promising direction is to develop a surrogate
model to approximate integral permeability of the area surrounding a well
based on a convolutional neural network (CNN) dealing with the input
permeability maps in form of digital images (pictures). Input images can be
grayscale, while the CNN can be either trained from a scratch or using the
transfer learning approach. This approach can be combined with generative
models of global reservoir permeability cube construction, in particular,
generative adversarial networks. They allow creating a permeability map
consistent with geological realism, which can be supplied by additional field
data (e.g., by results of seismic measurement interpretations).

We believe that the second limitation described above can be overcome if a
proper permeability parametrization in the vertical direction is applied, in
particular, Fourier series with a number of modes sufficient to obtain a
desired accuracy in current geological conditions.

As for the third issue, we plan to develop in-house semi-analytical model to
account for an effect of hydraulic fracture treatment on bottomhole pressure
dynamics during welltest made in vertical fractured or horizontal multistage
wells.

Application of the developed surrogate model approximating the integral
absolute permeability of a certain area surrounding vertical wells is not
limited by the global model of reservoir permeability cube construction as
described in this study. It can be used as a computationally efficient
submodel (“approximator”) evaluating the integral permeability of area
surrounding wells in a wide range of global permeability cube algorithms
provided the proper data flow is organized. For any permeability map generated
by the global algorithm, the developed surrogate model evaluates integral
permeabilities (and other integral parameters, e.g., well skin factor), which
can be used to modify the parameters of global permeability distribution
accordingly.

##  7 Summary and conclusions

In this paper, we proposed a novel method to construct the approximate two-
dimensional absolute permeability map $k(x,y)$. We use the available data
describing the actual absolute permeability distribution $k_{0}(x,y)$, namely,
the estimation of the absolute permeability near the well from the well-
logging ($k^{\mathrm{WL}}$), the integral permeability of the zone around the
well ($k^{\mathrm{WT}}$), and skin factor $S$ evaluated from the
interpretation of the well test measurements. The reservoir permeability map
approximation is based on Nadaraya-Watson kernel regression, and its
parameters are tuned via the solution of the optimization problem, in which we
minimize the differences between $k^{\mathrm{WL}}$, $k^{\mathrm{WT}}$, $S$ and
their predicted values. During the solution of minimization problem, the
integral permeability ($\widetilde{k}$) and skin factor ($\widetilde{S}$)
around each well corresponding to the approximate permeability field are
estimated by the surrogate model. It is based on an artificial neural network
machine learning algorithm trained on the physics-based synthetic dataset
generated with using the numerical hydrodynamic simulator (MUFITS) and in-
house semi-analytical reservoir model. We applied the developed approach for
the synthetic reservoir model (Egg Model). It is a two-phase oil-water
problem, in which oil is produced by 4 vertical wells, while 8 injectors
maintain the average pore pressure. Using the original absolute permeability
distribution, we calculated the values of $k^{\mathrm{WL}}$,
$k^{\mathrm{WT}}$, $S$ and constructed the approximate permeability map. By
running numerical reservoir simulations using original and approximate
absolute permeability maps, we obtained rather close results in terms of the
well flow rates, cumulative injected/produced volumes, pressure and water
saturation distributions at the end of the simulation period. As a result, the
proposed approach allows to generate the permeability distribution, which is
hydrodynamically similar to the actual one and eventually provides the
production and injection profiles with acceptable accuracy.

The proposed surrogate model as a submodel can be incorporated in a wide range
of existing methods of absolute reservoir permeability cube construction for
effective well test data fusion. This can be carried out if the developed
model is treated as an “approximator”, which at the input takes the
permeability map in a certain area surrounding the well and returns its
integral value and well skin factor.

There are several directions for future research related to the developed
models, namely, taking into account geological realism (i.e., complex
distributions of permeability either in the form of figures or using complex
parametrization), generalization of the proposed method to 3D taking into
account heterogeneity of the reservoir permeability in vertical direction and
taking into account a vertical fractured well or multiple fractures in
horizontal wells.

## Acknowledgements

The work was supported by the Analytical center under the RF Government
(subsidy 460 agreement
000000D730321P5Q0002, Grant No. 70-2021-00145 02.11.2021).

## Data availability

The original input file for the reservoir simulator MUFITS containing the Egg
Model is given in the website:
The Egg Model Study. The data files for the EGG Model can be found in the
website: The Egg Model - data files.

## Code availability

Repository name: Data_Fusion_HDM_ML

Program language: Python

Software required: Python libraries – numpy, pandas, scipy, sklearn, math,
pickle, joblib, matplotlib

Program size: 11 MB

The source codes are available for downloading at the link:
https://github.com/evgenii-kanin/Data_Fusion_HDM_ML.

## References

  * Afanasyev (2020) Afanasyev, A., 2020.  Mufits reservoir simulation software.  URL: www. mufits. imec. msu. ru .
  * Bai and Tahmasebi (2020) Bai, T., Tahmasebi, P., 2020\.  Hybrid geological modeling: Combining machine learning and multiple-point statistics.  Computers & geosciences 142, 104519\.
  * Bianco et al. (2007) Bianco, A., Cominelli, A., Dovera, L., Nævdal, G., Valles, B., 2007.  History matching and production forecast uncertainty by means of the ensemble kalman filter: A real field application, in: EUROPEC/EAGE Conference and Exhibition, OnePetro.
  * Coutinho et al. (2010) Coutinho, E.J., Emerick, A.A., Li, G., Reynolds, A.C., 2010\.  Conditioning multilayered geologic models to well-test and production-logging data using the ensemble kalman filter, in: SPE Annual Technical Conference and Exhibition, OnePetro.
  * Evensen and Eikrem (2018) Evensen, G., Eikrem, K.S., 2018\.  Conditioning reservoir models on rate data using ensemble smoothers.  Computational Geosciences 22, 1251–1270.
  * Evensen et al. (2007) Evensen, G., Hove, J., Meisingset, H., Reiso, E., Seim, K.S., Espelid, Ø., 2007\.  Using the enkf for assisted history matching of a north sea reservoir model, in: SPE reservoir simulation symposium, OnePetro.
  * He et al. (2000) He, N., Oliver, D.S., Reynolds, A.C., 2000.  Conditioning stochastic reservoir models to well-test data.  SPE Reservoir Evaluation & Engineering 3, 74–79.
  * Horne (1995) Horne, R.N., 1995.  Modern well test analysis.  Petroway Inc 926.
  * Hovanessian (1961) Hovanessian, S., 1961.  Pressure studies in bounded reservoirs.  Society of Petroleum Engineers Journal 1, 223–228.
  * Jansen et al. (2014) Jansen, J.D., Fonseca, R.M., Kahrobaei, S., Siraj, M., Van Essen, G., Van den Hof, P., 2014\.  The egg model–a geological ensemble for reservoir simulation.  Geoscience Data Journal 1, 192–195.
  * Kingma and Ba (2014) Kingma, D.P., Ba, J., 2014.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980 .
  * Kolesnikov et al. (2010) Kolesnikov, D., Shipenkov, R., Ignatov, I., Kostuchenko, S., Cherkas, E., Tsibizova, E., Babuhina, N., Mezhnova, N., Roschin, A., 2010.  Integration of production, geological and geophysical data while building geological models of oilfields, in: Society of Petroleum Engineers-SPE Russian Oil and Gas Technical Conference and Exhibition 2010, RO and G 10, pp. 1368–1371.
  * Martin (1959) Martin, J.C., 1959.  Simplified equations of flow in gas drive reservoirs and the theoretical foundation of multiphase pressure buildup analyses.  Transactions of the AIME 216, 321–323.
  * Nadaraya (1964) Nadaraya, E.A., 1964.  On estimating regression.  Theory of Probability & Its Applications 9, 141–142.
  * Ozkan (1988) Ozkan, E., 1988.  Performance of horizontal wells.  The university of Tulsa.
  * Ozkan (1994) Ozkan, E., 1994.  New solutions for well-test-analysis problems: part iii-additional algorithms, in: SPE Annual Technical Conference and Exhibition, OnePetro.
  * Ozkan and Raghavan (1991a) Ozkan, E., Raghavan, R., 1991a.  New solutions for well-test-analysis problems: part 1—analytical considerations.  SPE Formation Evaluation 6, 359–368.
  * Ozkan and Raghavan (1991b) Ozkan, E., Raghavan, R., 1991b.  New solutions for well-test-analysis problems: part 2—computational considerations and applications.  SPE Formation Evaluation 6, 369–378.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E., 2011\.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research 12, 2825–2830.
  * Perrine (1956) Perrine, R., 1956.  Analysis of pressure-buildup curves, in: Drilling and Production Practice, OnePetro.
  * Rosenblatt (1958) Rosenblatt, F., 1958.  The perceptron: a probabilistic model for information storage and organization in the brain.  Psychological review 65, 386\.
  * Stehfest (1970) Stehfest, H., 1970.  Algorithm 368: Numerical inversion of laplace transforms [d5].  Communications of the ACM 13, 47–49.
  * Thanh and Sugai (2021) Thanh, H.V., Sugai, Y., 2021\.  Integrated modelling framework for enhancement history matching in fluvial channel sandstone reservoirs.  Upstream Oil and Gas Technology 6, 100027.
  * Titus et al. (2022) Titus, Z., Heaney, C., Jacquemyn, C., Salinas, P., Jackson, M., Pain, C., 2022\.  Conditioning surface-based geological models to well data using artificial neural networks.  Computational Geosciences 26, 779–802.
  * Van Everdingen and Hurst (1949) Van Everdingen, A., Hurst, W., 1949\.  The application of the laplace transformation to flow problems in reservoirs.  Journal of Petroleum Technology 1, 305–324.
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, İ., Feng, Y., Moore, E.W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P., SciPy 1.0 Contributors, 2020.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods 17, 261–272.  doi:10.1038/s41592-019-0686-2.
  * Watson (1964) Watson, G.S., 1964.  Smooth regression analysis.  Sankhyā: The Indian Journal of Statistics, Series A , 359–372.
  * Zakirov et al. (2014) Zakirov, E., Indrupskiy, I., Lubimova, O., Shiriaev, I., 2014\.  Geostatistically-consistent history matching, in: ECMOR XIV-14th European Conference on the Mathematics of Oil Recovery, European Association of Geoscientists & Engineers. pp. 1–13.
  * Zakirov et al. (2016) Zakirov, E., Indrupskiy, I., Shiryaev, I., Lyubimova, O., Anikeev, D., 2016.  Advanced geologically-consistent history matching and uncertainty evaluation, in: ECMOR XV-15th European Conference on the Mathematics of Oil Recovery, European Association of Geoscientists & Engineers. pp. cp–494.
  * Zakirov et al. (2018) Zakirov, E., Shiryaev, I., Indrupskiy, I., Lyubimova, O., Arkhipova, E.Y., Anikeev, D., 2018\.  Geostatistically-consistent history matching of lithofacies and reservoir properties applied to synthetic and real-field cases, in: ECMOR XVI-16th European Conference on the Mathematics of Oil Recovery, European Association of Geoscientists & Engineers. pp. 1–21.
  * Zhang and Oliver (2011) Zhang, Y., Oliver, D.S., 2011\.  History matching using the ensemble kalman filter with multiscale parameterization: A field case study.  SPE Journal 16, 307–317.