Learning macroscopic parameters in nonlinear multiscale simulations using nonlocal multicontinua upscaling techniques

TL;DR

This paper introduces a machine learning-based nonlocal upscaling method for nonlinear multiscale simulations of flow in heterogeneous and fractured porous media, improving accuracy and computational efficiency.

Contribution

It develops a novel deep learning approach for computing macroscopic parameters in nonlocal multicontinua models, enhancing previous global machine learning methods.

Findings

High accuracy in coarse grid approximations

Fast computation of effective properties

Effective handling of complex nonlinear systems

Abstract

In this work, we present a novel nonlocal nonlinear coarse grid approximation using a machine learning algorithm. We consider unsaturated and two-phase flow problems in heterogeneous and fractured porous media, where mathematical models are formulated as general multicontinuum models. We construct a fine grid approximation using the finite volume method and embedded discrete fracture model. Macroscopic models for these complex nonlinear systems require nonlocal multicontinua approaches, which are developed in earlier works. These rigorous techniques require complex local computations, which involve solving local problems in oversampled regions subject to constraints. The solutions of these local problems can be replaced by solving original problem on a coarse (oversampled) region for many input parameters (boundary and source terms) and computing effective properties derived by nonlinear nonlocal multicontinua approaches. The effective properties depend on many variables (oversampled region and the number of continua), thus their calculations require some type of machine learning techniques. In this paper, our contribution is two fold. First, we present macroscopic models and discuss how to effectively compute macroscopic parameters using deep learning algorithms. The proposed method can be regarded as local machine learning and complements our earlier approaches on global machine learning. We consider a coarse grid approximation using two upscaling techniques with single phase upscaled transmissibilities and nonlocal nonlinear upscaled transmissibilities using a machine learning algorithm. We present results for two model problems in heterogeneous and fractured porous media and show that the presented method is highly accurate and provides fast coarse grid calculations.

Tables

Table 1. Table 1: Learning performance of machine learning algorithm for Test 1 (nonlinear flow). Errors for train and test sets

	MSE	RMSE (%)	MAE (%)
Train set (global data)
$N{N}_1$	0.012	1.101	1.072
$N{N}_2$	0.016	1.273	1.272
$N{N}_3$	0.013	1.156	0.838
Test set (global data)
$N{N}_1$	0.012	1.104	1.085
$N{N}_2$	0.016	1.286	1.280
$N{N}_3$	0.011	1.050	0.774
Train set (local data)
$N{N}_1^l$	0.081	2.861	2.060
$N{N}_2^l$	0.014	1.223	1.229
$N{N}_3^l$	0.042	2.058	1.791

Table 2. Table 2: Learning performance of machine learning algorithm for Test 2 (nonlinear flow and transport). Errors for train and test sets

	MSE	RMSE (%)	MAE (%)
$N{N}_1$	0.017	1.316	0.959
$N{N}_2$	0.043	2.092	1.507
$N{N}_3$	0.014	1.218	0.778
$N{N}_4$	0.052	2.301	1.328

Equations

$$\frac{\partial \Theta }{\partial t}-\nabla \cdot (k(x,p)\nabla (p+z))=f,x\in \Omega ,$$

$$\begin{array}{rl}& \frac{\partial {\Theta }^m}{\partial t}-\nabla \cdot (k^m(x,p^m)\nabla (p^m+z))+{\sigma }^{mf}(x,p^{mf})(p^m-p^f)=f^m,x\in \Omega ,\\ & \frac{\partial {\Theta }^f}{\partial t}-{\nabla }_{\gamma }\cdot (k^f(x,p^f){\nabla }_{\gamma }(p^f+z))-{\sigma }^{fm}(x,p^{mf})(p^m-p^f)=f^f,x\in \gamma ,\end{array}$$

$$\begin{array}{rl}& \frac{\partial {\Theta }^{\alpha }}{\partial t}-\nabla \cdot (k^{\alpha }(x,p^{\alpha })\nabla (p^{\alpha }+z))+\sum _{\beta }{\sigma }^{\alpha \beta }(x,p^{\alpha \beta })(p^{\alpha }-p^{\beta })=f^{\alpha },\end{array}$$

$$c^{\alpha }(x,p^{\alpha })\frac{\partial p^{\alpha }}{\partial t}-\nabla \cdot (k^{\alpha }(x,p^{\alpha })\nabla p^{\alpha })+\sum _{\beta }{\sigma }^{\alpha \beta }(x,p^{\alpha \beta })(p^{\alpha }-p^{\beta })=q^{\alpha },$$

$${\mathcal{T}}^h={\cup }_{i=1}^{N^{m,h}}{\varsigma }_i,{\mathcal{G}}^h={\cup }_{l=1}^{N^{f,h}}{\iota }_l,$$

$$\begin{array}{rl}& c_i^m\frac{p_i^m-{\check{p}}_i^m}{\tau }|{\varsigma }_i|+\sum _j{u}_{ij}^{mm}+\sum _l{u}_{il}^{mf}=q_i^m|{\varsigma }_i|,\forall i=1,...,N^{m,h}\\ & c_l^f\frac{p_l^f-{\check{p}}_l^f}{\tau }|{\iota }_l|+\sum _n{u}_{ln}^{ff}+\sum _i{u}_{il}^{fm}=q_l^f|{\iota }_l|,\forall l=1,..,N^{f,h},\end{array}$$

$$u^{mm}=-k^m(x,p^m)\nabla p^m,u^{ff}=-k^f(x,p^f){\nabla }_{\gamma }{p}^f,$$

$$u_{ij}^{mm}=u^{mm}\cdot n{|}_{E_{ij}}=T_{ij}^{mm}(p_i^m,p_j^m)(p_i^m-p_j^m),u_{ln}^{ff}=u^{ff}\cdot n{|}_{e_{ln}}=T_{ij}^{ff}(p_i^f,p_l^f)(p_i^f-p_j^f),$$

$$u^{mf}={\sigma }^{mf}(x,p^{mf})(p^m-p^f),u^{fm}={\sigma }^{fm}(x,p^{mf})(p^f-p^m),$$

$$u_{il}^{mf}=-u_{il}^{fm}=T_{il}^{mf}(p_i^m-p_l^f),$$

$$\begin{array}{rl}& c_i^m\frac{p_i^m-{\check{p}}_i^m}{\tau }|{\varsigma }_i|+\sum _j{T}_{ij}^{mm}(p_i^m,p_j^m)(p_i^m-p_j^m)+\sum _l{T}_{il}^{mf}(p_i^m,p_l^f)(p_i^m-p_l^f)=q_i^m|{\varsigma }_i|,\forall i=1,...,N^{m,h}\\ & c_l^f\frac{p_l^f-{\check{p}}_l^f}{\tau }|{\iota }_l|+\sum _n{T}_{ln}^{ff}(p_l^f,p_n^f)(p_l^f-p_n^f)-\sum _i{T}_{il}^{mf}(p_i^m,p_l^f)(p_i^m-p_l^f)=q_l^f|{\iota }_l|,\forall l=1,..,N^{f,h},\end{array}$$

$$c_i^{\alpha }\frac{p_i^{\alpha }-{\check{p}}_i^{\alpha }}{\tau }|{\varsigma }_i|+\sum _j{T}_{ij}^{\alpha \alpha }(p_i^{\alpha },p_j^{\alpha })(p_i^{\alpha }-p_j^{\alpha })+\sum _{\beta }\sum _l{T}_{il}^{\alpha \beta }(p_i^{\alpha },p_l^{\beta })(p_i^{\alpha }-p_l^{\beta })=q_i^{\alpha }|{\varsigma }_i|,\forall i=1,...,N^{\alpha ,h},$$

$$\begin{array}{r}\varphi \frac{\partial s}{\partial t}-\nabla \cdot ({\lambda }^w(s)k(x)\nabla p)=q^w,x\in \Omega ,\\ -\nabla \cdot (\lambda (s)k(x)\nabla p)=q,x\in \Omega ,\end{array}$$

$$\begin{array}{r}{\varphi }^m\frac{\partial s^m}{\partial t}-\nabla \cdot ({\lambda }^w(s^m)k^m(x)\nabla p^m)+{\lambda }^w(s^{mf}){\sigma }^{mf}(x)(p^m-p^f)=q^{w,m},x\in \Omega ,\\ -\nabla \cdot (\lambda (s^m)k^m(x)\nabla p^m)+\lambda (s^{mf}){\sigma }^{mf}(x)(p^m-p^f)=q^m,x\in \Omega ,\\ {\varphi }^f\frac{\partial s^f}{\partial t}-{\nabla }_{\gamma }\cdot ({\lambda }^w(s^f)k^f(x){\nabla }_{\gamma }{p}^f)-{\lambda }^w(s^{mf}){\sigma }^{fm}(x)(p^m-p^f)=q^{w,f},x\in \gamma ,\\ -{\nabla }_{\gamma }\cdot (\lambda (s^f)k^f(x){\nabla }_{\gamma }{p}^f)+\lambda (s^{mf}){\sigma }^{fm}(x)(p^m-p^f)=q^f,x\in \gamma ,\end{array}$$

$$\begin{array}{r}{\varphi }^{\alpha }\frac{\partial s^{\alpha }}{\partial t}-\nabla \cdot ({\lambda }^w(s^{\alpha })k^{\alpha }(x)\nabla p^{\alpha })+\sum _{\beta }{\lambda }^w(s^{\alpha \beta }){\sigma }^{\alpha \beta }(x)(p^{\alpha }-p^{\beta })=q^{w,\alpha },\\ -\nabla \cdot (\lambda (s^{\alpha })k^{\alpha }(x)\nabla p^{\alpha })+\sum _{\beta }\lambda (s^{\alpha \beta }){\sigma }^{\alpha \beta }(x)(p^{\alpha }-p^{\beta })=q^{\alpha },\end{array}$$

$$\begin{array}{rl}{\varphi }_i^m\frac{s_i^m-{\check{s}}_i^m}{\tau }|{\varsigma }_i|+& \sum _j{u}_{ij}^{w,mm}+\sum _l{u}_{il}^{w,mf}=(1-f_i^{w,m})q_i^{w,m}|{\varsigma }_i|,\\ & \sum _j{u}_{ij}^{mm}+\sum _l{u}_{il}^{mf}=q_i^m|{\varsigma }_i|,\\ {\varphi }_l^f\frac{s_l^f-{\check{s}}_l^f}{\tau }|{\iota }_l|+& \sum _n{u}_{ln}^{w,ff}+\sum _i{u}_{il}^{w,fm}=(1-f_l^{w,f})q_l^{w,f}|{\iota }_l|,\\ & \sum _n{u}_{ln}^{ff}+\sum _i{u}_{il}^{fm}=q_l^f|{\iota }_l|,\end{array}$$

$$u^{mm}=-\lambda (s^m)k^m(x)\nabla p^m,u^{w,mm}=-{\lambda }^w(s^m)k^m(x)\nabla p^m,$$

$$u^{ff}=-\lambda (s^f)k^f(x){\nabla }_{\gamma }{p}^f,u^{w,ff}=-{\lambda }^w(s^f)k^f(x){\nabla }_{\gamma }{p}^f,$$

$$u^{mf}=\lambda (s^{mf}){\sigma }^{mf}(x)(p^m-p^f),u^{w,mf}={\lambda }^w(s^{mf}){\sigma }^{mf}(x)(p^m-p^f),$$

$$u^{fm}=\lambda (s^{mf}){\sigma }^{fm}(x)(p^f-p^m),u^{w,fm}={\lambda }^w(s^{mf}){\sigma }^{fm}(x)(p^f-p^m),$$

$$u_{ij}^{mm}=u^{mm}\cdot n{|}_{E_{ij}}=T_{ij}^{mm}(p_i^m-p_j^m),u_{ij}^{w,mm}=u^{w,mm}\cdot n{|}_{E_{ij}}=T_{ij}^{w,mm}(p_i^m-p_j^m),$$

$$u_{ln}^{ff}=u^{ff}\cdot n{|}_{e_{ln}}=T_{ij}^{ff}(p_i^f-p_j^f),u_{ln}^{w,ff}=u^{w,ff}\cdot n{|}_{e_{ln}}=T_{ij}^{w,ff}(p_i^f-p_j^f),$$

$$u_{il}^{mf}=-u_{il}^{fm}=T_{il}^{mf}(p_i^m-p_l^f),u_{il}^{w,mf}=-u_{il}^{w,fm}=T_{il}^{w,mf}(p_i^m-p_l^f),$$

$$T_{ij}^{\alpha \beta }=T_{ij}^{\alpha \beta }(s_i^{\alpha },s_j^{\beta })=\lambda (s_{ij}^{\alpha \beta })W_{ij}^{\alpha \beta },T_{ij}^{w,\alpha \beta }=T_{ij}^{w,\alpha \beta }(s_i^{\alpha },s_j^{\beta })={\lambda }^w(s_{ij}^{\alpha \beta })W_{ij}^{\alpha \beta },\alpha ,\beta =m,f$$

$$\begin{array}{rl}{\varphi }_i^m\frac{s_i^m-{\check{s}}_i^m}{\tau }|{\varsigma }_i|+& \sum _j{T}_{ij}^{w,mm}({\check{s}}_i^m,{\check{s}}_j^m)(p_i^m-p_j^m)+\sum _l{T}_{il}^{w,mf}({\check{s}}_i^m,{\check{s}}_l^f)(p_i^m-p_l^f)=(1-f_i^{w,m})q_i^{w,m}|{\varsigma }_i|,\\ & \sum _j{T}_{ij}^{mm}({\check{s}}_i^m,{\check{s}}_j^m)(p_i^m-p_j^m)+\sum _l{T}_{il}^{mf}({\check{s}}_i^m,{\check{s}}_l^f)(p_i^m-p_l^f)=q_i^m|{\varsigma }_i|,\\ {\varphi }_l^f\frac{s_l^f-{\check{s}}_l^f}{\tau }|{\iota }_l|+& \sum _n{T}_{ln}^{w,ff}({\check{s}}_l^f,{\check{s}}_n^f)(p_l^f-p_n^f)-\sum _i{T}_{il}^{w,mf}({\check{s}}_i^m,{\check{s}}_l^f)(p_i^m-p_l^f)=(1-f_l^{w,f})q_l^{w,f}|{\iota }_l|,\\ & \sum _n{T}_{ln}^{ff}({\check{s}}_l^f,{\check{s}}_n^f)(p_l^f-p_n^f)-\sum _i{T}_{il}^{mf}({\check{s}}_i^m,{\check{s}}_l^f)(p_i^m-p_l^f)=q_l^f|{\iota }_l|,\end{array}$$

$${\lambda }^w(s_{ij}^{\alpha \beta })=\{\begin{array}{ll}{\lambda }^w(s_i^{\alpha }),& \text{mbox}if{T}_{ij}^{\alpha \beta }({\check{s}}_i^{\alpha },{\check{s}}_j^{\beta })(p_i^{\alpha }-p_j^{\beta })>0\\ {\lambda }^w(s_j^{\beta }),& else,\end{array}.$$

$$\begin{array}{rl}{\varphi }_i^{\alpha }\frac{s_i^{\alpha }-{\check{s}}_i^{\alpha }}{\tau }|{\varsigma }_i|+& \sum _j{T}_{ij}^{w,\alpha \alpha }({\check{s}}_i^{\alpha },{\check{s}}_j^{\alpha })(p_i^{\alpha }-p_j^{\alpha })+\sum _{\beta }\sum _l{T}_{il}^{w,\alpha \beta }({\check{s}}_i^{\alpha },{\check{s}}_l^{\beta })(p_i^{\alpha }-p_l^{\beta })=(1-f_i^{w,\alpha })q_i^{w,\alpha }|{\varsigma }_i|,\\ & \sum _j{T}_{ij}^{\alpha \alpha }({\check{s}}_i^{\alpha },{\check{s}}_j^{\alpha })(p_i^{\alpha }-p_j^{\alpha })+\sum _{\beta }\sum _l{T}_{il}^{\alpha \beta }({\check{s}}_i^{\alpha },{\check{s}}_l^{\beta })(p_i^{\alpha }-p_l^{\beta })=q_i^{\alpha }|{\varsigma }_i|,\end{array}$$

$${\mathcal{T}}^H={\cup }_{i=1}^{N^H}{K}_i,$$

$${\omega }_{ij}=K_i\cup K_j,K_i,K_j\in {\mathcal{T}}^H,$$

$$k^{\alpha }(x,p^{\alpha })=k_r(p^{\alpha })k_s^{\alpha }(x),{\sigma }^{\alpha \beta }(x,p^{\alpha \beta })={\sigma }_r(p^{\alpha \beta }){\sigma }_s^{\alpha \beta }(x),\alpha =m,f,$$

Full text

1 Introduction

Mathematical models of the flow and transport problems in heterogeneous and fractured porous media are required to solve large and complex nonlinear systems. Processes in fractured porous media are described by the mixed dimensional coupled system of equations [27, 13, 17, 24, 30]. Such models can be generalized as a general multicontinuum models similar to the dual porosity/dual permeability approaches [7, 32].

Solving problems in heterogeneous and fractured media requires constructing grids that resolve all small scale heterogeneity. Numerous model order reduction techniques have been developed to construct coarse grid approximations and reduce the computational time of the numerical simulations. Global model order reduction approaches rely on projection on the important modes space, where the full-order model is carefully study generating POD based basis functions to perform the fast online calculations stage [16, 40]. Local model reduction techniques are based on constructing local multiscale basis functions to represent the influence of small scale heterogeneity [30, 4, 34, 35]. One of the widely used ways is based on the numerical homogenization technique, where effective parameters are calculated in order to construct coarse grid approximations [2, 29, 37]. The coarse grid parameters are constructed by solving local problems with appropriate boundary conditions. For example, linear boundary conditions or periodicity can be used. The choice of boundary conditions have a strong impact on the accuracy of results. In [5], an interpolated global coarse grid solution is used for performing accurate construction of the upscaled transmissibilities, which involve iterations between global coarse grid model and local fine grid calculations. In standard upscaling methods, upscaled parameters are obtained independent of any global problem. However, these approaches lack several features, which are important for rigorous and accurate upscaling. These include the use of multiple continua and oversampled computations.

In the local model order reduction methods, an oversampling technique and multicontinua concepts are needed to achieve an accuracy independent of physical parameters, such as scales and contrast [11, 41]. For example, the oversampled domain is used for constructing multiscale basis and provide more accurate results in the Multiscale Finite Element Method [20]. In the Generalized Multiscale Finite Element Method [15, 9], a larger domain is used to construct a space of snapshots and solution of the local spectral problem to determine a dominant modes. Note that, the oversampled domain is used for local problem solution and only the interior information is used to define the basis functions.

In recently developed Constrained Energy Minimization and Nonlocal Multicontinuum methods [10, 12, 33, 36], multiscale basis functions are defined in the oversampled domains and constructed via solving local constrained energy minimization problems, where constraints are related to each continuum. Continuum plays a role of macroscopic parameter. In [5, 14], oversampling techniques have been developed in the context of the upscaling procedure, where an interesting local-global upscaling technique is presented for constructing coarse scale approximation for highly heterogeneous porous media. In this method, the coarse grid simulations are iterated with local calculations of the upscaled parameters, where the coarse grid solutions are used to determine the boundary conditions for the local calculation. The local-global upscaling method requires more computation than existing classic upscaling procedures. Therefore, in the upscaling and multiscale methods oversampled domains are used in two contexts: (1) as extended local domains for more accurate calculations of the coarse grid parameters with fine-scale information about heterogeneous properties and (2) for global or quasi-global information of the solutions that are used, for example, as boundary conditions in local calculations. The second context of the oversampling technique, due to the incorporation of solution information into the local problems leads to the nonlinear equations even in the case of the linear problems.

In this work, we consider flow and transport processes in heterogeneous multicontinuum media and construct coarse and fine grid approximations using a finite volume method with two-point flux approximation. Due to nonlinear nature of these flows, upscaled parameters are nonlinear functions, which depend on multiple coarse-grid variables defined in oversampled regions. Macroscopic equations use nonlinear nonlocal multicontinuum concept [8]. This framework first identifies macroscopic variables in each coarse-grid block and then solves local constraint problems in oversampled regions to compute macroscopic fluxes. The local oversampled computations require solving nonlinear problems with constraints that include the values of macroscopic variables. For example, for two-phase flow simulations, this requires solving two-phase flow problems with known values of pressures and saturations in each coarse-grid block. Each coarse-grid block may contain several macroscale pressures and saturations identified in the first step. Solving local nonlinear constraint problems can be challenging due to large number of nonlinear simulations in oversampled regions. Moreover, computing macroscale fluxes as a function of many variables as a look-up table is nearly impossible. In this work, we propose an efficient algorithm for solving the local problems consisting of original problems with various boundary conditions and using deep learning to train macroscale fluxes. This is a first step in designing computationally efficient and rigorous upscaling methods for nonlinear flows in porous media.

Constructing accurate upscaled transmissibilities for the coarse grid approximation is based on the information about solution (nonlinear transmissibilities). The presented method is based on the machine learning procedure for fast prediction of the nonlinear transmissibilities, where we construct neural networks that learn dependencies between the coarse grid quantities on the oversampled local domains and upscaled transmissibilities. We use a convolutional neural network and GPU training process to construct a machine learning algorithm [23, 22]. For constructing the datasets, we perform local or global calculations of the coarse grid quantities [5]. In the local approach, the upscaled transmissibilities are calculated on the local domain corresponding to the target face, where the fine-scale solution information is used to set boundary conditions. The global approach uses a global fine-scale solution for determining coarse scale parameters. For training the neural networks, we use a family of problem solutions for different input conditions. Note that, we should have many snapshots to capture all input condition variations because accuracy of the machine learning method depends on space of snapshots that is used as a train dataset. As soon as neural networks trained on the dataset, the fast and accurate calculations can be performed. To illustrate method construction and applicability, we considered two model problems: unsaturated flow problem and two-phase filtration problem in heterogeneous and fractured porous media. The presented method combines accuracy of the fine grid models with fast coarse grid calculations by constructing machine learning techniques for predicting accurate nonlinear upscaled transmissibilities.

The paper is organized as follows. In Section 2, we present the mathematical model and fine grid approximation. In Section 3, we consider single-phase upscaling for problems in multicontiuum media with variable separation for nonlinear and space dependent variables. Next, we present a novel nonlinear coarse grid approximation using a machine learning algorithm in Section 4. In Section 5, we consider two model problems in two - dimensional formulation and present numerical results, where we consider training of machine learning algorithm and relative errors between the reference fine grid solution and presented method. Finally, we present conclusions.

2 Model problem with reference fine grid approximation

As a model problem, we consider two nonlinear problems for fractured and heterogeneous porous media:

Unsaturated flow problem (nonlinear flow problem) 2. 2.

Two-phase flow problem (nonlinear transport and flow problem)

We start with the formulation of the mathematical model, where we formulate models for fractured media and generalize it for multicontinuum media. Next, we present a construction of the fine grid approximation using finite volume approximation and embedded fracture model.

2.1 Unsaturated flow problem

Mathematical model of the unsaturated flow in porous media described by the Richards’ equations [28]

[TABLE]

where $p$ is the pressure head, $k$ is the unsaturated hydraulic conductivity tensors, $z$ represent the influence of the gravity to the flow processes, $\Theta$ is the water content and $f$ refer to source and sink terms.

For fractured porous media, we consider a mixed dimensional formulation of the flow problem [27, 13, 17]. Let $\Omega\in\mathcal{R}^{d}$ is the $d$ - dimensional domain of the porous matrix, where $d=2,3$. Fracture network is considered as a $(d-1)$ - dimensional (lower dimensional) domain $\gamma\in\mathcal{R}^{d-1}$ due to small thickness of the fractures compared to the domain sizes. Then, for unsaturated flow in fractured porous media, we have the following coupled system of equations for $p^{m}$ and $p^{f}$:

[TABLE]

where $p^{m}$ and $p^{f}$ are the pressure head in matrix and fractures; $k^{m}$ and $k^{f}$ are the unsaturated hydraulic conductivity tensors for matrix and fractures; $z$ represent the influence of the gravity to the flow processes; $\nabla_{\gamma}$ contains partial derivativies along fracture $\gamma$; $\Theta^{m}$ and $\Theta^{f}$ are the water content for matrix and fracture; and $f^{m}$ and $f^{f}$ refer to source and sink terms. For transfer term, we have $\int_{V}\sigma^{mf}(p^{m}-p^{f})dx=\int_{A}\sigma^{fm}(p^{m}-p^{f})ds$, $\sigma^{mf}=\sigma/V$, $\sigma^{fm}=\sigma/A$ ($\sigma=CI\,k^{mf}$ and $k^{mf}(x,p^{mf})$ is the harmonic average between $k^{m}(x,p^{m})$ and $k^{f}(x,p^{f})$) for matrix volume $V$ intersecting with fracture surface and $CI$ is the connectivity index [4, 30]. As an initial condition, we set $p^{\alpha}=p_{0}$ ($\alpha=m,f$) and zero flux boundary conditions on $\partial\Omega$ and $\partial\gamma$.

In general, we have following multicontinuum model:

[TABLE]

where $\alpha=1,...,M$ and $M$ is the number of continuum.

Let $c^{\alpha}(x,p^{\alpha})=\partial\Theta^{\alpha}/\partial p$ therefore we have following coupled system of nonlinear parabolic equations

[TABLE]

where $q^{\alpha}=f^{\alpha}+\nabla\cdot(k^{\alpha}(x,p^{\alpha})z)$, $c^{\alpha}(p^{\alpha})$ and $k^{\alpha}(x,p^{\alpha})$ are the nonlinear coefficient ($\alpha=1,...,M$).

For the approximation on the fine grid, we use structured grids with embedded discrete fracture model (EDFM) [24, 30]. Let $\mathcal{T}^{h}$ denote a structured fine grid of the porous matrix domain $\Omega$ and $\mathcal{G}^{h}$ denote a fine grid of the fracture domain $\gamma$

[TABLE]

where $\varsigma_{i}$ and $\iota_{l}$ are the cell of the matrix and fractures fine grids, $N^{m,h}$ is the number of cells in $\mathcal{T}^{h}$, $N^{f,h}$ is the number of cell related to fracture mesh $\mathcal{G}^{h}$. Therefore, for finite volume approximation we have

[TABLE]

where $p_{i}^{m}$ and $p_{l}^{f}$ are pressure of matrix and fracture continuum in cells $\varsigma_{i}$ and $\iota_{l}$, $|\varsigma_{i}|$ and $|\iota_{l}|$ are the volume of cells. Here, we use an implicit scheme for time discretization, where $\check{p}^{m}_{i}$ and $\check{p}^{f}_{l}$ are the solutions from previous time step and $\tau$ is the given time step [34, 35].

For approximation of the flux in matrix ($u^{mm}$) and fracture ($u^{ff}$) continuum

[TABLE]

we use a classic two point flux approximation (TPFA)

[TABLE]

where $T_{ij}^{mm}=k^{m}_{ij}|E_{ij}|/d_{ij}$, $T_{ln}^{ff}=k^{f}_{ln}/\Delta_{ln}$, $k^{m}_{ij}=(k_{i}^{m}(p_{i}^{m})+k_{j}^{m}(p_{j}^{m}))/2$, $k^{f}_{ln}=(k_{l}^{f}(p_{l}^{f})+k_{n}^{f}(p_{n}^{f}))/2$, $E_{ij}$ and $e_{ln}$ are the interface between two cells, $|E_{ij}|$ is the length of face between cells $\varsigma_{i}$ and $\varsigma_{j}$, $d_{ij}$ is the distance between midpoint of cells $\varsigma_{i}$ and $\varsigma_{j}$, $\Delta_{ln}$ is the distance between midpoint of cells $\iota_{l}$ and $\iota_{n}$.

For the flux between matrix and fracture continuum

[TABLE]

we follow EDFM and have the following approximation

[TABLE]

where $T^{mf}_{il}=\sigma_{il}$ with $\sigma_{il}=CI_{il}k^{mf}_{il}$ if $\iota_{l}\subset\varsigma_{i}$ and zero else ($CI_{il}$ is the connectivity index from [24, 30, 18] that proportional to the distance and area of the intersection between the fracture cell $\iota_{l}$ and porous matrix cell $\varsigma_{i}$).

Therefore, we have following fine grid approximation

[TABLE]

or for general multicontinuum case [7, 32], we have

[TABLE]

where $\alpha=1,...,M$ and $M$ is the number of continua.

2.2 Two-phase flow problem

Mathematical model of the two-phase flow problem in porous media contains a conservation law and Darcy’s law [19]. For the case with incompressible fluid and rock and without gravitational and capillary forces, we have

[TABLE]

where $s=s^{w}$ is the saturation of the wetting phase, $p$ is the pressure, $q=q^{w}+q^{n}$, $q^{w}$ and $q^{n}$ are the source/sink of wetting and nonwetting phases, $\lambda^{i}={k_{ri}(s)}/{\mu_{i}}$, $\lambda=\lambda^{n}+\lambda^{w}$, $\phi$, $k$ are the porosity and permeability, $\mu_{i}$ and $k_{ri}$ are the viscosity and relative permeability for $i$-phase ($i=n,w$).

For the fractured porous media, we consider the mixed dimensional mathematical model for two-phase flow problem

[TABLE]

where $s^{m}$ and $s^{f}$ are the saturation in porous matrix and in fractures; $p^{m}$ and $p^{f}$ are the pressure in matrix and in fractures; $\phi^{m}$ and $\phi^{f}$ are the porosity for matrix and fracture continuum; $k^{m}$ and $k^{f}$ are the absolute permeability of matrix and fractures; and $q^{\alpha}=q^{w,\alpha}+q^{n,\alpha}$, $q^{w,\alpha}$ and $q^{n,\alpha}$ are the source/sink of wetting and nonwetting phases for continua $\alpha=m,f$. On the fine grid, we suppose $\lambda^{i,m}=\lambda^{i,f}=\lambda^{i}$ for $i=n,w$ and $\lambda=\lambda^{n}+\lambda^{w}$, where $\lambda^{i}={k_{ri}(s^{\alpha})}/{\mu_{i}}$, $\mu_{i}$ is the viscosity and $k_{ri}$ is relative permeability for $i$-phase that depends on flux direction. In general, realtive permeability functions can be different for each continua and for flux between them. For transfer term, similarly to the previous model for unsaturated flow, we have $\sigma^{mf}=\sigma/V$, $\sigma^{fm}=\sigma/A$ ($\sigma=CI\,k^{mf}$ and $k^{mf}$ is the harmonic average between $k^{m}$ and $k^{f}$) for matrix volume $V$ intersecting with fracture surface and $CI$ is the connectivity index [4, 30]. As an initial condition, we set $s^{\alpha}=s_{0}(x)$ ($\alpha=m,f$). For boundary conditions, we set zero flux on $\partial\Omega$ and $\partial\gamma$.

For the general multicontinuum case, we can write

[TABLE]

where $\alpha=1,...,M$ and $M$ is the number of continuum.

Similarly to the previous model problem of unsaturated flow, we use structured grids and construct a finite volume approximation with embedded discrete fracture model (EDFM) for approximation of the two-phase flow problem. We use same fine grid for the porous matrix domain ($\mathcal{T}^{h}=\cup_{i=1}^{N^{m,h}}\varsigma_{i}$) and the fracture domain ($\mathcal{G}^{h}=\cup_{l=1}^{N^{f,h}}\iota_{l}$) with cells $\varsigma_{i}$ and $\iota_{l}$, $N^{m,h}$ is the number of cells in $\mathcal{T}^{h}$, $N^{f,h}$ is the number of cell related to fracture mesh $\mathcal{G}^{h}$. For approximation by time, we use IMPES (implicit pressure explicit saturation) scheme and a finite volume approximation by space

[TABLE]

where $s_{i}^{m}$ and $s_{l}^{f}$ are saturation of matrix and fracture continuum in cells $\varsigma_{i}$ and $\iota_{l}$, $f^{w,\alpha}_{i}={\lambda^{w}(s^{\alpha}_{i})}/{\lambda(s^{\alpha}_{i})}$, $|\varsigma_{i}|$ and $|\iota_{l}|$ are the volume of cells. Here, we use an implicit scheme for time discretization, where $\check{s}^{m}_{i}$ and $\check{s}^{f}_{l}$ are the solutions from previous time step and $\tau$ is the given time step.

For the fluxes

[TABLE]

[TABLE]

[TABLE]

[TABLE]

we have following approximations

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

and $W_{ij}^{mm}=k^{m}_{ij}|E_{ij}|/d_{ij}$, $W_{ln}^{ff}=k^{f}_{ln}/\Delta_{ln}$, $W_{il}^{mf}=\sigma_{il}$ with $\sigma_{il}=CI_{il}k^{m}_{il}$ if $\iota_{l}\subset\varsigma_{i}$ and zero else. Here $|E_{ij}|$ is the length of face between cells $\varsigma_{i}$ and $\varsigma_{j}$, $d_{ij}$ is the distance between midpoint of cells $\varsigma_{i}$ and $\varsigma_{j}$, $\Delta_{ln}$ is the distance between points $l$ and $n$, $CI_{il}$ is the connectivity index from [24, 30, 18] that proportional to the distance and area of the intersection between the fracture cell $\iota_{l}$ and porous matrix cell $\varsigma_{i}$.

Therefore, we have following discrete problem on the fine grid

[TABLE]

where $i=1,...,N^{m,h}$ and $l=1,..,N^{f,h}$.

For approximation of the $\lambda^{w}(s^{\alpha\beta}_{ij})$, we use an upwind scheme

[TABLE]

and $\lambda(s^{\alpha\beta}_{ij})$ is the harmonic average between $\lambda(s^{\alpha}_{i})$ and $\lambda(s^{\beta}_{j})$.

For the general multicontinuum model, we have

[TABLE]

where $\alpha=1,...,M$.

3 Coarse grid upscaled model

Solution of the problems in heterogeneous and fractured media require the construction of the grids that resolve all small scale heterogeneity. One of the widely used way is based on the numerical homogenization technique, where effective parameters are calculated in order to construct a coarse grid approximations. The coarse grid parameters are constructed by a solution of the local problems with appropriate boundary conditions.

Let $\mathcal{T}^{H}$ be a structured coarse mesh of the computational domain $\Omega$

[TABLE]

where $N^{H}$ is the number of the coarse grid cells, $K_{i}$ is the quadrilateral coarse cell and $i$ is the coarse grid cell index [37, 31]. Form of the coarse grid upscaled model is similar to the fine grid model with finite volume approximation, where coarse grid transmissibilities $T_{ij}^{UP}$ are calculated by a solution of the local problems that take into account fine grid resolution of the heterogeneous permeability (see Figure 1).

We let $E_{ij}$ be the coarse grid face and we define the neighborhood (local domain) by

[TABLE]

where $\omega_{ij}$ is a union of two coarse cells, when $E_{ij}$ lies in the interior of the domain $\Omega$. For the edges on the boundary, we will use a no flux boundary conditions and therefore not need to calculate of the upscaled transmissibilities. For calculation of the upscaled transmissibilities $T_{ij}^{UP}$ for coarse face $E_{ij}$, we solve local problems for nonperiodic heterogeneous fractured media with linear boundary conditions in $\omega_{ij}$. For the fractured/multicontinuum media, we use a similar approach for calculation of the coarse grid transmissibilities $T_{ij}^{\alpha\beta,UP}$. Details of the calculations, we present below for each problem.

3.1 Nonlinear flow problem

We start with single-phase upscaling and suppose that

[TABLE]

where, in general, $k_{r}$ and $\sigma_{r}$ can be different for each continuum $\alpha$, but in this work, we assume similar relationships, for simplicity.

Let $\mathcal{T}^{H}$ denote a structured coarse grid of the porous matrix domain $\Omega$ and $\mathcal{G}^{H}$ denote a coarse grid of the fracture domain $\gamma$

[TABLE]

where $K_{i}$ and $\gamma_{l}$ are the cell of the matrix and fractures fine grids, $N^{m,H}$ is the number of coarse cells in $\mathcal{T}^{H}$, $N^{f,H}$ is the number of coarse cells related to $\mathcal{G}^{H}$. On the coarse grid for equation (4), we have the following discrete problem for $\overline{p}=(\overline{p}^{m},\overline{p}^{f})$

[TABLE]

where $l=1,..,N^{f,H}$ and $i=1,...,N^{m,H}$.

In general for multicontinuum model, we have

[TABLE]

where $\overline{c}^{\alpha}_{i}\approx\frac{1}{|K_{i}|}\int_{K_{i}}c^{\alpha}(\overline{p}_{i}^{\alpha})\,dx$ and

[TABLE]

and $W^{\alpha\beta,UP}_{ij}$ is the precalculated effective transmissibilities.

For the calculation of the upscaled transmissibilities for the porous matrix, we solve the following local problems in each $\omega_{ij}$ (see Figure 1, where local domain is depicted by a yellow color)

[TABLE]

with boundary conditions

[TABLE]

In this work, we consider two-dimensional problems with $x=(x_{1},x_{2})$. Therefore, we solve two local problems for $\psi^{l}$, $l=1,2$. For $\psi^{1}$, boundaries $\Gamma_{ij}^{1}$ and $\Gamma_{ij}^{2}$ are the left and right boundaries of the domain $\omega_{ij}$, respectively. For $\psi^{2}$, boundaries $\Gamma_{ij}^{1}$ and $\Gamma_{ij}^{2}$ are the top and bottom boundaries of the domain $\omega_{ij}$, respectively. Note that, another boundary conditions can be applied for local problems.

Therefore, for calculations $W_{ij}^{mm,UP}$ in (14), we solve following discrete problem for finite volume approximation up to fine grid resolution

[TABLE]

with appropriate boundary conditions.

After solution of the local problems in $\omega_{ij}$, we calculate upscaled transmissibility for the porous matrix (see Figure 1, where interface $E_{ij}$ is depicted by a red color)

[TABLE]

where $r,n$ are the fine cells around coarse face $E_{ij}$, $\overline{\psi}^{l}_{i}$ and $\overline{\psi}^{l}_{j}$ are the mean values in coarse cells $K_{i}$ and $K_{j}$. We use $\psi^{l}$ with $l=1$ for all vertical edges and $l=2$ for horizontal edges. For the fracture continuum, we suppose that $k^{f}=\mathop{\rm const}\nolimits$ and therefore set $W_{ln}^{ff,UP}=k^{f}/d_{ln}$ ($d_{ln}$ is the distance between points $l$ and $n$).

In this work, we suppose $k^{f}=const$, and therefore for the calculations of the coarse grid transmissibility between coarse grid fracture cells $\gamma_{l}$ and $\gamma_{n}$, we have $T^{ff,UP}_{ln}=k^{f}/\Delta_{ln}$, where $\Delta_{ln}$ is the distance between midpoint of cells $\gamma_{l}$ and $\gamma_{n}$.

Let $\omega^{mf}_{il}=\{K_{i}:K_{i}\cup\gamma_{l}\neq\emptyset,\,\gamma_{l}\in\mathcal{G}^{H},\,K_{i}\in\mathcal{T}^{H}\}$ be the local domain for calculation of the $W^{mf}_{il}$ in (14) (see Figure 1, where local domain is depicted by a orange color). For the calculation of the upscaled transmissibility between porous matrix and fracture, we solve local problems in $\omega^{mf}_{il}$

[TABLE]

where $\phi^{f}=1$ on $\gamma_{l}$ with zero flux boundary conditions on $\partial\omega^{mf}_{il}$. Therefore, we solve following discrete system for finite volume approximation up to fine grid resolution

[TABLE]

until $|\phi_{i}-\check{\phi}_{i}|>\epsilon$ and find upscaled matrix-fracture transmissibility using final time step solution

[TABLE]

where $r$ are the cell that contains fracture, $\overline{\phi}_{i}$ and $\overline{\phi}^{f}_{l}$ are the mean values in coarse cells $K_{i}$ and in fracture $\gamma_{l}$ (see Figure 1, where interface $\gamma_{l}$ is depicted by a red color).

Note that, there exist different approaches for calculation of the effective transmissibilities, for example, based on the different boundary conditions for local problems, using oversampled domains and the construction of look-up table for interpolation of the nonlinear dependence. In this work, for calculating the upscaled transmissibilities, we use the simplest classic approach. The main goal of the paper is the construction of the novel highly accurate nonlinear upscaled coarse grid approximations using machine learning techniques.

3.2 Nonlinear flow and transport problem

On the coarse grid for equation (8), we have following discrete problem for $\overline{p}^{m},\overline{p}^{f},\overline{s}^{m}$ and $\overline{s}^{f}$

[TABLE]

where $l=1,..,N^{f,H}$ and $i=1,...,N^{m,H}$. For approximation by time similarly to the fine grid approximation, the IMPES scheme is used.

For the general multicontinuum model, we have

[TABLE]

where

[TABLE]

with upwind scheme approximation of $\lambda^{w}$ and $W^{\alpha\beta,UP}_{ij}$ is the precalculated effective transmissibilities that is similar to the previous problem and based on the single phase upscaling.

The choice of boundary conditions have a strong impact on the accuracy of results. In the presented standard upscaling method, coarse grid parameters are obtained independently to global problem solution information. More accurate approaches can be based on the information about the fine scale flow in the local domains up to fine grid resolution and without variable separation of nonlinear coefficients. For example, an interpolated global coarse grid solution is used for performing accurate construction of the upscaled transmissibilities in [5], which involve iterations between global coarse grid model and local fine grid calculations with updating of the upscaled transmissibilities. The local-global upscaling method requires extra computations than existing classic upscaling procedures.

In this work, the construction of the accurate upscaled transmissibilities for the coarse grid approximation is also based on the information about global solution (nonlinear transmissibilities). Moreover, the presented method is based on the machine learning procedure for fast prediction of the nonlinear transmissibilities, where we construct neural network that learn dependencies between the coarse grid quantities on the oversampled local domains and upscaled transmissibilities. We use a convolutional neural network and GPU training process to construct a machine learning algorithm.

4 Machine learning for nonlinear nonlocal upscaled transmissibilities

We consider a machine learning approach for prediction of the upscaled nonlinear nonlocal transmissibilities for accurate and fast coarse grid approximation. We have following main steps:

Generate dataset to train, validate and test of the neural network. 2. 2.

Neural networks training, validation and testing. 3. 3.

Calculation of the nonlinear upscaled transmissibilities on the fly using constructed neural networks during coarse system construction, fast and accurate solution of the upscaled system.

For construction of the datasets, we perform local or global calculations of the coarse grid quantities [5]. In local approach, the upscaled transmissibilities are calculated on the local domain corresponding to the target face, where the fine-scale solution information is used to set boundary conditions. Global approach uses a global fine-scale solution for the determination of coarse scale parameters. For training of the neural networks, we use a family of problem solutions for different input conditions (snapshots). Note that, we should have many snapshots to capture all input condition variations because the accuracy of the machine learning method depends on snapshot space that is as train dataset. Next, we consider dataset generation and network construction in detail.

4.1 Dataset

The most accurate case can be based on the fine grid solution, at the same time for upscaled model, we would like to use only coarse-grid information. For possible applicability of this, we construct a novel coarse grid model, using machine learning algorithms and construct neural network that learn dependency between coarse grid functions $\overline{p}^{\alpha}$ in local domains and upscaled nonlinear transmissibilities.

For constructing accurate neural network for prediction of the transmissibilities, we should train network on the highly accurate dataset. One of the most accurate approach for calculating upscaled transmissibilities based on the direct calculation from the fine scale solution. We use following coarse grid approximation (similar to previous section)

• •

Unsaturated flow problem (nonlinear flow):

[TABLE]

with nonlinear upscaled transmissibilities

[TABLE]

• •

Two-phase flow problem (nonlinear transport and flow):

[TABLE]

with nonlinear upscaled transmissibilities

[TABLE]

Because ${T}^{w,\alpha\beta,NL}$ and ${T}^{\alpha\beta,NL}$ ($\alpha,\beta=m,f$) are nonlinear and depend on the fine grid solutions $p^{\alpha},p^{\beta},s^{\alpha},s^{\beta}$, we cannot directly use such transmissibilities on the coarse grid model. For possible applicability of this, we will use a machine learning algorithms and construct a neural network that learn dependence between coarse grid functions $\overline{p}^{\alpha},\overline{p}^{\beta},\overline{s}^{\alpha},\overline{s}^{\beta}$ in local domains (oversampled) and upscaled nonlinear transmissibilities.

Let $E_{l}$ is the interface, where we define upscaled transmissibility, and $X_{l}$ and $Y_{l}$ are the input data and output data for machine learning algorithm and

[TABLE]

For constructing neural network for upscaled transmissibilities, based on the (23) and (25), we use a nonlocal upscaled transmissibilities ${T}_{l}^{w,\alpha\beta,NL}$ and ${T}_{l}^{\alpha\beta,NL}$ ($\alpha,\beta=m,f$) as output data $Y_{l}$. Input data $X_{l}$ is contains information about fine scale permeabilities, fracture position in local domain $\omega_{E_{l}}$, coarse grid functions $\overline{p}$ and $\overline{s}$ in the oversampled local domains. For this purpose, we use a local multi-input data for training neural network

[TABLE]

where $X_{l}^{k}$ and $X_{l}^{f}$ are the local heterogeneous permeabilities and local fracture position markers in local domain $\omega_{E_{l}}$; $X_{l+}^{\overline{p}^{\alpha}}$ and $X_{l+}^{\overline{s}^{\alpha}}$ are the coarse grid nonlocal mean values for pressure and saturation for continuum $\alpha$ in oversampled local domain $\omega^{H+}_{E_{l}}$. Each of the input fields is represented as two-dimensional array for two-dimensional problem. The scale of each array in dataset is re-scaled to fall within the range [math] to $1$.

In Figure 2, we present an illustration of the local domains $\omega_{E_{l}}$ and $\omega^{H+}_{E_{l}}$, $l=1,...,N_{E}$ ($N_{E}$ – number local domains). Local domain $\omega_{E_{l}}$ is the domain for edge $E_{l}$ up to fine grid resolution that is similar to the classic (single phase) upscaling presented in previous section. In local domain $\omega_{E_{l}}$, we define $X_{l}^{k}\in\omega_{E_{l}}$ and $X_{l}^{f}\in\omega_{E_{l}}$. Oversampled local domain $\omega^{H+}_{E_{l}}$ is the domain around $E_{l}$ up to coarse grid resolution, where we define $X_{l+}^{\overline{p}^{\alpha}}\in\omega^{H+}_{E_{l}}$ and $X_{l+}^{\overline{s}^{\alpha}}\in\omega^{H+}_{E_{l}}$. To ensure same size and structure of the input data, we divide all local data into four types: matrix-matrix flow through horizontal edge ($T_{l}^{mm,NL}$), matrix-matrix flow through vertical edge ($T_{l}^{mm,NL}$, fracture -matrix flow for $T_{l}^{mf,NL}$ and fracture-fracture flow for $T_{l}^{ff,NL}$.

The output is the normalized array of the upscaled transmissibilities

[TABLE]

for edge $E_{l}$. Dataset is divided into train, validation and test sets with sizes $N_{train}$, $N_{val}$ and $N_{test}$ ($N=N_{train}+N_{val}+N_{test}$). For each type of local domain as a test set, we take 50 % of data, another 50 % divided between train and validation set in 80/20 proportion.

We use dataset $(X_{l},Y_{l})$ for training of the neural network. To ensure a good learning rate and for obtaining a wide coverage of data, we generate several solution snapshots by varying of the source term in global fine grid model. Another approach is related to the localization of the dataset generation, where we can use local domains calculations for fine grid solutions and calculations of the $X_{l}$ and $Y_{l}$. Note that, this machine learning approach for the learning of the nonlocal nonlinear upscaled transmissibilities can be also applied for the linear problems and has a recap with nonlinear finite volume methods, where transmissibilities are also depends on solution.

4.2 Network

In machine learning algorithm, we use a multi-input deep neural network (convolutional neural network). Let

[TABLE]

where $X_{l}=(X_{l}^{1},...,X_{l}^{s})$ ($s$ is the number of the input data for $E_{l}$, see (26)). Each input data $X_{l}^{i}$ is defined in $\omega^{i}_{E_{l}}$ and represented as two-dimensional array for two-dimensional problems. The architecture of the multi-input deep neural network for prediction of the nonlinear nonlocal upscaled transmissibilities is presented in Figure 3. For each input data $X^{i}$, we use a convolutional neural network [23, 22]. Several convolutional and pooling layers with rectified linear units activation layer are stacked with a several fully-connected layers with dropout. Several layers of convolutions and pooling are alternated in order to detect higher order features for better accuracy of the method. After performing convolutions, pooling, activation and dropout layers for each $X^{i}$ ($i=1,...,s$), we add a fully connected layers, where we compose all outputs on CNN together. By a training process, a machine learning algorithm solve the optimization problem to find model weights that best describe the train set by minimization of the loss function.

We train a convolutional neural network by a dataset of local multi-input data ($X_{l}$) and upscaled transmissibilities ($Y_{l}$). As a loss function, we use the mean square error (MSE)

[TABLE]

For solution of the minimization problem, we use gradient-based optimizer Adam [21]. Implementation of the machine learning method is based on the open source library Keras [6] with TensorFlow backend [1] and performed on the GPU. Constructed machine learning algorithm will efficiently determine dependence between coarse grid functions in local domains and upscaled transmissibilities.

5 Numerical result

In this section, we present numerical results for the proposed method. We consider following model problems in fractured and heterogeneous porous media:

• Test 1: Nonlinear flow problem (unsaturated flow problem)
• Test 2: Nonlinear transport and flow problem (two-phase flow problem)

We solve model problem in $\Omega=[0,1]\times[0,1]$ with no flux boundary conditions. We use $10\times 10$ coarse grid. Location of the source terms and fracture position are depicted in Figure 4. In Figure 4, we show a heterogeneous porous matrix permeability for both test problems. The numerical calculations of the effective properties has been implemented with the open-source finite element software PETSc and FEniCS [25, 26, 3].

To measure difference between reference solution and coarse grid solution, we compute relative $L_{2}$ error

[TABLE]

where $u=p,s$, $\overline{u}^{fine}$ is the reference solution (mean value on coarse grid of the fine grid solution) and $\overline{u}$ is the solution on the coarse grid.

For each test problem, we present results of the fine scale solution, for upscaling technique presented in Section 3 and new method from Section 4. Computational algorithm for single-phase upscaling method with $T^{\alpha\beta,UP}$ (Section 3):

Loading of the precalculated effective transmissibilities $W^{\alpha\beta,UP}$. 2. 2.

Solution of the multicontinuum model:

• Test 1: Nonlinear flow problem with

[TABLE]

• Test 2: Nonlinear transport and flow problem with

[TABLE]

[TABLE]

with upwind approximation of $\lambda^{w}$ on the coarse grid.

For the new nonlocal nonlinear machine learning technique with $T^{\alpha\beta,NL}$ (Section 4), we have:

Loading of the machine learning models, $NN_{i}$ ($i=1,2,...$). 2. 2.

Solution of the multicontinuum model:

• Test 1: Nonlinear flow problem with

[TABLE]

where ${T}_{ij}^{\alpha\beta,ML}$ is the value predicted using machine learning algorithm.

• Test 2: Nonlinear transport and flow problem with

[TABLE]

[TABLE]

where ${T}_{ij}^{\alpha\beta,ML}$ and ${T}_{ij}^{w,\alpha\beta,ML}$ is the value predicted using machine learning algorithm.

Note that, the loss of positivity of upscaled transmissibilities can happen, and we use a threshold value $\varepsilon$ for the pressure difference to guarantee a good values of the coarse grid parameters, where linear upscaling is used for the faces with small pressure difference. Moreover, we used predicted transmissibilities adaptively with parameter $\varepsilon_{s}$ in the coarse grid model with machine learning approach.

We will show results of the learning process of deep neural network for nonlocal nonlinear upscaled transmissibilities and calculate errors for a given datasets. Finally, we consider a coarse grid solution of the problem, where nonlocal nonlinear upscaled transmissibilities are calculated using constructed machine learning method. Finally, we discuss the computational time of the neural networks construction and solution of the coarse grid system using classic upscaling and machine learning approaches. We divide calculation on the offline and online stages. On the online stage, we train neural network on the GPU by a given train and validation datasets. On the offline stage, we have two steps: loading of the preconstructed neural network and prediction of the upscaled coarse grid transmissibilities on each time iteration or/and nonlinear iteration.

5.1 Nonlinear flow problem

We consider the solution of the nonlinear equation in fractured and heterogeneous porous media. We set source terms $f^{\pm}=\pm q$, $q=10^{5}$. For the nonlinear coefficient, we use $k^{\alpha\beta}(x,u)=k_{s}(x)k_{r}(u)$ with $k_{r}(u)=\exp(-a|u|)$, $a=0.1$ ($\alpha,\beta=m,f$). In Figure 4 (second column), we show a heterogeneous porous matrix permeability $k^{m}_{s}(x)$ and fracture position. We set $c^{m}=1$, $c^{f}=0$, $k_{s}^{f}=10^{6}$ and $T_{max}=10^{-3}$ with 20 time steps. Coarse grid is $10\times 10$ and fine grid is $640\times 640$ for domain $\Omega$.

We present results for the machine learning algorithm and calculate errors for train and test datasets. For the training of the neural networks, we investigate two datasets: local and global. For the global dataset, we extract local information from the fine grid calculations on the global domain $\Omega$. For the local dataset, we calculate each data by solution of the local problem up to fine grid resolution with different boundary conditions for generation of the possible set of solutions (snapshots). We use six random values of the source term to generate datasets ($N_{r}=6$). We train three neural networks for each type of transmissibility: $NN_{1}$ for horizontal coarse edges for matrix-matrix flow, $NN_{2}$ for vertical coarse edges s for matrix-matrix flow and $NN_{3}$ for matrix - fracture flow. For $10\times 10$ coarse mesh, we have $N_{E}=90$ horizontal and $N_{E}=90$ vertical coarse edges (without boundary edges due to no flux boundary conditions), furthermore, we have $N_{E}=5$ coarse cells with fracture.Therefore, the train dataset for neural network contains $N=N_{r}\cdot N_{E}\cdot N_{t}$ samples for learning process, where $N_{t}$ is the number of time steps. We have $N=10800$ for $NN_{1}$ and $NN_{2}$; and $N=600$ for $NN_{3}$. Each sample $X_{l}$ contains information about heterogeneous permeability and fracture position up to fine grid resolution in local domain, coarse grid mean value of the solution in oversampled local domain

[TABLE]

Each dataset is divided into training and validation sets with $80:20$ ratio. For testing, we calculate another six solution snapshots.

For calculations, we use 500 epochs with a batch size $N_{b}=90$ and Adam optimizer with learning rate $\epsilon=0.001$. For accelerating of the training process of the multi-input CNN, we use GPU. We use $3\times 3$ convolutions and $2\times 2$ maxpooling layers with RELU activation for $X^{k}$ and $X^{f}$, and $3\times 3$ convolutions with RELU activation for $X^{\overline{p}^{m}}$. For each input data, we have 2 layers of CNN with one final fully connected layer. Convolution layer contains 8 and 16 feature maps for $X^{k}$ and $X^{f}$; and 4 and 8 feature maps for $X^{\overline{p}^{m}}$. We use dropout with rate 10 % in each layer in order to prevent over-fitting. Finally, we combine CNN output and perform two additional fully connected layers with size 200 and 1(one final output). Presented algorithm is used to learn dependence between multi-input data and upscaled nonlinear transmissibilities.

For error calculation on the train and test dataset, we use mean square errors, relative mean absolute and relative root mean square errors

[TABLE]

where $Y_{i}$ and $\tilde{Y}_{i}$ denotes reference and predicted values for sample $X_{i}$

Convergence of the loss functions for three neural networks for Test 1 are presented in Figure 5, where we plot the MSE loss function vs epoch number for train and validation sets. In Figure 6, we present a parity plots comparing reference values against predicted using trained neural networks for train and test datasets (green and blue colors). Learning performance for neural networks are presented in Tables 7 for global and local datasets. We observe good convergence of the relative errors for train and test sets with $1-2\%$ of RMSE.

Next, we consider errors between solution of the coarse grid problem with reference and predicted upscaled transmissibilities. In Figure 7, we present results for 50 test problems with random value of the source term. We show a relative $L_{2}$ errors for pressure head on the coarse mesh with classic upscaling algorithm and using new nonlocal nonlinear transmissibilities. We observe small errors ($1-2\%$) for predicted nonlocal nonlinear transmissibilities compared with classical upscaling technique, where we have $\approx 15\%$ of relative $L_{2}$ errors for pressure head. Furthermore, we see that local calculation of the dataset provide similar results as a globally caclulated data.

In Figure 8, we depict solution of the problem on the fine grid, coarse grid upscaled solution using classic approach from Section 3 and for new method presented in Section 4 ($u^{fine}$, $\overline{u}^{fine}$, $\overline{u}^{UP}$ and $\overline{u}^{NL}$). For $\overline{u}^{UP}$, we apply presented upscaling method 16, 18 and 12. We have $e(\overline{u}^{UP})=14.772\%$ and $e(\overline{u}^{NL})=1.463\%$ at final time. For the nonlinear nonlocal transmissibilities, we set $\varepsilon=0.5\cdot 10^{-1}$ for $NN_{1}$ and $NN_{2}$, $\varepsilon=10^{-20}$ for $NN_{3}$. Note that, we didn’t construct $NN_{4}$ of data ($T^{ff,NL}$) because for our test problem we observe almost constant pressure on the fracture and set $T^{ff,NL}=T^{ff,UP}$ on the coarse grid.

We perform training of the neural networks on the GPU, where we train three neural networks: $NN_{1}$, $NN_{2}$ and $NN_{3}$. Online stage (neural network training) time is 25 minutes for $NN_{1}$, 28 minutes for $NN_{2}$ and 6 minutes for $NN_{3}$ on GPU (GeForce GTX 1060). Note that, the training time depends on size of the dataset and GPU model. Here we didn’t consider time of the dataset construction which depends on type of calculations (global or local) and number of solution snapshots, that we used for training. Number of snapshots ($N_{r}$) is also effects to the algorithm errors because we should have sufficient number of snapshots to capture all variations of the input data to know how it effects to the output.

Time of the online stage contains 6.6 seconds of loading three neural networks and 13.0 seconds for calculations on the $10\times 10$ coarse grid with prediction of the nonlinear nonlocal transmissibilities. Fine grid calculations time is 454 seconds for 20 time steps on $640\times 640$ fine grid. We have approximately 35 time faster calculations for a new method with small error of the coarse grid solution.

5.2 Nonlinear flow and transport problem

We consider solution of the two-phase flow problem in fractured and heterogeneous porous media. For nonlinear coefficient, we set $\lambda^{w}(s)=s^{2}$ and $\lambda^{n}(s)=(1-s)^{2}$. In Figure 4 (third column), we show the heterogeneous porous matrix permeability $k^{m}(x)$ and fracture position. We set $\phi^{\alpha}=1$ ($\alpha=m,f$), $k^{f}=10^{3}$ and $T_{max}=25\cdot 10^{-3}$ with 250 time steps. Coarse grid is $10\times 10$ and fine grid is $160\times 160$ for domain $\Omega$.

For the training of the neural networks, we use a global dataset, where we extract local information from the fine grid calculations on the global domain $\Omega$. For generation of the train datasets, we use a three random shapshots ($N_{r}=3$) with $T_{max}=40\cdot 10^{-3}$ and 400 time steps. We train four neural networks for each type of transmissibility: $NN_{1}$ for horizontal coarse edges for matrix-matrix flow, $NN_{2}$ for vertical coarse edges s for matrix-matrix flow, $NN_{3}$ for matrix - fracture flow and $NN_{4}$ for fracture - fracture flow. The train dataset for first and second neural networks contains $N=108000$; $N=6000$ for $NN_{3}$ and $N=4800$ for $NN_{4}$, where dataset is randomly divided into training and validation sets with $80:20$ ratio. Each sample $X_{l}$ contains information about heterogeneous permeability and fracture position up to fine grid resolution in local domain, mean value of the solution in oversampled local domain (coarse grid)

[TABLE]

and output

[TABLE]

For calculations, we use 150 epochs with a batch size $N_{b}=90$ and perform calculations on GPU. Architecture of the neural networks are similar to the previous test problem but as output for this case, we obtain two values, $T$. Learning performance for neural networks are presented in Tables 9 and 2 for train datasets. We observe a good convergence with small error for each neural network.

In Figure 10, we present results for 20 test problems with random value of the source terms. We show a relative mean square error in percentages for pressure and for saturation on the coarse mesh with classic upscaling algorithm and using new nonlocal nonlinear transmissibilities.

In Figure 11, we depict solution of the problem using different methods. On the first conlumn, we depict a reference fine grid solution ($s^{fine}$, $p^{fine}$), mean value on coarse grid of the fine grid solution ($\overline{s}^{fine}$, $\overline{p}^{fine}$) on the second column, coarse grid solution using upscaling method ($\overline{s}^{UP}$, $\overline{p}^{UP}$) on the third column and coarse grid solution using nonlinear nonlocal machine learning method ($\overline{s}^{NL}$, $\overline{p}^{NL}$) on the fourth column. On the first, second and third rows, we show a saturation for time $t_{m}$, $m=50,150,250$ and on fourth row, we have pressure for time $t_{m}$, $m=250$. For solution on the coarse grid ($\overline{p}^{UP}$ and $\overline{s}^{UP}$), we applied classic upscaling method (see Section 3). Fine grid (reference) solution is performed using finite volume approximation with embedded discrete fracture model, where for error calculations we used a mean values of the reference solution on the coarse grid, $\overline{p}^{fine}$ and $\overline{s}^{fine}$. On the last column of the Figure 11, we depict a coarse grid solution using nonlinear nonlocal transmissibilities that calculate based on the machine learning approach. For machine learning approach, we have $e(\overline{p}^{NL})=0.920\%$, $e(\overline{s}^{NL})=3.957\%$, and for upscaling $e(\overline{p}^{UP})=13.519\%$, $e(\overline{s}^{UP})=13.227\%$ at final time $t_{m}$, $m=250$. For the nonlinear nonlocal transmissibilities, we set $\varepsilon_{S}=10^{-2}$ for transport and $\varepsilon=0.5\cdot 10^{-2}$ for $NN_{1}$, $\varepsilon=10^{-4}$ for $NN_{2}$, $\varepsilon=10^{-3}$ for $NN_{3}$ and $\varepsilon=10^{-20}$ for $NN_{4}$ for flow.

We perform training of the neural networks on the GPU, where we train four neural networks: $NN_{1}$, $NN_{2}$, $NN_{3}$ and $NN_{4}$. Online stage (neural network training) time is 80 minutes for $NN_{1}$, 59 minutes for $NN_{2}$, 2 minutes for $NN_{3}$ and 4 minutes for $NN_{4}$ on GPU (GeForce GTX 1060). Note that, the training time depends on size of the dataset and GPU model. Time of the online stage contains 16.7 seconds of loading four neural networks and 46.9 seconds for calculations on the $10\times 10$ coarse grid with prediction of the nonlinear nonlocal transmissibilities. Fine grid calculations time is 812 seconds for 250 time steps on $160\times 160$ fine grid for transport and flow model. We observe a good results with fast calculations using a machine learning algorithm for presented method.

6 Conclusion

In this work, we consider two nonlinear problems in heterogeneous and fractured porous media. Mathematical models are formulated as a general multicontinuum models, where fine grid approximations are constructed using finite volume method. For the accurate solution of the nonlinear problems on the coarse grid, a novel machine learning algorithm combined with nonlinear nonlocal multicontinua approach for calculating nonlocal nonlinear transmissibilities is presented and investigated. We presented the construction of the dataset for training deep neural networks. The construction of the neural network is based on the multi-input convolutional neural networks, where GPU is used for performing a fast learning process. To illustrate the applicability of the presented method, we presented numerical results for two test problems. Numerical results showed that presented algorithm provides fast and accurate calculations of the nonlocal nonlinear transmissibilities.

7 Acknowledgements

MV’s work is supported by the mega-grant of the Russian Federation Government N14.Y26.31.0013 and RSF N17-71-20055. EC’s work is partially supported by the Hong Kong RGC General Research Fund (Project numbers 14304217 and 14302018) and the CUHK Faculty of Science Direct Grant 2018-19.