Predicting flow in porous media: a comparison of physics-driven neural network approaches

1. Introduction

Over the years, significant progress has been made in understanding multiscale physics through the numerical solution of partial differential equations. However, traditional numerical methods face challenges when dealing with real-world physics problems that have missing or noisy boundary conditions. In contrast, machine learning can explore large design spaces, identify correlations in multi-dimensional data, and manage ill-posed problems.

Deep learning tools [5] are particularly effective at extracting useful information from enormous amounts of data and linking important features with approximate models. They are especially useful for multi-dimensional subsurface flow problems because of their inverse nature. Nevertheless, current machine learning approaches often lack the ability to extract interpretable information and knowledge from the data, leading to poor generalization and physically inconsistent predictions.

To overcome these issues, physics-informed machine learning techniques [1-3] can incorporate physical constraints and domain knowledge into machine learning models. By “teaching” ML models about the governing physical rules, this approach can provide informative priors, such as strong theoretical constraints and inductive biases, that improve the performance of learning algorithms [4]. This process requires physics-based learning, which refers to the use of prior knowledge derived from observational, empirical, physical, or mathematical understanding of the world to improve machine learning algorithms' performance.

Machine learning methods can be categorized into supervised learning, unsupervised learning, and reinforcement learning. Supervised learning involves training a model using labeled input and output data and can be further divided into regression and classification. Examples of regression models include linear regression [6], Lasso [7], support vector regression [8], and random forest regression [9], while examples of classification models include logistic regression [10], naive Bayes classifier, and KNN classifier [11]. Unsupervised learning is used for unlabeled output data, with examples including principal component analysis [12] and singular value decomposition [13]. Reinforcement learning, on the other hand, learns from the environment based on reward and penalty.

In this paper we aim to apply Machine learning models to solve Single-phase Darcy’s flow equation. The Darcy flow equation is given as:

where $K$ is permeability or conductivity of the porous medium, $\upsilon$ – Darcy’s flow velocity, $p$ – pressure, $\rho$ – fluid density, $\mu$ – dynamic viscosity of the fluid, $g$ – gravitational constant, $z$ – spatial direction. The given equation can be solved by assuming constant porosity of the porous media domain and incompressibility, which reduces the equation into an elliptic pressure equation. The elliptic pressure equation’s solution can be approximated via machine learning models. The equation’s solution is relevant for understanding flow in subsurface hydrogeology, hydrocarbon reservoirs, and geothermal systems. In recent times researchers have tried to use machine learning models for modelling porous media flow, a comparative analysis of recent research work is presented in Table 1.

This paper addresses the challenges of numerical simulation for subsurface flow and transport problems in porous media, emphasizing the complexity of domain modeling required for accurate representation of real-world entities and relationships. Traditional numerical simulations involve solving PDEs on a network of networks, which increases computational complexity and necessitates the creation of large matrices. To address this, the paper proposes the use of machine learning as a cost-effective and time-saving alternative for solving these problems. The finite-volume approach generates data for the study, and Fig. 1 shows the permeability distribution and pressure in a 2×8 grid.

The strength of this approach lies in the comparative analysis of several machine learning models for solving single-phase Darcy’s flow equation. The objective is to develop a machine learning model that replicates numerical reservoir simulator outputs for single-phase flow problems, utilizing deep learning tools such as the Perceptron [14], Convolution Neural Network [15-22], and Artificial Neural Network models. The presented models are robust enough to predict pressure variation in the reservoir based on given permeability variations. The proposed approach removes challenges associated with numerical methods such as numerical diffusion, dispersion, and computational cost.

The paper is organized as follows: subsurface flow problem is described in Section 2. Flow equations are described in Section 3. Machine learning methodology used in this paper is described in Section 4. Results are presented in Section 5 and conclusions follow in Section 6.

Fig. 1a) Shows 1-D linear permeability distribution, b) pressure solution on a 2×8 grid

a)

b)

2. Single-phase Darcy flow equation

To derive an equation for flow in 1D, we assume that the cross-section area $A$ for flow as well as the depth $D$, are functions of the variable $x$ in our 1D space. Additionally, we introduce a term for the injection $q$ of fluid, which is equal to the mass rate of injection per unit volume of reservoir. Consider a mass balance in a small box shown in Fig. 2. The length of the box is $∆x$, the left side has area $A\left(x\right)$, the right side has area $A\left(x+∆x\right)$.

Fig. 2Differential elements of volume for a one-dimensional flow

Rate at which fluid mass enters the box at the left face is given by:

Rate at which fluid mass leaves at the right face:

If we define the average value of $A$ and $q$ between $x$ and $x+∆x$ as $A^{*}$and $q^{*}$, we get, that the volume of the box is $A^{*}∆x$.

Rate at which fluid mass is injected into the box is:

The mass contained in the box is ${\varphi }^{*}{\rho }^{*}{A}^{*}∆x$. So, we get the rate of accumulation of mass in the box:

From conservation of mass:

Putting values from Eq. (2-5) in Eq. (6), we yield:

Divided by $\mathrm{\Delta }x$, assuming area is constant, and taking limit $x\to 0$, we will get:

or

Here the source term $q$ models sources and sinks, that is, outflow and inflow per volume at designated well locations. For the incompressible single-phase flow assume that the porosity $\varphi$ of the rock is constant in time and that the fluid is incompressible, which means density is constant. Then the time dependent derivative vanishes, and we obtain the elliptic equation for the water pressure. Putting value of $v$ from Eq. (1) into Eq. (9), we will get:

When solving fluid flow problems, it is common practice to specify boundary conditions to ensure that the system is well-defined, and that the solution is physically meaningful. In this case, a no-flow boundary condition is imposed on the reservoir boundary, which means that the fluid velocity vector $v_w$ must be perpendicular to the boundary, as indicated by the dot product $v_w$, $n=$ 0, where n is the normal vector pointing out of the boundary. This condition ensures that no fluid can enter or exit the reservoir through the boundary, and the system is effectively isolated. Other types of boundary conditions, such as constant pressure or velocity, could also be used depending on the specific problem being solved.

2.1. Numerical solution: finite volume approach

In finite-difference methods, the domain is discretized into a grid of points and the partial derivatives are approximated using difference formulas at each grid point. The solution is then obtained by solving a system of algebraic equations obtained from discretizing the PDE at each grid point. In contrast, the finite-volume method (FVM) divides the computational domain into a finite number of small control volumes, or cells. The fluxes of fluid through the boundaries of each cell are then calculated, and the equation is discretized over the entire domain, see Fig. 3. The resulting set of algebraic equations can be solved using iterative methods to obtain a numerical solution. There are several variations of the FVM, such as the two-point flux method (TPFM) [29] and the multi-point flux method (MPFM) [29], which differ in the way they approximate the fluxes. In the TPFM, the fluxes across each face of a control volume are approximated using a two-point approximation, which assumes that the flux at a face is proportional to the pressure difference between the neighboring control volumes. The TPFM is based on the principle of conservation of flux, which states that the total flux into a control volume must be equal to the total flux out of the control volume.

We derive an equation that models flow of a fluid, say, water ($w$), through a porous medium characterized by a permeability field $K$ and $a$ corresponding porosity distribution. So, Eq. (10) can be written as:

11

$\nabla \left(v_w\right)=\frac{q_w}{{\rho }_w}.$

To derive a set of FV (finite volume) mass-balance equations for Eq. (11), denote a grid cell by ${\mathrm{\Omega }}_x$ and consider the following integral over ${\mathrm{\Omega }}_x$:

12

${\int }_{{\mathrm{\Omega }}_x}^{}\left(\frac{q_w}{{\rho }_w}\mathrm{}-\mathrm{}\nabla \left(v_w\right)\right)dx=\mathrm{}0.$

Fig. 3Notation of finite volume method

Applying divergence theorem:

13

${\int }_{{\mathrm{\Omega }}_x}^{}{v}_{w}ndv=\mathrm{}{\int }_{{\mathrm{\Omega }}_x}^{}\frac{q_w}{{\rho }_w}dx,$

or

14

$-{\int }_{\Omega }\nabla K\left(x,y\right)\nabla \left(p\right)ds={\int }_{\Omega }qds=M,\mathrm{\Omega }=\left[0\le x\le 1,0\le y\le 1\right],$

where ${\mathrm{\Omega }}_x$ is the interior, $\partial {\mathrm{\Omega }}_x$ is the boundary of control volume $x$ and $n$ denotes the outward-pointing unit normal on $\partial {\mathrm{\Omega }}_x$. On writing Eq. (12) in more compact form, we get:

15

$-\nabla \lambda \nabla u_w=\mathrm{}\frac{q_w}{{\rho }_w\mathrm{}},$

where $\lambda$ is mobility of water:

16

$\lambda =\mathrm{}\frac{K}{{\mu }_w}.$

Flow potential:

17

$u_w=\nabla p+{\rho }_{w}gh.$

Approximation of flux between $x$ and $y$:

18

$v_{xy}=\mathrm{}-{\int }_{{\gamma }_{xy}}^{}\left(\lambda \nabla u\right)ndv.$

The gradient $\nabla u$ on ${\gamma }_{xy}$ in the TPFA method have been replaced with:

19

$\delta u_{xy}=\mathrm{}\frac{2(u_y-\mathrm{}{u}_x)}{\mathrm{\Delta }{x}_x+\mathrm{}\mathrm{\Delta }{x}_y},$

where $\mathrm{\Delta }{x}_x$ and $\mathrm{\Delta }{x}_y$ denote the respective cell dimensions in the $x$-coordinate direction.

Permeability $K$ can be calculated using distance weighted harmonic average of the respective directional cell permeability, ${\lambda }_{x,\mathrm{}xy}=n_{xy}{\lambda }_x{n}_{xy}$ and${\mathrm{}\lambda }_{y,\mathrm{}xy}=n_{xy}{\lambda }_y{n}_{xy}$.

$n_{xy}$ directional permeability ${\lambda }_{xy}$ on ${\gamma }_{xy}$:

20

${\lambda }_{xy}=\left(\mathrm{\Delta }{x}_x+\mathrm{}\mathrm{\Delta }{x}_y\right){\left(\frac{\mathrm{\Delta }{x}_x}{{\lambda }_{x,\mathrm{}xy}}+\mathrm{}\frac{\mathrm{\Delta }{x}_y}{{\lambda }_{y,xy}}\right)}^{-1}.$

So:

21

$v_{xy}=\mathrm{}-{\int }_{{\gamma }_{xy}}^{}\left(\lambda \nabla u\right)n\mathrm{}dv={\gamma }_{xy}|\mathrm{}{\lambda }_{xy}\mathrm{}\delta u_{xy}=2|\mathrm{}{\gamma }_{xy}|\left(u_x-\mathrm{}{u}_y\right){\left(\frac{\mathrm{\Delta }{x}_x}{{\lambda }_{x,\mathrm{}xy}}+\mathrm{}\frac{\mathrm{\Delta }{x}_y}{{\lambda }_{y,xy}}\right)}^{-1}.$

Interface Transmissibility:

22

$t_{xy}=2\left|\mathrm{}{\gamma }_{xy}\right|{\left(\frac{\mathrm{\Delta }{x}_x}{{\lambda }_{x,\mathrm{}xy}}+\mathrm{}\frac{\mathrm{\Delta }{x}_y}{{\lambda }_{y,xy}}\right)}^{-1}.$

We get that Eq. (14) can be written as:

23

$\sum _y{t}_{xy}\left(u_x-u_y\right)=q.$

We will get linear equation in form of:

24

$\left[A\right]\left[P\right]=\left[Q\right],$

where $A$ is coefficient matrix, $A=\left[a_{xk}\right]$ where:

25

$a_{xk}=\left\{\begin{array}{c}\sum _y{t}_{xy}\mathrm{},\mathrm{}k=x,\\ -t_{xk},k\ne x,\end{array}\right.$

where $\left[P\right]$ is the vector of pressure at each node, and [$Q$] is the vector of flow rate at each node.

3. Machine learning method: solution methodology

The objective of this paper is to predict the pressure variation in reservoir with permeability using machine learning. It can be done using various methods, but ML is one of the best methods to do this task as conventional methods are time consuming. Fig. 4 shows workflow of the approach that we used to get Machine Learning result.

3.1. Example problem investigated using machine learning

In this paper, the given problem can be formulated as an elliptic Darcy PDE bounded by the 2-dimensional domain $\Omega$. Eq. (14) can be written into a given form:

$-{\int }_{\mathrm{\Omega }}\nabla K\left(x,y\right)\nabla \left(p\right)ds={\int }_{\mathrm{\Omega }}qds=M,\mathrm{\Omega }=\left[0\le x\le 1,0\le y\le 1\right].$

Fig. 4Methodology workflow diagram

Table 2Parameter and values

Parameter	Values
Pressure	Bar, 1 bar = 14.7 psi
Inflow rate	1 m3/sec
Outflow rate	1m3/sec
Permeability	Millidarcy
Length	1 m
Height	1 m
Width	1 m
Viscosity of water	1 cp

Table 2 presents the parameter values used in the problems discussed in this paper, while Table 3 shows the corresponding boundary conditions. The specific inflow and outflow rates used were both 1 m³/sec. The focus of the study was on investigating changes in pressure, which was specified at the corners of a simple 2D flow model using injection and production wells. The main objective was to solve equation 2 for pressure p over an arbitrary domain $\mathrm{\Omega }$, subject to appropriate boundary conditions (either Neumann or Dirichlet) on boundary $\partial \mathrm{\Omega }$. The term $M$ represents a specified flow rate, while $\nabla \mathrm{}=\left(\partial x,\partial y\right)$. The matrix $K$ can be a diagonal or a full cartesian tensor with a general form that can be expressed as:

26

$K=\left(\begin{array}{cc}{K}_{11}& K_{12}\\ K_{12}& K_{22}\end{array}\right).$

The full tensor pressure equation is assumed to be elliptic such that $K_{12}^2\le K_{11}{K}_{22}$. The tensor can be discontinuous across internal boundaries of $\mathrm{\Omega }$ ([25]).

Table 3Boundary condition

$p\left(\mathrm{0,0}\right)=$ 1	$p\left(\mathrm{1,1}\right)=$ –1
$\frac{dp(0,y)}{dx}=0$	$\frac{dp(1,y)}{dx}=0$	0 $<y<$ 1
$\frac{dp(x,0)}{dy}=0$	$\frac{dp(x,1)}{dy}=0$	0 $<x<$ 1

3.5. Feature scaling

Data scientists often assume that their data follows a Gaussian (normal) distribution, which forms the basis of many machine learning models. Therefore, it is essential to convert the dataset into a normal distribution form. This is achieved by applying appropriate transformations to the dataset. To verify normality in the dataset, we can use a histogram or QQ plot. Feature scaling is a technique that standardizes the independent or dependent features of the data to a fixed range. If feature scaling is not applied, the machine learning algorithm may assign higher weights to greater values and lower weights to smaller values, regardless of their units. Feature scaling is performed during data preprocessing to handle data with varying units, magnitudes or values.

Fig. 5Checking for missing values in grid 8×8

Fig. 6Workflow for Handling missing values

4. Results

Next, we present the results of our machine learning approaches. Table 4 shows the architecture of CNN model for 2×8 grid, 8×8 grid, 32×32 grid and 64×64 grid dimensions. We used Relu. Activation function in hidden layers for all dimensions of grid and linear activation in output layers for all dimension of grid. For 2×8 grid, which means that grid has 2 layers in $y$ direction and 8 cells in $x$ direction, support vector regression model works best with a mean absolute error of 0.69. SVR works better than linear regression because it uses kernel function, so our model can work on non-linear datasets too. Kernel function which we used in the model is radial basis function (RBF). Without the hyper-tuning Decision tree regression model having an over-fitting problem on 2×8 Grid, we get very low errors in the training dataset of order 10-16 and we get extremely high errors on validation and testing datasets. Best Parameters which minimize over-fitting are max depth = 3, min samples leaf = 1, and min samples split = 2. We found thar our model for comparison in higher dimension grid gives larger errors.

Fig. 7a) 2D view of permeability distribution using FVM, b) CNN predicted pressure, c) actual pressure using FVM and d) Lasso predicted pressure for 8×8 grid

4.1. Mean absolute and mean squared errors

The results are compared for different models on grid dimensions of 2×8, 8×8, 32×32, and 64×64 using mean absolute errors (presented in Figs. 8-11) and mean squared errors (presented in Figs. 11-15). Results indicated that the models performed well, and there was not a significant difference in mean absolute error between validation and testing. The 2D results obtained from machine learning models are compared to actual pressure on 8×8, 32×32, and 64×64 grids (shown in Fig. 7 and Fig. 16-17). The (a) part shows input permeability, the (c) part shows the actual pressure output obtained from the Finite-Volume approach, and the (b) and (d) parts show the pressure predicted through machine learning model.

Table 4Architecture of CNN model

Grid	No. of hidden layers	Activation		Kernel size	kernel_regularizer	Learning rate
		Hidden layer	Output layers
2×8	5 (4 CNN layers (filters=512,256, 128,64) and 1 Dense layers (filters = 32))	Relu	Linear	3	L2 Regularizer with an alpha value of .01 at the Dense Layer	1e-06
8×8	6 (5 CNN layers (filters=1024,512,256 ,128,64) and 1 Dense layers (filters = 32))	Relu	Linear	3	L2 Regularizer with an alpha value of .01 at the last layer of CNN layers	1e-03
32×32	4 (3 CNN layers (filters = 256, 128, 64) and 1 Dense layers)	Relu	Linear	3	L2 Regularizer with an alpha value of .01 at the last layer of CNN layers	1e-06
64×64	3 (2 CNN layers (filters = 128, 64) and 1 Dense layers (filters = 32))	Relu	Linear	3	L2 Regularizer with an alpha value of .01 at the last layer of CNN layers	1e-03

Fig. 8Mean absolute error on 2×8 grid

Fig. 9Mean absolute error on 8×8 grid

Fig. 10Mean absolute error on 32×32 grid

Fig. 11Mean absolute error on 64×64 grid

Fig. 12Mean squared error on 2×8 grid

Fig. 13Mean squared error on 8×8 grid

Fig. 14Mean squared error on 32×32 grid

Fig. 15Mean squared error on 64×64 grid

Fig. 16a) 2D view of permeability distribution, b) CNN predicted pressure, c) actual pressure, d) lasso predicted pressure for 32×32 grid

a)

b)

c)

d)

Fig. 17a) 2D view of permeability distribution, b) ANN predicted pressure, c) actual pressure, d) CNN predicted pressure for 64×64 grid

a)

b)

c)

d)

4.2. Residual errors

To verify our result, we used residual plot. The residual plot shows the difference between the observed response and the fitted response values. Simply if our model can predict actual pressure then it shows zero difference, in our plot this can be shown by white color. The residual plots are shown for 8×8 grid, 32×32 and 64×64 grids for different methods that are presented in Figs. 18-21.

Fig. 18Residual plot for CNN, ANN, random forest and decision tree regression for 8×8 on a grid

Fig. 19Residual plot for perceptron, support vector regression, linear regression, Lasso regression, light weight gradient, eXtreme Gradient Boosting on an 8×8 grid

Fig. 20Residual plot for a) random forest, b) eXtreme gradient boosting, c) Lasso regression, d) CNN model, e) ANN Model, f) support vector regression model for 32×32 grid

Fig. 21Residual plot for CNN, ANN, perceptron, linear regression, random forest regression, eXtreme gradient boosting, Lasso regression, light weight gradient boosting for 64×64 grid

5. Conclusions

The paper begins by generating a material balance equation and defining a finite volume approach mathematically for Darcy's equation. To predict pressure variation in a reservoir based on permeability variation in different grid sizes, several machine learning approach have been tested. Ten different machine learning models, including linear regression, support vector regression, random forest regression, lasso regression, light weight gradient boosting, extreme gradient boosting, decision tree regression, perceptron, artificial neural network (ANN), and convolutional neural network (CNN) were used to predict pressure. The analysis shows that CNN, linear regression, and lasso regression had the lowest mean absolute error (MAE) and mean squared error (MSE) values, while the decision tree regression model had the highest error. The CNN model demonstrated the best performance on the validation dataset for most grid dimensions, and both the CNN and ANN models performed well on a 64×64 grid. In conclusion, the lower the error values, the better the model performance, and the CNN and ANN models demonstrated good performance for this particular prediction task.