Physics-Informed Neural Network for the Transient Diffusivity Equation
  in Reservoir Engineering
Physics-Informed machine learning models have recently emerged with some
interesting and unique features that can be applied to reservoir engineering.
In particular, physics-informed neural networks (PINN) leverage the fact that
neural networks are a type of universal function approximators that can embed
the knowledge of any physical laws that govern a given data-set in the learning
process, and can be described by partial differential equations. The transient
diffusivity equation is a fundamental equation in reservoir engineering and the
general solution to this equation forms the basis for Pressure Transient
Analysis (PTA). The diffusivity equation is derived by combining three physical
principles, the continuity equation, Darcy's equation, and the equation of
state for a slightly compressible liquid. Obtaining general solutions to this
equation is imperative to understand flow regimes in porous media. Analytical
solutions of the transient diffusivity equation are usually hard to obtain due
to the stiff nature of the equation caused by the steep gradients of the
pressure near the well. In this work we apply physics-informed neural networks
to the one and two dimensional diffusivity equation and demonstrate that
decomposing the space domain into very few subdomains can overcome the
stiffness problem of the equation. Additionally, we demonstrate that the
inverse capabilities of PINNs can estimate missing physics such as permeability
and distance from sealing boundary similar to buildup tests without shutting in
the well.
Daniel Badawi, Eduardo Gildin
# Physics-Informed Neural Network for the Transient Diffusivity Equation in Reservoir Engineering

###### Abstract

Physics-Informed machine learning models have recently emerged with some interesting and unique features that can be applied to reservoir engineering. In particular, physics-informed neural networks (PINN) leverage the fact that neural networks are a type of universal function approximators that can embed the knowledge of any physical laws that govern a given data-set in the learning process, and can be described by partial differential equations. The transient diffusivity equation is a fundamental equation in reservoir engineering and the general solution to this equation forms the basis for Pressure Transient Analysis (PTA). The diffusivity equation is derived by combining three physical principles, the continuity equation, Darcy's equation, and the equation of state for a slightly compressible liquid. Obtaining general solutions to this equation is imperative to understand flow regimes in porous media. Analytical solutions of the transient diffusivity equation are usually hard to obtain due to the stiff nature of the equation caused by the steep gradients of the pressure near the well. In this work we apply physics-informed neural networks to the one and two dimensional diffusivity equation and demonstrate that decomposing the space domain into very few subdomains can overcome the stiffness problem of the equation. Additionally, we demonstrate that the inverse capabilities of PINNs can estimate missing physics such as permeability and distance from sealing boundary similar to buildup tests without shutting in the well.

Physics-informed Neural Networks (PINNs) Physics-informed Machine Learning Diffusivity Equation Reservoir Engineering Buildup Tests

## 1 Introduction

The diffusivity equation is a stiff spatiotemporal non-linear partial differential equation (PDE) that describes fluids flow in porous media. It results from combining three physical principles, the material balance equation, Darcy's law and the equation of state of slightly compressible liquid. Thus it is regarded the most important equation in reservoir engineering. Analytical solutions for some simplified forms of the diffusivity equation exist however their practical applications are very limited especially for history matching applications. For single-phase infinite-acting circular reservoirs analytical solution is obtained using Boltzmann's variable transformation which transforms the PDE to an ordinary differential equation (ODE). For circular bounded reservoirs analytical solutions are available in Laplace domain and thus require inversion algorithms (Stehlest, 1970, Murli and Rizzardi, 1990, de Hoog et al., 1982) for the PDE solution to be expressed in the temporal domain.

The stiffness of the PDE which is caused by the steep gradients near the well poses a major challenge even to numerical methods where very fine meshes need to be considered. Numerical methods are very good tools to solve the diffusivity equation, in fact it is the essence of reservoir simulations. However for well testing techniques which are inverse problems, numerical methods are very limited. Another drawback of numerical methods is that they only offer solutions at discrete specified times which limits our ability to calculate the flow rates at all times. For that, we prefer to work with continuous analytical solutions or any good approximation of analytical solutions.
Recently, there is noticeable and growing increase in the use of machine learning (ML) applications in various areas in science and engineering which most of them are fully data-driven. One major shortcoming of data-driven machine learning is the failure to honor any physical relationship between inputs and outputs. Consequently, they suffer to extrapolate predictions outside the training domain and therefore falls short when applied to time-dependent reservoir engineering applications where temporal extrapolation is essential.

A PINN is a scientific machine learning technique used to solve problems involving PDEs (Raissi et al., 2019). PINNs approximate PDE solutions by training a neural network to minimize a loss function; it includes terms reflecting the initial and boundary conditions along the space-time domain's boundary and the PDE residual at selected points in the domain (called collocation points). Incorporating a residual network that encodes the governing physics equations is a significant novelty with PINNs. The basic concept behind PINN training is that it can be thought of as an unsupervised strategy that does not require labelled data, such as results from prior simulations or experiments. The PINN algorithm is essentially a mesh-free technique that finds PDE solutions by converting the problem of directly solving the governing equations to a loss function optimization problem. It works by integrating the mathematical model into the network and reinforcing the loss function with a residual term from the governing equation (Cuomo et al., 2022). Additionally, using neural network back-propagation, when a PINN model is trained, not only the PDE solution is obtained but also the solution derivatives w.r.t inputs are also obtained. Therefore, PINNs can be regarded as alternatives to analytical solutions.

In this paper, our work can be broadly split into two sections. In section 1, we apply the PINN methodology to the 1D cartesian diffusivity equation and show how to utilize the inverse capabilities of PINNs to predict some unknown reservoir properties which previously were estimated using well testing techniques. In section 2, we apply PINNs to the radial diffusivity equation in a circular bounded reservoir under both constant bottom-hole pressure (BHP) and under constant production flow rate. We also present the domain-decomposition method to tackle the stiffness of the radial diffusivity equation.

## 2 Background

### Diffusivity Equation

The diffusivity equation describes the flow in porous media and can be derived using the material balance equation and Darcy's law. Taking an infinitesimal volume element as in Fig. 1, the material balance states that the difference between the mass flow rate entering the volume element and the mass flow rate leaving the volume element must equal the rate of change of mass accumulating in the volume element due expansion of the pore volume and change of mass contributed by sources/sinks. Assuming 2D case and putting in equation form gives the following:

\[\underbrace{\big{(}q\rho\big{|}_{x}+q\rho\big{|}_{y}\big{)}}_{\text{Mass flow rate}}-\underbrace{\big{(}q\rho\big{|}_{x+dx}+q\rho\big{|}_{y+dy}\big{)}}_{\text{OUT}} =\underbrace{\frac{\partial m}{\partial t}}_{\text{Rate of change}}+ \underbrace{q_{ss}\cdot dV}_{\text{sources}\,/\,\text{sinks}} \tag{1}\]

where \(q\rho\big{|}_{x}\) is the flux at \(x\). \(q\) is the volumetric flux in \([m^{3}/sec]\), \(\rho\) is the density of the fluid in \([kg/m^{3}]\) and \(q_{ss}\) is the sources/sinks term in \([kg/m^{3}]\). In this work we don't have sources or sinks and therefore the term will be neglected. \(m\) is the accumulated mass inside the element volume and is given by:

\[m=\rho dV\]

with \(dV\) is the volume element in \([m^{3}]\). It is important to note that the fluid can occupy a fraction of the element volume called the pore volume and is given by:

\[dV_{p} =\phi\cdot dV\] \[dV =dx\cdot dy\cdot dz\]

where \(\phi\) is the porosity of the element volume \(dV\) and defined as the volume that can be occupied by the fluid divided by the total geometric volume. As a result of the definition \(\phi\in[0,1]\). Therefore, the element volume \(dV\) can be divided into two volumes: pore volume \(dV_{p}\) and rock volume \(dV_{r}\).

\[dV=dV_{p}+dV_{r}\]

\[-\big{(}q\rho\big{|}_{x+dx}-q\rho\big{|}_{x}\big{)}-\big{(}q\rho\big{|}_{y+ dy}-q\rho\big{|}_{y}\big{)}=dV\frac{\partial\big{(}\phi\rho\big{)}}{\partial t} \tag{2}\]
\[-\frac{1}{A_{x}}\frac{\left(q\rho\big{|}_{x+dx}-q\rho\big{|}_{x}\right)}{dx}- \frac{1}{A_{y}}\frac{\left(q\rho\big{|}_{y+dy}-q\rho\big{|}_{y}\right)}{dy}= \frac{\partial\big{(}\phi\rho\big{)}}{\partial t} \tag{3}\]

where \(A_{x}=dy\cdot dz\) and \(A_{y}=dx\cdot dz\).

Simplifying eq. 3 gives:

\[-\frac{1}{A_{x}}\frac{\partial\big{(}q\rho\big{)}_{x}}{\partial x}-\frac{1}{A_ {y}}\frac{\partial\big{(}q\rho\big{)}_{y}}{\partial y}=\frac{\partial\big{(} \phi\rho\big{)}}{\partial t} \tag{4}\]

Darcy's law [20] can be mathematically written as the following:

\[q=-\frac{kA}{\mu}(\nabla p+\rho g\nabla z) \tag{5}\]

where \(q\) is the volumetric flux in \([m^{3}/sec]\), \(k\) is the permeability \([m^{2}]\), \(A\) is the flux area in \([m^{2}]\), \(\mu\) is the fluid viscosity in [\(pa\cdot sec\)], \(\nabla p\) is the pressure gradient, \(\rho\) is the fluid density in \([kg/m^{3}]\), \(g\) is the gravity constant in \([m/sec^{2}]\), and \(\nabla z\) is the gradient of the elevation (\(z\) is in the gravity direction). Neglecting the gravity term, the 2D Darcy's law can be expressed as:

\[q_{x}=-\frac{k_{x}A_{x}}{\mu}\frac{\partial p}{\partial x} \tag{6a}\] \[q_{y}=-\frac{k_{y}A_{y}}{\mu}\frac{\partial p}{\partial y} \tag{6b}\]

where \(A_{x}\) and \(A_{y}\) is the flux area perpendicular to \(x\) direction and \(y\) direction respectively. \(k_{x}\) and \(k_{y}\) is the permeability in \(x\) and \(y\) direction respectively.

Substituting eqs. 6a, 6b into eq. 4 and expanding the right hand side gives:

\[\frac{\partial}{\partial x}\bigg{(}\rho\frac{k_{x}}{\mu}\frac{\partial p}{ \partial x}\bigg{)}+\frac{\partial}{\partial y}\bigg{(}\rho\frac{k_{y}}{\mu} \frac{\partial p}{\partial y}\bigg{)}=\phi\frac{\partial\rho}{\partial t}+ \rho\frac{\partial\phi}{\partial t} \tag{7}\]

The time derivative of the density appearing on the right of eq. 7 can be expressed in terms of a time derivative of the pressure by using the basic thermo-dynamic definition of isothermal compressibility of fluids:

\[c_{f}=-\frac{1}{V_{f}}\frac{\partial V_{f}}{\partial p} \tag{8}\]

where \(c_{f}\) is the fluid compressibility in \([1/pa]\), \(V_{f}\) is the fluid volume in \([m^{3}]\) and \(p\) is the pressure acting on the fluid in the pore volume (pore pressure) in \([pa]\), and since

\[\rho=\frac{m}{V_{f}}\]

then the compressibility can be alternatively expressed as:

\[c_{f}=-\frac{\rho}{m}\frac{\partial}{\partial p}\bigg{(}\frac{m}{\rho}\bigg{)} =\frac{1}{\rho}\frac{\partial\rho}{\partial p} \tag{9}\]

and differentiating with respect to time gives:

\[\rho c_{f}\frac{\partial p}{\partial t}=\frac{\partial\rho}{\partial t} \tag{10}\]

The pore volume also changes with pressure, therefore similar to the fluid compressibility the rock compressibility is defined as:

\[c_{r}=\frac{1}{V_{p}}\frac{\partial V_{p}}{\partial p} \tag{11}\]

Notice that eq. 11 has no negative sign, this because the rock compressibility is expressed in terms of pore volume \(V_{p}\) and not rock volume \(V_{r}\), and since

\[V_{p}=\phi V\]
eq. 11 can be expressed as:

\[c_{r}=\frac{1}{\phi}\frac{\partial\phi}{\partial p} \tag{12}\]

and differentiating with respect to time gives:

\[\phi c_{r}\frac{\partial p}{\partial t}=\frac{\partial\phi}{\partial t} \tag{13}\]

finally, substituting eqs. 10 and 13 in eq. 7 reduces to the latter:

\[\frac{\partial}{\partial x}\bigg{(}\rho\frac{k_{x}}{\mu}\frac{\partial p}{ \partial x}\bigg{)}+\frac{\partial}{\partial y}\bigg{(}\rho\frac{k_{y}}{\mu} \frac{\partial p}{\partial y}\bigg{)}=\phi c_{t}\rho\frac{\partial p}{\partial t} \tag{14}\]

[2pt] where \(c_{t}\) is the total compressibility and is given by:

\[c_{t}=c_{f}+c_{r}\]

Eq. 14 is the basic 2D, partial differential equation of any single phase fluid in a porous medium. The equation of course is highly non-linear due to the implicit pressure dependence of the density, viscosity, and compressibility. For that, it is almost hopeless to find simple analytical solutions of the equation. One option to reduce the non-linearity of the equation is to assume constant compressibility, viscosity, and density and this yields to the following:

\[\frac{\partial}{\partial x}\bigg{(}k_{x}\frac{\partial p}{ \partial x}\bigg{)}+\frac{\partial}{\partial y}\bigg{(}k_{y}\frac{\partial p }{\partial y}\bigg{)}=\phi\mu c_{t}\frac{\partial p}{\partial t} \tag{15}\] \[x\in[x_{e1},x_{e2}],\ \ \ y\in[y_{e1},y_{e2}],\ \ \ t\in[0,t_{f}]\]

In this paper we will work with the 1D Cartesian case, thus eq. 15 reduces to:

\[\frac{\partial}{\partial x}\bigg{(}k_{x}\frac{\partial p}{\partial x}\bigg{)} =\phi\mu c_{t}\frac{\partial p}{\partial t} \tag{16}\] \[x\in[x_{e1},x_{e2}],\ \ \ t\in[0,t_{f}]\]

Eq. 15 can be expressed in cylindrical coordinates as follows:

\[\frac{1}{r}\frac{\partial}{\partial r}\bigg{(}k_{r}r\frac{ \partial p}{\partial r}\bigg{)}=\phi\mu c_{t}\frac{\partial p}{\partial t} \tag{17}\] \[r\in[r_{w},r_{e}],\ \ \ t\in[0,t_{f}]\]

where \(k_{r}\) is the permeability in the radial direction in \([m^{2}]\), \(r_{w}\) and \(r_{e}\) are the wellbore radius and the outer reservoir radius in \([m]\) respectively. As mentioned above, eq. 17 is a stiff PDE and that is due to the steep pressure gradients near the well. One reason for the steep gradients is the small wellbore radius \(r_{w}\) which usually ranges from \(3[in]\) to \(5[in]\). Equations 16 17 are two-point boundary value problems that require an initial condition and two boundary conditions.

## 3 Buildup Test

Buildup tests are perhaps the most popular tests in well testing. They are utilized to estimate reservoir properties such as average reservoir permeability, distance from boundary (fault), wellbore storage, skin factor, and many more. A Traditional build up test is performed to a well producing under constant flow rate. The test starts by shutting in the well and recording the bottom-hole pressure change as a function of time. Then the data is plotted on a diagnostic plot. A diagnostic plot is a log-log plot of the pressure change and pressure derivative (vertical axis) versus elapsed time (horizontal axis). It is qualitative plot used to identify flow regimes at different time periods. For each flow regime period we fit an analytical solution and from that we can obtain the reservoir properties. A typical diagnostic plot looks similar to Fig. 2. At early times the pressure buildup is dominated by wellbore storage, this is when \(\Delta p_{wf}\) and \(\Delta p^{\prime}_{wf}\) are equal with slopes equal to unity, next we have a transition period, followed by infinite acting period, this is when \(\Delta p^{\prime}_{wf}\) is constant, and finally boundary dominated period this is when the boundaries felt the shut-in.
As mentioned previously, under infinite-acting conditions and constant permeability \(k_{r}\), eq. 17 has analytical solution in terms of exponential integral

\[p_{wf}=p_{i}-\frac{q\mu}{2\pi k_{r}h}ei\bigg{(}\frac{\phi\mu c_{t}r_{w}^{2}}{4k_{ r}t}\bigg{)} \tag{18}\]

and since

\[ei(x)\approx-ln(\gamma x)\hskip 28.452756pt\text{for }x<0.01\]

where \(\gamma=1.781\). Therefore, eq. 18 can be approximated as

\[p_{wf}=p_{i}-\frac{q\mu}{2\pi k_{r}h}ln\bigg{(}\frac{4k_{r}t}{\gamma\phi\mu c_ {t}r_{w}^{2}}\bigg{)} \tag{19}\]

\[\Delta p_{wf}\equiv p_{i}-p_{wf}=\frac{q\mu}{2\pi k_{r}h}ln\bigg{(}\frac{4k_{r }t}{\gamma\phi\mu c_{t}r_{w}^{2}}\bigg{)} \tag{20}\]

Figure 1: Infinitesimal Cartesian volume element for the diffusivity equation analysis.

Figure 2: Diagnostic plot. The blue points are the bottom-hole pressure recording after well shut-in. The red points are the derivative of the recorded bottom-hole pressure with respect to shut-in time.
where \(p_{i}\) is the reservoir initial pressure, \(p_{wf}\) is the bottom-hole pressure, \(q\) is the flow rate prior to shut-in. \(\gamma=1.781\), \(t\) is time, and \(h\) is the producing layer thickness in \([m]\).

Diagnostic plots are qualitative plots used to identify different flow periods. For quantitative analysis specialized plots are used. Each flow period has its own specialized plot, for example, the specialized plot associated with the infinite-acting period is a semi-log plot with pressure change and its derivative on the linear vertical axis and elapsed time on the log horizontal axis. The reason for this is simply because it is evident from eq. 19 that the pressure change and elapsed time have a linear relationship on the semi-log plot. Therefore if we plot the pressure data on a semi-log specialized plot we can quantify the value \(m\) which is the slope and subsequently we can estimate \(k_{r}\) using eq. 22.

\[m =\frac{\partial\Delta p_{wf}}{\partial ln(t)}=\frac{q\mu}{2\pi k_ {r}h} \tag{21}\] \[k_{r} =\frac{q\mu}{2\pi mh} \tag{22}\]

### Initial and Boundary Conditions

This work is divided into two sections. Section 1 is associated with eq. 16 and Section 2 is associated with eq. 17 and both of these equations are two-point boundary problems thus, need an initial condition and two boundary conditions. Each section has two cases in which they differ by their boundary conditions. For all cases of both sections the initial condition is defined as:

\[p(r,t=0)=p_{0}=25.0\ [MPa] \tag{2.20}\]

**Section 1:**

* **Case 1.1**, the well is controlled by constant BHP at \(x_{w}=0\), and a no-flow boundary at \(x_{e}\), thus: \[p(x_{w},t)=p_{wf}=3.0\ [MPa]\] \[x_{e}\frac{\partial p}{\partial x}\Big{|}_{x_{e}}=0\]
* **Case 1.2**, the well is placed inside the domain \(x_{w}\in[-500,750]\) at \(x_{w}=0\), and controlled by constant BHP. Also, two no-flow boundaries are set at \(x_{e1}=-500\) and \(x_{e2}=750[m]\), thus: \[p(x_{w}=0,t)=p_{wf}=3.0\ [MPa]\] \[x_{e1}\cdot\frac{\partial p}{\partial x}\Big{|}_{x_{e1}}=x_{e2}\cdot\frac{ \partial p}{\partial x}\Big{|}_{x_{e2}}=0\]

**Section 2:**

* For both cases of section II, the outer boundary is a no-flow boundary, that is: \[q_{e}=0\]

Figure 3: Schematic representation of cases 1.1 and 1.2. The green circle is the well and the brick walls are the no-flow boundaries.
Using Darcy's law: \[q(r,t)=\frac{kA}{\mu}\frac{\partial p}{\partial r}\] \(A\) is the flux area. The flux area at the boundary is given by: \[A_{e}=2\pi r_{e}h\] \(h\) is the producing layer thickness. Combining the above yields to: \[r_{e}\frac{\partial p}{\partial r}\Big{|}_{r_{e}}=0\]
* **Case 2.1**, the well is producing at a constant BHP: \[p(r_{w},t)=p_{wf}=3.0\ [MPa]\]
* **Case 2.2**, the well is producing at a constant rate: \[r_{w}\frac{\partial p}{\partial r}\Big{|}_{r_{w}}=\frac{q_{w}\mu}{2\pi kh}=p_{ch}\]

### Physics-Informed Neural Networks

PINNs can be used to approximate solutions of PDEs. For a general PDE of the form:

\[\frac{\partial u(\mathbf{x},t)}{\partial t}+\mathcal{N}_{x}(u(\mathbf{x},t))=0,\ \ \ \ \ \mathbf{x}\in\Omega,t\in[0,T] \tag{23a}\] \[u(\mathbf{x},t)=g(\mathbf{x},t),\ \ \ \mathbf{x}\in\partial\Omega,t\in[0,T]\] (23b) \[u(\mathbf{x},0)=u_{0}(\mathbf{x}),\ \ \ \mathbf{x}\in\Omega \tag{23c}\]

where \(u(\mathbf{x},t)\) is the solution to eq. 23a, \(\Omega\) is the space domain, \(x\) is a spatial vector variable, \(t\) is time, and \(\mathcal{N}_{x}(.)\) is a differential operator. Equations 23b, 23c are the boundary and initial conditions respectively. Physics-informed neural network \(\tilde{u}(\mathbf{x},t,\mathbf{w})\) can approximate the solution of eq. 23a, where \(\mathbf{w}\) are the weights and biases of the PINN ([Raissi et al., 2019]). The PINN is trained by minimizing the loss function which is given as follows:

\[\mathcal{L}(\mathbf{w})=\lambda_{r}\mathcal{L}_{r}(\mathbf{w})+\lambda_{b}\mathcal{L} _{b}(\mathbf{w})+\lambda_{0}\mathcal{L}_{0}(\mathbf{w}) \tag{24}\]

where \(\lambda_{r},\lambda_{b},\lambda_{0}\in\mathbb{R}\) are weights and \(\mathcal{L}_{r}(w),\mathcal{L}_{b}(w),\mathcal{L}_{0}(w)\) correspond to the residual, boundary conditions, and initial conditions accordingly, and are given by:

\[\mathcal{L}_{r}(\mathbf{w})= \frac{1}{N_{r}}\sum_{i=1}^{N_{r}}r(\mathbf{x}_{r}^{i},t_{r}^{i},\mathbf{ w}), \tag{25a}\] \[\mathcal{L}_{b}(\mathbf{w})= \frac{1}{N_{b}}\sum_{i=1}^{N_{b}}\big{(}\tilde{u}(\mathbf{x}_{b}^{i}, t_{b}^{i},\mathbf{w})-g(\mathbf{x}_{b}^{i},t_{b}^{i})\big{)}^{2},\] (25b) \[\mathcal{L}_{0}(\mathbf{w})= \frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\big{(}\tilde{u}(\mathbf{x}_{0}^{i}, 0^{i},\mathbf{w})-u_{0}(\mathbf{x}_{0}^{i})\big{)}^{2}, \tag{25c}\]

Figure 4: Schematic of the 2D cases of Section II.
where \(\{x_{r}^{i},t_{r}^{i}\}_{i=1}^{N_{r}}\) is a collocation points list sampled from the spatial and temporal domain 23a, \(\{x_{b}^{i},t_{b}^{i}\}_{i=1}^{N_{b}}\) is a boundary points list sampled from the boundaries 23b, and \(\{x_{0}^{i},0^{i}\}_{i=1}^{N_{0}}\) is an initial condition points list sampled from the initial domain 23c. These points is where the loss function components \(\mathcal{L}_{r}(w),\mathcal{L}_{b}(w),\mathcal{L}_{0}(w)\) are minimized respectively. \(r(\mathbf{x},t;\mathbf{w})\) is the residual of the PDE and is defined as:

\[r(\mathbf{x},t;\mathbf{w})=\frac{\partial\tilde{u}(\mathbf{x},t;\mathbf{w})}{ \partial t}+\mathcal{N}_{x}(\tilde{u}(\mathbf{x},t,\mathbf{w})) \tag{26}\]

The training of the PINN is achieved by gradient descent, where the gradients are calculated via automatic differentiation [Atilim Gunes Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, Jeffrey Mark Siskind, 2018]. The gradient descent minimization algorithm is expressed by eq.27 where \(\alpha\) is the learning rate. The architecture of a physics-informed neural network can be schematically visualized as in Fig. 5.

\[\mathbf{w}_{j+1}=\mathbf{w}_{j}-\alpha\frac{\partial\mathcal{L}}{\partial \mathbf{w}_{j}} \tag{27}\]

## 4 Physics-Informed Neural Network Setup

### Equation Scaling and Domain Generation

In this section we show how to reform the diffusivity equation to fit the requirements of PINNs. In general, for efficient training of neural networks the inputs and outputs are scaled between \([0,1]\) or \([-1,1]\). In our study both inputs (space-time) and the output (pressure) are scaled between \([0,1]\) unless stated otherwise.

First, we write eq. 16 and eq. 17 in a dimensionless form. Note that in this work the permeabilities \(k_{x}\) and \(k_{r}\) in eqs. 16,17 are assumed to be constant and therefore can be moved to the right hand side of the equation by division.

Figure 5: Schematic of a physics-informed neural network.
\[\frac{\partial p_{D}}{\partial t_{D}}=\frac{\partial^{2}p_{D}}{ \partial x_{D}^{2}}\] (28a) \[x_{D}=\frac{x}{x_{e}};\ \ \ t_{D}=\frac{k_{x}t}{\phi\mu c_{t}x_{e}^{2}}; \ \ \ p_{D}=\frac{p-p_{wf}}{p_{0}-p_{wf}}\] (28b) \[\frac{\partial p_{D}}{\partial t_{D}}=\frac{1}{r_{D}}\frac{ \partial p_{D}}{\partial r_{D}}+\frac{\partial^{2}p_{D}}{\partial r_{D}^{2}}\] \[r_{D}=\frac{r}{r_{e}};\ \ \ t_{D}=\frac{k_{r}t}{\phi\mu c_{t}r_{e}^{2}}\] \[p_{D}=\frac{p-p_{wf}}{p_{0}-p_{wf}}\ \ \ \ \ \ \text{for case 2.1}\] \[p_{D}=\frac{p_{0}-p}{p_{ch}}\
### Loss Components

For case 1.1 the components of the loss function are as follows:

\[\mathcal{L}_{r}(\mathbf{w})= \frac{1}{N_{r}}\sum_{i=1}^{N_{r}}r(\mathbf{x}_{r,}^{i}t_{r,}^{i},\mathbf{w}), \tag{31}\] \[\mathcal{L}_{w}(\mathbf{w})= \frac{1}{N_{w}}\sum_{i=1}^{N_{w}}\big{(}\tilde{u}(\mathbf{x}_{w}^{i},t _{w}^{i},\mathbf{w})-g(\mathbf{x}_{w}^{i},t_{w}^{i})\big{)}^{2},\] (32) \[\mathcal{L}_{b}(\mathbf{w})= \frac{1}{N_{b}}\sum_{i=1}^{N_{b}}\big{(}\tilde{u}(\mathbf{x}_{b}^{i},t _{b}^{i},\mathbf{w})-g(\mathbf{x}_{b}^{i},t_{b}^{i})\big{)}^{2},\] (33) \[\mathcal{L}_{0}(\mathbf{w})= \frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\big{(}\tilde{u}(\mathbf{x}_{0}^{i}, 0^{i},\mathbf{w})-u_{0}(\mathbf{x}_{0}^{i})\big{)}^{2}, \tag{34}\]

where \(\mathcal{L}_{r}\), \(\mathcal{L}_{w}\), \(\mathcal{L}_{b1}\), \(\mathcal{L}_{b2}\), and \(\mathcal{L}_{0}\) are the loss function components associated with PDE residual, well, right no-flow boundary, left no-flow boundary, and initial condition respectively.

Figure 6: Domain points: (a) Input points for 1D diffusivity PINN model with a single no-flow boundary. (b) Input points for 1D diffusivity PINN model with two no-flow boundary conditions. (c) Input points for the radial diffusivity PINN model with a single no-flow boundary. Red sets are for the PDE (collocation points), green sets are for the well condition, black sets are for the no-flow boundary conditions, and blue sets are for the initial condition
Finally for case 2.1 and 2.2 the components of the loss function are as follows:

\[\mathcal{L}_{r}(\mathbf{w})= \frac{1}{N_{r}}\sum_{i=1}^{N_{r}}r(\mathbf{x}_{r}^{i},t_{r}^{i},\mathbf{w}), \tag{35}\] \[\mathcal{L}_{w}(\mathbf{w})= \frac{1}{N_{w}}\sum_{i=1}^{N_{w}}\big{(}\tilde{u}(\mathbf{x}_{w}^{i},t _{w}^{i},\mathbf{w})-g(\mathbf{x}_{w}^{i},t_{w}^{i})\big{)}^{2},\] (36) \[\mathcal{L}_{b1}(\mathbf{w})= \frac{1}{N_{b1}}\sum_{i=1}^{N_{b1}}\big{(}\tilde{u}(\mathbf{x}_{b1}^{i },t_{b1}^{i},\mathbf{w})-g(\mathbf{x}_{b1}^{i},t_{b1}^{i})\big{)}^{2},\] (37) \[\mathcal{L}_{b2}(\mathbf{w})= \frac{1}{N_{b2}}\sum_{i=1}^{N_{b2}}\big{(}\tilde{u}(\mathbf{x}_{b2}^{i },t_{b2}^{i},\mathbf{w})-g(\mathbf{x}_{b2}^{i},t_{b2}^{i})\big{)}^{2},\] (38) \[\mathcal{L}_{0}(\mathbf{w})= \frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\big{(}\tilde{u}(\mathbf{x}_{0}^{i },0^{i},\mathbf{w})-u_{0}(\mathbf{x}_{0}^{i})\big{)}^{2}, \tag{39}\]

where \(\mathcal{L}_{r}\), \(\mathcal{L}_{w}\), \(\mathcal{L}_{b1}\), \(\mathcal{L}_{b2}\), and \(\mathcal{L}_{0}\) are the loss function components associated with PDE residual, well, right no-flow boundary, left no-flow boundary, and initial condition respectively.

Finally for case 2.1 and 2.2 the components of the loss function are as follows:

\[\mathcal{L}_{r}(\mathbf{w})= \frac{1}{N_{r}}\sum_{i=1}^{N_{r}}r(\mathbf{r}_{r}^{i},t_{r}^{i},\mathbf{w }), \tag{40}\] \[\mathcal{L}_{w}(\mathbf{w})= \frac{1}{N_{w}}\sum_{i=1}^{N_{w}}\big{(}\tilde{u}(\mathbf{r}_{w}^{i}, t_{w}^{i},\mathbf{w})-g(\mathbf{r}_{w}^{i},t_{w}^{i})\big{)}^{2},\] (41) \[\mathcal{L}_{b}(\mathbf{w})= \frac{1}{N_{b}}\sum_{i=1}^{N_{b}}\big{(}\tilde{u}(\mathbf{r}_{b}^{i}, t_{b}^{i},\mathbf{w})-g(\mathbf{r}_{b}^{i},t_{b}^{i})\big{)}^{2},\] (42) \[\mathcal{L}_{0}(\mathbf{w})= \frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\big{(}\tilde{u}(\mathbf{r}_{0}^{i}, 0^{i},\mathbf{w})-u_{0}(\mathbf{r}_{0}^{i})\big{)}^{2}, \tag{43}\]

where \(\mathcal{L}_{r}\), \(\mathcal{L}_{w}\), \(\mathcal{L}_{b}\), and \(\mathcal{L}_{0}\) are the loss function components associated with PDE residual, well, no-flow boundary, and initial condition respectively.

## 5 Results

### 1D Forward Problems:

* **Case 1.1:** In case 1.1, a well producing at constant BHP and a no-flow boundary. The results are shown in Fig. 6(a). Once the PINN is trained production rates can be easily calculated using Darcy's law: \[\tilde{q}=\frac{kA}{\mu}\frac{d\tilde{p}}{dx}\] where \(A\) is flow surface area, and \(\frac{d\tilde{p}}{dx}\) is the pressure gradient which can be obtained using back-propagation. Note that this is a one-dimensional problem, therefore, the flow area \(A\) is set to unity. In Fig. 7 notice that PINNs predictions for both pressure and flow rates are continuous in time whereas the high-fidelity are discrete and user-specified. In other words, the user must specify at which specific times they wish to obtain solutions for pressures and flow rates.
* **Case 1.2:** In case 1.2 a well producing at constant BHP between two no-flow boundaries. The results are shown in Fig. 8.

### Dimensionless Variables

It is important to note that working with scaled dimensionless form of the PDE enables making instantaneous predictions for all reservoir properties. Any change in reservoir properties \((\phi,\mu,k,c_{t},x_{e})\) is equivalent to change in the dimensionless time \(t_{D}\). Unlike numerical methods where any change in the reservoir properties requires restarting the method from the beginning, in PINNs there is no need to retrain the model since predicting for a different reservoir property is equivalent to predicting for a different dimensionless time.

For practical purposes we can generate interactive-plots with sliders for instantaneous predictions as in Fig. 9. Additionally, if we have production rates data placed on the right plot of Fig. 9, we can change reservoir properties until achieving a good match and with that we can estimate reservoir properties associated with the production rates data.

Figure 8: PINN pressure prediction (red dashed) versus high-fidelity numerical solution (solid blue). The relative L2 error is less than 0.1%

Figure 7: (a) - PINN prediction for pressure (red dashed) versus high-fidelity numerical solution (solid blue). (b) - flow rates prediction using PINN (red dashed) versus high-fidelity numerical solution (blue dots). The relative L2 error is less than 0.1%.
### 1D Inverse Problems:

In this section we test the inverse capabilities of PINNs. The problem setup is the following. We introduce a production flow rate history data \(q_{data}\) and we treat one of the reservoir properties as an unknown. We tackle two problems in this section, the first, we treat the permeability \(k_{x}\) as an unknown and the second we treat the distance from the boundary \(x_{e}\) as unknown. For the first problem, eq. 28a now has two unknowns, the pressure and the permeability with one additional boundary condition is introduced that is the flow rates history data shown in Fig. 10. The loss function components now are as follows:

\[\mathcal{L}_{r}(\mathbf{w})= \frac{1}{N_{r}}\sum_{i=1}^{N_{r}}r(\mathbf{x}_{r}^{i},t_{r}^{i},\mathbf{w }), \tag{44}\] \[\mathcal{L}_{w}(\mathbf{w})= \frac{1}{N_{w}}\sum_{i=1}^{N_{w}}\big{(}\tilde{u}(\mathbf{x}_{w}^{i},t _{w}^{i},\mathbf{w})-g(\mathbf{x}_{w}^{i},t_{w}^{i})\big{)}^{2},\] (45) \[\mathcal{L}_{b}(\mathbf{w})= \frac{1}{N_{b}}\sum_{i=1}^{N_{b}}\big{(}\tilde{u}(\mathbf{x}_{b}^{i},t _{b}^{i},\mathbf{w})-g(\mathbf{x}_{b}^{i},t_{b}^{i})\big{)}^{2},\] (46) \[\mathcal{L}_{0}(\mathbf{w})= \frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\big{(}\tilde{u}(\mathbf{x}_{0}^{i},0 ^{i},\mathbf{w})-u_{0}(\mathbf{x}_{0}^{i})\big{)}^{2},\] (47) \[\mathcal{L}_{d}(\mathbf{w})= \frac{1}{N_{d}}\sum_{i=1}^{N_{d}}\big{(}\tilde{q}(0^{i},t_{d}^{i },\mathbf{w})-q_{d}(t_{d}^{i})\big{)}^{2}, \tag{48}\]

Where \(\tilde{q}\) in eq. 48 is the predicted flow rates by the PINN and \(q_{d}\) is the production flow rates data. As shown in Fig. 11, a PINN model was able to learn the unknown permeability that satisfies the PDE, the initial condition, boundary conditions and the production history data with a relative error of 0.16%. Similarly, for the second problem with unknown boundary distance, a PINN model is trained to satisfy the PDE, the initial condition, boundary conditions and the production history data with a relative error of 0.79% as shown in Fig. 12. It is important to note that both models were trained with the noisy data (blue curve) in Fig. 10.

In order for the model to learn an unknown reservoir property, the property is learnt by gradient descent similar to network weights and biases. For the case with unknown permeability we rewrite eq. 16:

\[\begin{split}&\frac{\phi\mu c_{t}}{k_{x}}\bigg{(}\frac{k_{m}}{k_{ m}}\bigg{)}\frac{\partial p}{\partial t}=\frac{\partial^{2}p}{\partial x^{2}}\\ & x\in[0,x_{e}],\ \ \ t\in[0,t_{f}]\end{split} \tag{49}\]

Figure 9: Interactive plot for permeability \(k\) and boundary distance from well \(x_{e}\). As we change the values with the sliders the plots are instantaneously updated with the new values.
where \(k_{m}\) is an upper estimation of the permeability field and \(k_{x}\) is the unknown permeability. Here we take \(k_{m}=300[mD]\).

\[\begin{split}&\frac{\phi\mu c_{t}}{k_{m}}\bigg{(}\frac{k_{m}}{k_{x}} \bigg{)}\frac{\partial p}{\partial t}=\frac{\partial^{2}p}{\partial x^{2}}\\ &\eta=\frac{k_{x}}{k_{m}}\\ &\frac{\phi\mu c_{t}}{k_{m}}\frac{\partial p}{\partial t}=\eta \frac{\partial^{2}p}{\partial x^{2}}\end{split} \tag{50}\]

where \(\eta\) is called the permeability factor. Now we can rewrite eq. 28a as follows:

\[\begin{split}&\frac{1}{t_{D,max}}\frac{\partial p_{D}}{\partial t _{D}}=\eta\frac{\partial^{2}p_{D}}{\partial x_{D}^{2}}\\ & x_{D}=\frac{x}{x_{e}};\ \ \ t_{D}=\frac{k_{m}t}{\phi\mu c_{t}x_{e} ^{2}};\ \ \ p_{D}=\frac{p-p_{wf}}{p_{0}-p_{wf}}\end{split} \tag{51}\]

In addition to the previous, notice now that the residual of the PDE of eq. 51 is also a function of \(\eta\), and the loss function of the residual can be written as follows:

\[\begin{split}\mathcal{L}_{r}(\mathbf{w},\eta)=&\frac{1 }{N_{r}}\sum_{i=1}^{N_{r}}r(\mathbf{x}_{r}^{i},t_{r}^{i},\mathbf{w},\eta)\end{split} \tag{52}\]

The production flow rates history data \(q_{data}\) is also a function of \(\eta\) for the following reason:

\[\begin{split}& q_{d}=\frac{kA}{\mu}\frac{dp}{dx}=\frac{kA}{\mu} \frac{k_{m}}{k_{m}}\frac{dp}{dx}=\frac{k_{m}A}{\mu}\eta\frac{dp}{dx}=\eta\frac {k_{m}A}{\mu}\bigg{(}\frac{p_{0}-p_{wf}}{x_{e}}\bigg{)}\frac{dp_{D}}{dx_{D}} \\ &\underbrace{\frac{q_{d}\mu}{k_{m}A}\bigg{(}\frac{x_{e}}{p_{0}-p _{wf}}\bigg{)}}_{\text{known quantity}}\ \ \ =\ \ \ \ \underbrace{\eta}_{\text{learnt quantity}}\ \ \ \ \times\ \ \ \ \underbrace{\frac{dp_{D}}{dx_{D}}}_{\text{back propagation}}\end{split} \tag{53}\]

\[\begin{split}\mathcal{L}_{d}(\mathbf{w},\eta)=\frac{q_{d}\mu}{k_{m}A} \bigg{(}\frac{x_{e}}{p_{0}-p_{wf}}\bigg{)}\ -\ \eta\ \times\ \frac{dp_{D}}{dx_{D}}\end{split} \tag{54}\]

The total loss function of the network is therefore a function of \(\eta\) and given by:

\[\begin{split}\mathcal{L}(\mathbf{w},\eta)=\mathcal{L}_{r}(\mathbf{w}, \eta)+\mathcal{L}_{0}(\mathbf{w})+\mathcal{L}_{w}(\mathbf{w})+\mathcal{L}_{b}(\mathbf{w})+ \mathcal{L}_{d}(\mathbf{w},\eta)\end{split} \tag{55}\]

Similar analogy can be applied for the case with unknown boundary distance. The results for both cases are shown in Fig. 11 and Fig. 12.
To summarize this section, we showed that from production data we are able to estimate the missing permeability and distance from sealing boundary without having to shut-in the well, whereas, in buildup tests to estimate these properties a well shut-in is required.

### Forward 2D case

To this end, we dealt with the 1D diffusivity equation which is not a stiff PDE, and we saw how a single physics-informed neural network was able to learn the solution. When working with the 2D diffusivity equation which is a very stiff

Figure 11: Permeability prediction for a given production flow rates history data.

Figure 12: Boundary distance prediction for a given production flow rates history data.

Figure 10: Production flow rates history
PDE (especially close to the well) the problem become much more challenging. For several traditional neural network architectures that we tried, a single neural network was not able to learn the solution for the 2D case. A possible reason for this is because the solution of the equation exhibits highly logarithmic nature close to the well, and linear nature very far from the well with a transition nature in between. Applying the conservative PINNs (cPINNs) methodology [20] to a circular reservoir with well in the center, we decompose the reservoir into multiple subdomains and each subdomain is associated with a completely independent neural network, and at the interface between the subdomains we require the pressures predictions and pressure gradients predictions (fluxes) to be equal.

As an example let's take Fig. 13. For domain No.1, the loss function components are:

\[\mathcal{L}_{r}^{(1)}\big{(}\mathbf{w_{1}}\big{)}= \frac{1}{N_{r}}\sum_{i=1}^{N_{r}}r^{(1)}(\mathbf{x}_{r}^{i},t_{r}^{i}, \mathbf{w_{1}})^{2},\] \[\mathcal{L}_{0}^{(1)}\big{(}\mathbf{w_{1}}\big{)}= \frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\bigg{(}\tilde{u}^{(1)}(\mathbf{x}_{0 }^{i},0^{i},\mathbf{w_{1}})-u_{0}(\mathbf{x}_{0}^{i})\bigg{)}^{2},\] \[\mathcal{L}_{w}^{(1)}\big{(}\mathbf{w_{1}}\big{)}= \frac{1}{N_{w}}\sum_{i=1}^{N_{w}}\bigg{(}\tilde{u}^{(1)}(\mathbf{x}_{ w}^{i},t_{w}^{i},\mathbf{w_{1}})-g(\mathbf{x}_{w}^{i},t_{w}^{i})\bigg{)}^{2},\] \[\mathcal{L}_{pc,12}^{(1)}\big{(}\mathbf{w_{1}},\mathbf{w_{2}}\big{)}= \frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\bigg{(}\tilde{u}^{(1)}\big{(}\bm {x}_{c,1}^{i},t_{c,1}^{i},\mathbf{w_{1}}\big{)}-\tilde{u}^{(2)}\big{(}\mathbf{x}_{c,2} ^{i},t_{c,2}^{i},\mathbf{w_{2}}\big{)}\bigg{)}^{2},\] \[\mathcal{L}_{fc,12}^{(1)}\big{(}\mathbf{w_{1}},\mathbf{w_{2}}\big{)}= \mathcal{L}_{r}^{(1)}+\mathcal{L}_{0}^{(1)}+\mathcal{L}_{w}^{(1)}+ \mathcal{L}_{pc,12}^{(1)}+\mathcal{L}_{fc,12}^{(1)}\]

Domain No.2, the loss function components are:

\[\mathcal{L}_{r}^{(2)}\big{(}\mathbf{w_{2}}\big{)}= \frac{1}{N_{r}}\sum_{i=1}^{N_{r}}r^{(2)}(\mathbf{x}_{r}^{i},t_{r}^{i}, \mathbf{w_{2}})^{2},\] \[\mathcal{L}_{0}^{(2)}\big{(}\mathbf{w_{2}}\big{)}= \frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\bigg{(}\tilde{u}^{(2)}(\mathbf{x}_{0 }^{i},0^{i},\mathbf{w_{2}})-u_{0}(\mathbf{x}_{0}^{i})\bigg{)}^{2},\] \[\mathcal{L}_{pc,21}^{(2)}\big{(}\mathbf{w_{1}},\mathbf{w_{2}}\big{)}= \frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\bigg{(}\tilde{u}^{(1)}\big{(}\bm {x}_{c,1}^{i},t_{c,1}^{i},\mathbf{w_{1}}\big{)}-\tilde{u}^{(2)}\big{(}\mathbf{x}_{c,2} ^{i},t_{c,2}^{i},\mathbf{w_{2}}\big{)}\bigg{)}^{2},\] \[\mathcal{L}_{fc,21}^{(2)}\big{(}\mathbf{w_{1}},\mathbf{w_{2}}\big{)}= \frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\bigg{(}\tilde{u}_{x}^{(1)}\big{(} \mathbf{x}_{c,1}^{i},t_{c,1}^{i},\mathbf{w_{1}}\big{)}-\tilde{u}_{x}^{(2)}\big{(}\mathbf{x}_ {c,2}^{i},t_{c,2}^{i},\mathbf{w_{2}}\big{)}\bigg{)}^{2},\] \[\mathcal{L}_{pc,23}^{(2)}\big{(}\mathbf{w_{2}},\mathbf{w_{3}}\big{)}= \frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\bigg{(}\tilde{u}_{x}^{(2)}\big{(} \mathbf{x}_{c,2}^{i},t_{c,2}^{i},\mathbf{w_{2}}\big{)}-\tilde{u}^{(3)}\big{(}\mathbf{x}_{c, 3}^{i},t_{c,3}^{i},\mathbf{w_{3}}\big{)}\bigg{)}^{2},\] \[\mathcal{L}^{(1)}\big{(}\mathbf{w_{1}},\mathbf{w_{2}},\mathbf{w_{3}}\big{)}= \mathcal{L}_{r}^{(2)}+\mathcal{L}_{0}^{(2)}+\mathcal{L}_{pc,21}^{(2 )}+\mathcal{L}_{fc,21}^{(2)}+\mathcal{L}_{pc,23}^{(2)}+\mathcal{L}_{fc,23}^{(2)}\]

Domain No.3, the loss function components are:
\[\mathcal{L}_{r}^{(3)}(\mathbf{w_{3}})= \frac{1}{N_{r}}\sum_{i=1}^{N_{r}}r^{(3)}(\mathbf{x}_{r}^{i},t_{r}^{i}, \mathbf{w_{3}})^{2},\] \[\mathcal{L}_{0}^{(3)}(\mathbf{w_{3}})= \frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\bigg{(}\tilde{u}^{(3)}(\mathbf{x}_{0 }^{i},0^{i},\mathbf{w_{3}})-u_{0}(\mathbf{x}_{0}^{i})\bigg{)}^{2},\] \[\mathcal{L}_{b}^{(3)}(\mathbf{w_{3}})= \frac{1}{N_{b}}\sum_{i=1}^{N_{b}}\bigg{(}\tilde{u}^{(3)}(\mathbf{x}_{ b}^{i},t_{b}^{i},\mathbf{w_{3}})-g(\mathbf{x}_{b}^{i},t_{b}^{i})\bigg{)}^{2},\] \[\mathcal{L}_{pc,32}^{(1)}(\mathbf{w_{3}},\mathbf{w_{2}})= \frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\bigg{(}\tilde{u}^{(3)}(\mathbf{x}_ {c,3}^{i},t_{c,3}^{i},\mathbf{w_{3}})-\tilde{u}^{(2)}\big{(}\mathbf{x}_{c,2}^{i},t_{c,2}^{i},\mathbf{w_{2}}\big{)}\bigg{)}^{2},\] \[\mathcal{L}_{fc,32}^{(1)}(\mathbf{w_{3}},\mathbf{w_{2}})= \frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\bigg{(}\tilde{u}^{(3)}_{x}(\mathbf{ x}_{c,3}^{i},t_{c,3}^{i},\mathbf{w_{3}})-\tilde{u}^{(2)}_{x}\big{(}\mathbf{x}_{c,2}^{i},t _{c,2}^{i},\mathbf{w_{2}}\big{)}\bigg{)}^{2},\]

\[\mathcal{L}^{(3)}(\mathbf{w_{3}},\mathbf{w_{2}})= \mathcal{L}_{r}^{(3)}+\mathcal{L}_{0}^{(3)}+\mathcal{L}_{b}^{(3) }+\mathcal{L}_{pc,32}^{(3)}+\mathcal{L}_{fc,32}^{(3)}\]

For case 2.1, the well is controlled by bottom-hole pressure the results are displayed in Fig. 14, and for case 2.2, with well being controlled by constant production rate the results are displayed in Fig. 16. Fig. 15 shows the flow rates predictions for case 2.1 calculated using Darcy's equation 3.1 and the back-propagation to calculate \(\frac{d\tilde{p}}{dx}\).

Figure 14: Pressure predictions using the decomposition method for a circular reservoir with a well in the center controlled by constant BHP and a no-flow outer boundary versus high-fidelity (HF) numerical solution. The left figure is linear-linear and the right figure is semi-log. The relative L2 error is less than 0.2%

Figure 13: Decomposition of circular reservoir into 3 subdomains. The reservoir has a a well (black) in the center and no-flow outer boundary.
## 6 Conclusions and Future Work

In this work we showed that physics-informed neural network (PINN) can approximate partial differential equations to very good accuracy both for pressure predictions and flow rates predictions. In relation to reservoir engineering, a PINN model can be regarded as a good replacement for the analytical solution and therefore enables us to perform instantaneous predictions. Utilizing the inverse features of PINNs and a production flow rates history data we managed to estimate reservoir properties such as permeability and boundary distance without having to shut in the well. Alternatively, such quantities are estimated using pressure transient or rate transient tests which both require shutting-in the well for number of days.

For the 2D case, we proposed a decomposition method which is essential to overcome the stiffness of the PDE. With that we were able to train a model that can accurately predict the solution of the diffusivity equation. Finally, it is imperative to mention that the inverse applications that were applied to the 1D case can be applied to the 2D case.

## 7 Acknowledgements

Part of this work was done during my internship at OXY in summer 2022. I would like to thank OXY for providing all the reservoir simulations support needed to achieve these results.
