## Generalized Multiscale Finite Element method for multicontinua unsaturated flow problems in fractured porous media

Denis Spiridonov ∗

Maria Vasilyeva †

Eric T. Chung ‡

August 2, 2019

## Abstract

In this paper, we present a multiscale method for simulations of the multicontinua unsaturated flow problems in heterogeneous fractured porous media. The mathematical model is described by the system of Richards equations for each continuum that coupled by the specific transfer term. To illustrate the idea of our approach, we consider a dual continua background model with discrete fractures networks that generalized as a multicontinua model for unsaturated fluid flow in the complex heterogeneous porous media. We present fine grid approximation based on the finite element method and Discrete Fracture Model (DFM) approach. In this model, we construct an unstructured fine grid that take into account a complex fracture geometries for two and three dimensional formulations. Due to construction of the unstructured grid, the fine grid approximation leads to the very large system of equations. For reduction of the discrete system size, we develop a multiscale method for coarse grid approximation of the coupled problem using Generalized Multiscale Finite Element Method (GMsFEM). In this method, we construct a coupled multiscale basis functions that used to construct highly accurate coarse grid approximation. The multiscale method allowed us to capture detailed interactions between multiple continua. The adaptive approach is investigated, where we consider two approaches for multiscale basis functions construction: (1) based on the spectral characteristics of the local problems and (2) using simplified multiscale basis functions. We investigate accuracy of the proposed method for the several test problems in two and three dimensional formulations. We present a comparison of the relative error for different number of basis functions and for adaptive approach. Numerical results illustrate that the presented method provide accurate solution of the unsaturated multicontinua problem on the coarse grid with huge reduction of the discrete system size.

∗ Multiscale model reduction laboratory, North-Eastern Federal University, Yakutsk, Republic of Sakha (Yakutia), Russia, 677980. Email: d.stalnov@mail.ru .

† Institute for Scientific Computation, Texas A&amp;M University, College Station, TX 77843-3368 &amp; Department of Computational Technologies, North-Eastern Federal University, Yakutsk, Republic of Sakha (Yakutia), Russia, 677980. Email: vasilyevadotmdotv@gmail.com .

‡ Department of Mathematics, The Chinese University of Hong Kong (CUHK), Hong Kong SAR. Email: tschung@math.cuhk.edu.hk .

## 1 Introduction

Prediction of fluid flow in unsaturated soils is an important problem in science and engineering, for example, in agriculture and environmental engineering, oil and gas production, groundwater hydrology etc. Mathematical models of the unsaturated filtration processes are described by the Richards equation [1, 2, 3], where porous media characterized by a complex heterogeneous rock properties, complex fracture distribution, multiple scales and high contrast of the properties. Due to the high permeability of the fractures, they have a significant effect on the flow processes and requires a special approach in the construction of a mathematical model and computational algorithms. For example, a hierarchical model is used to describe the multiscale fractures [4, 5, 6, 7], where flow in the small scale highly connected fracture networks (natural fractures) can be described by the dual continua approach with additional lower dimensional discrete fracture for simulation of the flow in large scale fractures.

In this work, we construct a multicontinua models for unsaturated flow in porous media, that based on the coupled system of Richards equations. We consider fractured porous media, where one of the continua is describes a flow in the fracture networks. For accurate simulations in the fractured and heterogeneous porous media, the unstructured grids should be used to resolve features in the level of mesh construction. Such grids lead to the large number of unknowns of the discrete system and computationally expensive. To solve such problems, a homogenization techniques or multiscale methods are used to reduce size of the discrete systems.

Several multiscale methods are developed to solve problems with heterogeneous properties, for example, multiscale finite element method (MsFEM) [8], multiscale finite volume method (MsFVM) [9], heterogeneous multiscale methods (HMM) [10], generalized multiscale finite element method (GMsFEM) [11], constraint energy minimizing generalized multiscale finite element method (CEM-GMsFEM) [12]. In [13, 14, 15], the authors present construction of the coarse grid approximation based on the MsFEM for solving the unsaturated flow problems with heterogeneous coefficients. Upscaling method for the Richards equation is presented in [16]. Multiscale methods for solution of the flow problems in fractured porous media are presented in [17, 18, 19, 20, 21]. In our previous works, we developed multiscale model reduction techniques based on the Generalized Multiscale Finite Element Method (GMsFEM) for flow in fractured porous media [17, 6, 22]. Recently in [23, 24, 25], we introduced a nonlocal multicontinuum (NLMC) method for problems in fractured porous media, where we construct multiscale basis functions based on the solution of some local constrained energy minimization problems as in the CEM-GMsFEM.

In this work, we present the Generalized Multiscale Finite Element method [26, 27] for solution of the unsaturated multicontinua flow problem in fractured heterogeneous porous media. Previously in [6], we presented construction of the GMsFEM for solution of the saturated flow problem in fractured media, in presented work we extend proposed framework for unsaturated flow case problem in two and three dimensional formulations. In GMsFEM, we construct multiscale basis functions that automatically identify each continuum via solution of the local spectral problems [28, 29, 30]. For coupled system of equations in multicontinua models, we solve local coupled system of equations to construct highly accurate basis functions that capture complex interaction between continua in the coarse grid model. We study an adaptive approach with simplified basis functions construction and compare with adaptive choosing based on the local eigen-

values of the spectral problem. The simplified approach is effective for the cases, when the fracture networks are known and have simplified geometries. To illustrate robustness and accuracy of the proposed method, we consider a dual continuum background model that coupled with the discrete fracture networks (triple continua model). We present several numerical examples for two and three dimensional test problems.

We organize the paper as follows. In Section 2, we present a mathematical model for unsaturated multicontinua flow problem in fractured porous media as a general multicontinia formulation of the unsaturated flow. In section 3, we consider fine-scale approximation using Discrete Fracture Model. Next in Section 4, we describe the construction of the coarse grid approximation and describe construction of the multiscale basis functions. Furthermore, we present an adaptive approach with simplified multiscale basis functions and using adaptive approach based on the spectral characteristics of the spectral problems. We present numerical results in Sections 5 for two and three dimensional model problems. Finally, we present Conclusion.

## 2 Mathematical model

Mathematical model of the unsaturated flow in porous media described by the Richard's equations. For fractured porous media, we consider a mixed dimensional formulation of the flow problem [31, 32, 33, 34]. Let Ω ∈ 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 γ ∈ 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 following coupled system of equations

$$\dots \in \Sigma \text{wasm} \\ \frac { \partial \Theta _ { m } } { \partial t } + \nabla \cdot q _ { m } + L _ { m f } = f _ { m }, \quad x \in \Omega, \\ q _ { m } = - k _ { m } ( x, p _ { m } ) \nabla ( p _ { m } + z ), \quad x \in \Omega, \\ \frac { \partial \Theta _ { f } } { \partial t } + \nabla \cdot q _ { f } - L _ { f m } = f _ { f }, \quad x \in \gamma, \\ q _ { f } = - k _ { f } ( x, p _ { f } ) \nabla ( p _ { f } + z ), \quad x \in \gamma, \\ \cdot \ \sim \quad \cdot \ \dots \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdots$$

where q m and q f are the Darcy velocities in matrix and fractures; 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; Θ m and Θ f are the water content for matrix and fracture; and f m and f f refer to source and sink terms. Here L mf and L fm are the transfer term between matrix-fracture and fracture-matrix, L mf = η m mf q , L fm = η f q mf , q mf ≈ σ mf ( x, p m , p f )( p m -p f ) and ∫ L mf dx -∫ L fm ds = 0.

## Ω γ

After substitution the Darcy's Law into the mass conservation equation, we obtain following system of equations for p m and p f :

$$& \frac { \partial \Theta _ { m } } { \partial t } - \nabla \cdot ( k _ { m } ( x, p _ { m } ) \nabla ( p _ { m } + z ) ) + L _ { m f } = f _ { m }, \quad x \in \Omega, \\ & \frac { \partial \Theta _ { f } } { \partial t } - \nabla \cdot ( k _ { f } ( x, p _ { f } ) \nabla ( p _ { f } + z ) ) - L _ { f m } = f _ { f }, \quad x \in \gamma.$$

Similarly, for the dual continuum background model with lower-dimensional fracture model, we have

following system of equations for p 1 , p 2 and p f :

$$\cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \\ \quad \frac { \partial \Theta _ { 1 } } { \partial t } - \nabla \cdot ( k _ { 1 } ( x, p _ { 1 } ) \nabla ( p _ { 1 } + z ) ) + L _ { 1 2 } + L _ { 1 f } = f _ { 1 }, \quad x \in \Omega, \\ \frac { \partial \Theta _ { 2 } } { \partial t } - \nabla \cdot ( k _ { 2 } ( x, p _ { 2 } ) \nabla ( p _ { 2 } + z ) ) - L _ { 2 1 } + L _ { 2 f } = f _ { 2 }, \quad x \in \Omega, \\ \frac { \partial \Theta _ { f } } { \partial t } - \nabla \cdot ( k _ { f } ( x, p _ { f } ) \nabla ( p _ { f } + z ) ) - L _ { f 1 } - L _ { f 2 } = f _ { f }, \quad x \in \gamma. \\ \text{kindoc} 1 9 \text{ and } \text{$row$ of $not$ and $condition$ and $column$ and $column$ for in allowed to $two$}$$

where subindices 1 2 are , related for first and second continua, and subindex f is related to the lowerdimensional fracture model. System of equations is coupled by the mass exchange terms between continuum and ∫ Ω L αf dx -∫ γ L fα ds = 0, ∫ Ω L αβ dx -∫ Ω L βα dx = 0 for α, β = 1 2 and suppose following linear relation , L 12 = L 21 ≈ σ 12 ( x, p 1 , p 2 )( p 1 -p 2 ).

We consider system of equations (3) unsaturated flow model with following initial conditions

$$p _ { 1 } = p _ { 2 } = p _ { f } = p ^ { 0 },$$

and boundary conditions

$$p _ { 1 } = p _ { 2 } = p _ { f } & = g, \quad x \in \Gamma _ { D }, \\ q _ { 1 } \cdot n = q _ { 2 } \cdot n & = q _ { f } \cdot n = 0, \quad x \in \Gamma _ { N },$$

where ∂ Ω = Γ D ∪ Γ N , Γ D is the top boundary of the computational domain and Γ N denotes left, right and bottom boundaries of the Ω. On the fracture boundaries that inside domain Ω, we set a zero Neumann boundary conditions.

In general, we have following multicontinuum model:

glyph[negationslash]

$$\frac { \partial \Theta _ { \alpha } } { \partial t } - \nabla \cdot ( k _ { \alpha } ( x, p _ { \alpha } ) \nabla ( p _ { \alpha } + z ) ) + \sum _ { \beta \neq \alpha } L _ { \alpha \beta } = f _ { \alpha },$$

where α = 1 , ..., L and L is the number of continuum.

As constitutive relations, we use Haverkamp model [3]

$$\Theta _ { \alpha } ( p _ { \alpha } ) = \frac { A _ { \alpha } ( \Theta _ { \alpha, s } - \Theta _ { \alpha, r } ) } { A _ { \alpha } + | p _ { \alpha } | ^ { B _ { \alpha } } } + \Theta _ { \alpha, r }, \ \ k _ { \alpha } ( x, p _ { \alpha } ) = k _ { \alpha, s } ( x ) \frac { C _ { \alpha } } { C _ { \alpha } + | p _ { \alpha } | ^ { D _ { \alpha } } },$$

where A α B α , C α , D α , Θ α,r and Θ α,s are the Haverkamp model coefficients, k α,s ( x ) is the heterogeneous saturated hydraulic conductivity for α continuum.

## 3 Fine grid approximation

For numerical solution, we construct unstructured fine grid that explicitly resolve fractures in the level of mesh. We construct a discrete system based on the finite element method and Discrete Fracture Model (DFM) for fracture networks. Let T h denote a finite element partition of the domain Ω ∈ R d and T h = ⋃ i K i , where K i is the triangular element for d = 2 and tetrahedron for d = 3. Let E h the set of all the faces between the elements T h and E γ ⊂ E h be the subset of all faces that represent fractures. Furthermore, fracture facets

E γ represent the lower dimensional fracture grid. Then finite element approximation for system of equations (3) with boundary conditions (5) can be written as follows: find p = ( p , p 1 2 , p f ) ∈ V 1 × V 2 × V f such that

$$\mathfrak { o } _ { J } \text{ with boundary continuous } \left ( \mathfrak { o } \right ) \text{ can we written as follows: } \text{ in } p = \left \langle p _ { 1 }, p _ { 2 }, p _ { f } \right \rangle & \in v _ { 1 } \times v _ { 2 } \times v _ { f } \text{ such } \text{ in} \\ m _ { 1 } \left ( \frac { \partial \Theta _ { 1 } } { \partial t }, v _ { 1 } \right ) + a _ { 1 } ( p _ { 1 }, v _ { 1 } ) + q _ { 1 2 } ( p _ { 1 } - p _ { 2 }, v _ { 1 } ) + q _ { 1 f } ( p _ { 1 } - p _ { f }, v _ { 1 } ) & = l ( v _ { 1 } ), \\ m _ { 2 } \left ( \frac { \partial \Theta _ { 2 } } { \partial t }, v _ { 2 } \right ) + a _ { 2 } ( p _ { 2 }, v _ { 2 } ) - q _ { 2 1 } ( p _ { 1 } - p _ { 2 }, v _ { 2 } ) + q _ { 2 f } ( p _ { 2 } - p _ { f }, v _ { 2 } ) & = l ( v _ { 2 } ), \\ m _ { f } \left ( \frac { \partial \Theta _ { f } } { \partial t }, v _ { f } \right ) + a _ { f } ( p _ { f }, v _ { f } ) - q _ { 1 f } ( p _ { 1 } - p _ { f }, v _ { f } ) - q _ { 2 f } ( p _ { 2 } - p _ { f }, v _ { f } ) & = l ( v _ { f } ), \\ \dots \vee \sqrt { \dots } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{\ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text { } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text
\dots \vee \sqrt { \dots } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{\ } \text{ } \text{ } \text{\ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{\ } \text{ } \text{ } \text
\end{$$

for ∀ ( v , v 1 2 , v f ) ∈ ˆ V 1 × ˆ V 2 × ˆ V f and

$$V _ { 1 } = V _ { 2 } = \{ v \in H ^ { 1 } ( \Omega ) \, \colon v = g, \, x \in \Gamma _ { D } \}, \ \ V _ { f } = \{ v \in H ^ { 1 } ( \gamma ) \, \colon v = g, \, x \in \Gamma _ { D } \},$$

$$\hat { V } _ { 1 } = V _ { 2 } = \{ v \in H ^ { 1 } ( \Omega ) \, \colon v = 0, \, x \in \Gamma _ { D } \}, \quad \hat { V } _ { f } = \{ v \in H ^ { 1 } ( \gamma ) \, \colon v = 0, \, x \in \Gamma _ { D } \}.$$

For bilinear and linear forms, we have

$$m _ { \alpha } \left ( \frac { \partial \Theta _ { \alpha } } { \partial t }, v _ { \alpha } \right ) = \int _ { \Omega } \frac { \partial \Theta _ { \alpha } } { \partial t } \, v _ { \alpha } \, d x, \, \ m _ { f } \left ( \frac { \partial \Theta _ { f } } { \partial t }, v _ { f } \right ) = \int _ { \gamma } \frac { \partial \Theta _ { f } } { \partial t } \, v _ { f } \, d s,$$

$$& \text{For bilinear and linear forms, we have} \\ & \quad m _ { \alpha } \left ( \frac { \partial \Theta _ { \alpha } } { \partial t }, v _ { \alpha } \right ) = \int _ { \Omega } \frac { \partial \Theta _ { \alpha } } { \partial t } \, v _ { \alpha } \, d x, \quad m _ { f } \left ( \frac { \partial \Theta _ { f } } { \partial t }, v _ { f } \right ) = \int _ { \gamma } \frac { \partial \Theta _ { f } } { \partial t } \, v _ { f } \, d s, \\ & \quad a _ { \alpha } ( p _ { \alpha }, v _ { \alpha } ) = \int _ { \Omega } k _ { \alpha } ( x, p _ { \alpha } ) \nabla p _ { \alpha } \cdot \nabla v _ { \alpha } \, d x, \quad a _ { f } ( p _ { f }, v _ { f } ) = \int _ { \gamma } k _ { f } ( x, p _ { f } ) \nabla p _ { f } \cdot \nabla v _ { f } \, d s, \\ & \quad l ( v _ { \alpha } ) = \int _ { \Omega } f _ { \alpha } \, v _ { \alpha } \, d x + \int _ { \Omega } \frac { \partial k _ { \alpha } ( x, p _ { \alpha } ) } { \partial z } v _ { \alpha } \, d x, \quad l ( v _ { f } ) = \int _ { \gamma } f _ { f } \, v _ { f } \, d s + \int _ { \gamma } \frac { \partial k _ { f } ( x, p _ { f } ) } { \partial z } v _ { f } \, d s, \\ & \quad q _ { 1 2 } ( p _ { 1 } - p _ { 2 }, v _ { 1 } ) = q _ { 2 1 } ( p _ { 1 } - p _ { 2 }, v _ { 2 } ) = \int _ { \Omega } \sigma _ { 1 2 } ( p _ { 1 } - p _ { 2 } ) \, v _ { 2 } \, d x, \\ & \quad q _ { \alpha f } ( p _ { \alpha } - p _ { f }, v _ { \alpha } ) = q _ { \alpha f } ( p _ { \alpha } - p _ { f }, v _ { f } ) = \int _ { \gamma } \sigma _ { \alpha f } ( \mathcal { P } ^ { \gamma } ( p _ { \alpha } ) - p _ { f } ) \, v _ { f } \, d s \\ & \text{for $\alpha,\beta=1,2.$} \\ & \text{$\dots$} \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \end{cases}$$

for α, β = 1 2. ,

Similarly, we can generalize model for multicontinuum case: find p = ( p , ..., p 1 L ) ∈ V 1 × ... × V L such that glyph[negationslash]

$$m _ { \alpha } \left ( \frac { \partial \Theta _ { \alpha } } { \partial t }, v _ { \alpha } \right ) + a _ { \alpha } ( p _ { \alpha }, v _ { \alpha } ) + \sum _ { \beta \neq \alpha } q _ { \alpha \beta } ( p _ { \alpha } - p _ { \beta }, v _ { \alpha } ) = l ( v _ { \alpha } ), \quad \alpha = 1, \dots, L,$$

for ∀ ( v , ..., v 1 L ) ∈ ˆ V 1 × ... × ˆ V L .

To solve the coupled nonlinear system of equations (9), we apply implicit approximation by time with a Picard iteration scheme. For Θ n +1 ,m +1 α we use Taylor series with respect to p α

$$\lambda _ { \alpha } ^ { n + 1, m + 1 } & \approx \Theta _ { \alpha } ^ { n + 1, m } + \frac { d \Theta _ { \alpha } } { d p _ { \alpha } } \Big | ^ { n + 1, m } ( p _ { \alpha } ^ { n + 1, m + 1 } - p _ { \alpha } ^ { n + 1, m } ) \\ & = \Theta _ { \alpha } ^ { n + 1, m } + C _ { \alpha } ^ { n + 1, m } ( p _ { \alpha } ^ { n + 1, m + 1 } - p _ { \alpha } ^ { n + 1, m } ),$$

where C α = d Θ α /dp α , superscript n indicate time iteration and m is the nonlinear iteration number. Therefore, we obtain following variational formulation glyph[negationslash]

$$c _ { \alpha } ^ { n + 1, m } & \left ( \frac { p _ { \alpha } ^ { n + 1, m + 1 } - p _ { \alpha } ^ { n } } { \tau }, v _ { \alpha } \right ) + a _ { \alpha } ^ { n + 1, m } ( \nabla p _ { \alpha } ^ { n + 1, m + 1 } \cdot \nabla v _ { \alpha } ) \\ & + \sum _ { \beta \neq \alpha } q _ { \alpha \beta } ^ { n + 1, m } ( p _ { \alpha } ^ { n + 1, m + 1 } - p _ { \beta } ^ { n + 1, m + 1 }, v _ { \alpha } ) = \eta ^ { n + 1, m } ( v _ { \alpha } ) - m _ { \alpha } \left ( \frac { \Theta _ { \alpha } ^ { n + 1, m } - \Theta _ { \alpha } ^ { n } } { \tau }, v _ { \alpha } \right ),$$

where α = 1 , ..., L , τ is the given time step and

$$c _ { \alpha } ^ { n + 1, m } ( p _ { \alpha }, v _ { \alpha } ) = \int _ { \Omega _ { \alpha } } C ^ { m + 1, m } p _ { \alpha } \, v _ { \alpha } \, d x, \quad l ^ { n + 1, m } ( v _ { \alpha } ) = \int _ { \Omega _ { \alpha } } f _ { \alpha } \, v _ { \alpha } \, d x + \int _ { \Omega _ { \alpha } } \frac { \partial k _ { \alpha } ( x, p _ { \alpha } ^ { n + 1, m } ) } { \partial z } v _ { \alpha } \, d x,$$

a n +1 ,m α ( p α , v α ) = ∫ Ω α k α ( x, p n +1 ,m α ) ∇ p α · ∇ v α dx, q n +1 ,m αβ ( p α -p β , v α ) = ∫ Ω α σ n +1 ,m αβ ( p α -p β ) v α dx.

For triple continuum model (8) with the Galerkin finite element method, we write the solution as

$$p _ { 1 } = \sum _ { i = 1 } ^ { N _ { f } ^ { \Omega } } p _ { 1, i } \phi _ { i }, \ \ p _ { 2 } = \sum _ { i = 1 } ^ { N _ { f } ^ { \Omega } } p _ { 2, i } \phi _ { i }, \ \ p _ { f } = \sum _ { i = 1 } ^ { N _ { f } ^ { \gamma } } p _ { f, i } \psi _ { i },$$

where φ i are the standard d - dimensional linear element basis functions defined on T h , and ψ i are ( d -1) dimensional linear basis function on the fracture mesh, N Ω f and N γ f denotes the number of the nodes for T h and for fracture mesh.

We write approximation for (8), in the matrix form as follows for p = ( p , p 1 2 , p f ) T

$$C ^ { n + 1, m } \frac { p ^ { n + 1, m + 1 } - p ^ { n } } { \tau } + A ^ { n + 1, m } p ^ { n + 1, m + 1 } = F ^ { n + 1, m },$$

where

$$C ^ { n + 1, m } = \begin{pmatrix} C _ { 1 } ^ { n + 1, m } & 0 & 0 \\ 0 & C _ { 2 } ^ { n + 1, m } & 0 \\ 0 & 0 & C _ { f } ^ { n + 1, m } \end{pmatrix}, \quad F ^ { n + 1, m } = \begin{pmatrix} F _ { 1 } ^ { n + 1, m } \\ F _ { 2 } ^ { n + 1, m } \\ F _ { f } ^ { n + 1, m } \end{pmatrix},$$

$$A ^ { n + 1, m } = \begin{pmatrix} A _ { 1 } ^ { n + 1, m } + Q _ { 1 2 } ^ { n + 1, m } + Q _ { 1 f } ^ { n + 1, m } & - Q _ { 1 2 } ^ { n + 1, m } & & - Q _ { 1 f } ^ { n + 1, m } \\ - Q _ { 1 2 } ^ { n + 1, m } & A _ { 2 } ^ { n + 1, m } + Q _ { 1 2 } ^ { n + 1, m } + Q _ { 2 f } ^ { n + 1, m } & & - Q _ { 2 f } ^ { n + 1, m } \\ - Q _ { 1 f } ^ { n + 1, m } & - Q _ { 2 f } ^ { n + 1, m } & A _ { f } ^ { n + 1, m } + Q _ { 1 f } ^ { n + 1, m } + Q _ { 2 f } ^ { n + 1, m } \end{pmatrix}.$$

This fine-scale discretization yields matrices of the size N f × N f with N f = N Ω f + N Ω f + N γ f .

Remark 1. In this paper, for simplification of the matrix construction, we use a modified DFM approach and consider the case when σ 2 f = 0 . We assume that p 1 = p f and using superposition principle [35, 17, 6], we eliminate p f from equations (12) and obtain following coupled system of equations for p = ( p , p 1 2 ) T

$$C ^ { n + 1, m } \frac { p ^ { n + 1, m + 1 } - p ^ { n } } { \tau } + A ^ { n + 1, m } p ^ { n + 1, m + 1 } = F ^ { n + 1, m },$$

where

$$\text{where} \\ C ^ { n + 1, m } = \begin{pmatrix} C ^ { n + 1, m } _ { 1 } + C ^ { n + 1, m } _ { f } & 0 \\ 0 & C ^ { n + 1, m } _ { 2 } \end{pmatrix}, \quad F ^ { n + 1, m } = \begin{pmatrix} F ^ { n + 1, m } _ { 1 } + F ^ { n + 1, m } _ { f } \\ F ^ { n + 1, m } _ { 2 } \end{pmatrix}, \\ A ^ { n + 1, m } = \begin{pmatrix} A ^ { n + 1, m } _ { 1 } + A ^ { n + 1, m } _ { f } + Q ^ { n + 1, m } _ { 1 2 } & - Q ^ { n + 1, m } _ { 1 2 } \\ - Q ^ { n + 1, m } _ { 1 2 } & A ^ { n + 1, m } _ { 2 } + Q ^ { n + 1, m } _ { 1 2 } \end{pmatrix}. \\ \dots \colon \dots \dots \dots \dots \tau \dots \tau \dots \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \tau \dots \tau \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \ tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau \tau \dots \tau$$

with matrices of the size N f × N f with N f = N Ω f + N Ω f .

## 4 Coarse grid approximation using GMsFEM

Let T H is the coarse grid for computational domain Ω and ω i is the local domains, where i = 1 , ..., N Ω c and N Ω c is the number of coarse grid nodes. For construction of the coarse grid solver, we use a Generalized Multiscale Finite Element Method (GMsFEM) [30, 36]. For the construction of the multiscale space for coarse grid approximation, we solve spectral problems in each local domain ω i , i = 1 , ..., N Ω c to identify the most important characteristics of the problem. Let p C = ∑ i,k p ψ i k i k ( x ), where ψ i k are the multiscale basis functions that supported in local domain ω i , and the index k represents the numbering of the basis functions, k = 1 , .., M i .

Figure 1: Illustration of the GMsFEM algorithm

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-eee5b0c1-97fd-42e4-95f6-a8757fcc0175.jpg)

In GMsFEM, we have offline and online steps (see Figure 1):

- · Offline steps:
- -Coarse grid generation and local domains construction, ω i .
- -Construction of the local multiscale basis functions, ψ i k in local domain ω i ( k = 1 , .., M i ).
- -Generation of the projection matrix, R T = ( ψ , ..., ψ 1 1 1 M 1 , ..., ψ N Ω c 1 , ..., ψ N Ω c M N c Ω ).
- · Online steps: For each time step and nonlinear iteration:
- -Construct coarse grid system using projection matrix.

- -Solve coarse grid system and reconstruct fine-scale solution.

Here, for a given configuration of heterogeneous properties and fracture geometry, we precompute a multiscale basis functions on the offline step and generate projection matrix. On the online step, we project nonlinear matrices into the multiscale space and solve coarse grid systems. Note that, we focus method on the structured coarse grids in this work for illustration of the method, but in general, the method in presented form can be extended to unstructured coarse grids.

In (7), we consider the case, when the nonlinearity and heterogeneity of k α ( x, p α ) is separable,

$$k _ { \alpha } ( x, p _ { \alpha } ) = k _ { \alpha, r } ( p _ { \alpha } ) k _ { \alpha, s } ( x ).$$

Therefore, we can use a linear multiscale space and precompute multiscale basis functions.

For general multicontinuum case, let p = ( p , ..., p 1 L ) ∈ V 1 × , ..., × V L and glyph[negationslash]

$$\alpha _ { \alpha } ^ { s } ( p _ { \alpha }, v _ { \alpha } ) + \sum _ { \beta \neq \alpha } q _ { \alpha \beta } ^ { s } ( p _ { \alpha } - p _ { \beta }, v _ { \alpha } ) = 0, \quad \alpha = 1,.., L,$$

where

$$\ w i c { a } _ { \alpha } ^ { s } ( p _ { \alpha }, v _ { \alpha } ) = \int _ { \Omega _ { \alpha } } k _ { \alpha, s } ( x ) \nabla p _ { \alpha } \cdot \nabla v _ { \alpha } \, d x, \ \ q _ { \alpha \beta } ^ { s } ( p _ { \alpha } - p _ { \beta }, v _ { \alpha } ) = \int _ { \Omega _ { \alpha } } \sigma _ { \alpha \beta } ^ { s } ( p _ { \alpha } - p _ { \beta } ) \, v _ { \alpha } \, d x$$

and V α = H 1 (Ω α ).

Next, we describe construction of the coupled mutiscale basis functions for the multicontinua flow problem. We start with the construction of the snapshot space in the local domain ω i . The constructed local snapshots contain the information about local heterogeneities. After, the multiscale spaces are obtained from the snapshot spaces by a dimension reduction via local spectral problems [6, 17, 37, 35].

The snapshot space is constructed by the solution of the local problems in ω i with all possible boundary conditions. Let V α = { v ∈ H 1 ( ω i ) : v = δ l on ∂ω i } ( α = 1 , ..., L ), we find φ i,l = ( φ i,l 1 , ..., φ i,l L ) ∈ V 1 × ... × V L such that glyph[negationslash]

$$\alpha _ { \alpha } ^ { s } ( \phi _ { \alpha } ^ { i, l }, v _ { \alpha } ) + \sum _ { \beta \neq \alpha } q _ { \alpha \beta } ^ { s } ( \phi _ { \alpha } ^ { i, l } - \phi _ { \beta } ^ { i, l }, v _ { \alpha } ) = 0, \quad \alpha = 1, \dots, L,$$

for ∀ v α ∈ V 0 α = { v ∈ H 1 ( ω i ) : v = 0on ∂ω i } . Here δ j is the function, which takes the value one at node x l ∈ ∂ω i and zero elsewhere, l = 1 , ..., J i ( J i is the number of nodes on the boundary of ω i ). Therefore, we define

$$V ^ { s n a p, i } = \text{span} \{ \phi _ { 1 } ^ { i, l }, \dots, \phi _ { J _ { i } } ^ { i, l } \}, \quad \text{and} \quad R _ { s n a p, i } ^ { T } = ( \phi _ { 1 } ^ { i, l }, \dots, \phi _ { J _ { i } } ^ { i, l } ).$$

For multiscale basis functions Ψ i = (Ψ i 1 , ..., Ψ ) i L T in ω i , we solve the local spectral problem on the snapshot space for

$$A _ { s n a p } \Psi ^ { s n a p, i } = \lambda S _ { s n a p } \Psi ^ { s n a p, i },$$

where Ψ i α,k = R snap,i Ψ snap,i α,k . We choose the smallest M i eigenvalues and use them for the construction of multiscale basis functions for k = 1 2 , .., M i .

Here

$$A ^ { s n a p } = R _ { s n a p, i } A ^ { s } R _ { s n a p, i } ^ { T }, \ \ A ^ { s n a p } = R _ { s n a p, i } S ^ { s } R _ { s n a p, i } ^ { T },$$

$$S ^ { s } = \begin{pmatrix} S ^ { s } _ { 1 } & 0 & \dots & 0 \\ 0 & S ^ { s } _ { 2 } & \dots & 0 \\ \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & S ^ { s } _ { L } \end{pmatrix}, \quad A ^ { s } = \begin{pmatrix} A ^ { s } _ { 1 } + \sum Q ^ { s } _ { 1 \beta } & - Q ^ { s } _ { 1 2 } & \dots & - Q ^ { s } _ { 1 L } \\ - Q ^ { s } _ { 1 2 } & A ^ { s } _ { 2 } + \sum Q ^ { s } _ { 2 \beta } & \dots & - Q ^ { s } _ { 2 L } \\ \dots & \dots & \dots & \dots \\ - Q ^ { s } _ { 1 L } & - Q ^ { s } _ { 2 L } & \dots & A ^ { s } _ { L } + \sum Q ^ { s } _ { L \beta } \end{pmatrix}.$$

and A s α = { a s α,lk } , S s α = { s s α,lk } for α = 1 , ..., L with

$$a _ { \alpha, l k } ^ { s } = \int _ { \omega _ { i } } k _ { \alpha, s } ( x ) \nabla \varphi _ { \alpha, l } \cdot \nabla \varphi _ { \alpha, k } d x, \ \ s _ { \alpha, l k } ^ { s } = \int _ { \omega _ { i } } k _ { \alpha, s } ( x ) \, \varphi _ { \alpha, l } \, \varphi _ { \alpha, k } \, d x.$$

where ϕ α,l are the linear basis functions. Note that, the adaptive approach can be applied based on the calculated eigenvalues λ k of the spectral problem (15).

Remark 2 (Simplified multiscale basis functions) . For the simplified heterogeneities cases, where fracture networks are known and have simplified geometries, the simplified approach for multiscale basis functions construction can be implemented. Simplified multiscale basis functions are constructed without solution of the local spectral problems on the snapshot spaces and defined by a direct setting specific boundary conditions in local problem (see [6] for detailed explanations).

For construction of the continuous multiscale space, let χ i are the linear basis functions in the local domains ω i on the coarse grid. The multiscale space is defined as

$$V _ { c } = s p a n ( \psi _ { 1 } ^ { 1 }, \dots, \psi _ { M _ { 1 } } ^ { 1 }, \dots, \psi _ { 1 } ^ { N _ { c } ^ { \Omega } }, \dots, \psi _ { M _ { N _ { c } ^ { \Omega } } ^ { \Omega } } ^ { \Omega } ),$$

with following projection matrix

$$R ^ { T } = ( \psi _ { 1 } ^ { 1 }, \dots, \psi _ { M _ { 1 } } ^ { 1 }, \dots, \psi _ { 1 } ^ { N _ { c } ^ { \Omega } }, \dots, \psi _ { M _ { N _ { c } ^ { \Omega } } ^ { \Omega } } ^ { N _ { c } ^ { \Omega } } ).$$

where ψ i k = χ i Ψ i k and χ i is the standard linear partition of unity functions.

Therefore, we have following coarse grid approximation for the multicontinua problem

$$C _ { c } ^ { n + 1, m } \frac { p _ { c } ^ { n + 1, m + 1 } - p _ { c } ^ { n + 1, m } } { \tau } + A _ { c } ^ { n + 1, m } p _ { c } ^ { n + 1, m + 1 } = F _ { c } ^ { n + 1, m },$$

where C n +1 ,m c = RC n +1 ,m R T , A n +1 ,m c = RA n +1 ,m R T , F n +1 ,m c = RF n +1 ,m and the fine grid solution can be reconstructed, p ms = R p T c .

## 5 Numerical results

We present numerical results for several model problems:

Test 1. Two-dimensional problem with Ω = [0 1] , , 2 p = ( p m , p f ) with T max = 66 · 10 -4 , k m,s = 10 and k f,s = 10 9 (homogeneous porous matrix).

Test 2. Two-dimensional problem with Ω = [0 1] , , 2 p = ( p m , p f ) with T max = 11 · 10 -4 , k m,s = k m,s ( x ) and k f,s = 10 9 (heterogeneous porous matrix).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-d15df541-9b5a-48f6-a4c3-e10c3af15fee.jpg)

Figure 2: Computational grids and heterogeneous properties for Test 2 and or Test 3 (two-dimensional problem). Left: coarse grid (blue color), fine grid (green color) and fracture geometry (red color). Middle: Heterogeneous coefficient k 1 ,s ( x ) for Test 3 . Right: Heterogeneous coefficient k 2 ,s ( x ) for Test 3 ( k m,s ( x ) for Test 2 ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-1f9f1263-b591-45d7-a453-be1e5af3fc3c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-50305402-1bf3-4b18-8950-5a4a8e6bf7e7.jpg)

Test 3. Two-dimensional problem with Ω 1 = Ω 2 = [0 1] , 2 , p = ( p , p 1 2 , p f ) with T max = 10 -4 , k 1 ,s = k 1 ,s ( x ), k 2 ,s = k 2 ,s ( x ) and k f,s = 10 . 9

Test 4. Three-dimensional problem with Ω = [0 10] , , 3 p = ( p m , p f ) with T max = 31 · 10 -4 , k m,s = 10 and k f,s = 10 9 (homogeneous porous matrix).

Test 5. Three-dimensional problem with Ω = [0 10] , , 3 p = ( p m , p f ) with T max = 1 25 . · 10 -4 , k m,s = k m,s ( x ) and k f,s = 10 9 (heterogeneous porous matrix).

Test 6. Three-dimensional problem with Ω 1 = Ω 2 = [0 10] , , 3 p = ( p , p 1 2 , p f ) with T max = 3 · 10 -4 , k 1 ,s = k 1 ,s ( x ), k 2 ,s = k 2 ,s ( x ) and k f,s = 10 . 9

In Figures 2 for and 3, we show computational domains and heterogeneous coefficients for test problems are presented, by green and blue color, we depict the fine and coarse grids. The coarse grids are uniform with 121 nodes (100 cells) for two-dimensional test problems and 726 nodes (500 cells) for three-dimensional tests. We use DOF f (Degree of Freedom) to denote fine grid system size and DOF c to denote problem size of the multiscale system using GMsFEM.

In Haverkamp model, we set A α = 1 511 . · 10 , 6 B α = 3 96, Θ . α,s = 0 287, Θ . α,r = 0 075, . C α = 1 175 . · 10 6 and D α = 4 74 for . α = m,f for Tests 1,2,4,5 and α = 1 2 , , f for Tests 3, 6 . We simulate for T max = 66 10 · -4 with 200 time steps and set initial condition p 0 = -61 5. . As boundary conditions, we set p g = -20 7 on the . top boundary (Γ D ) and zero flux on other boundaries (Γ N ).

To compare the results, we use the relative L 2 error, relative energy error and relative error for multicontinuum system between multiscale and fine-scale solutions for continuum α

$$e _ { \alpha } ^ { L _ { 2 } } = \frac { | | p _ { \alpha } - p _ { m s, \alpha } | | _ { L _ { 2 } } } { | | p _ { \alpha } | | _ { L _ { 2 } } }, \ \ e _ { \alpha } ^ { a } = \frac { | | p _ { \alpha } - p _ { m s, \alpha } | | _ { a } } { | | p _ { \alpha } | | _ { a } }, \ \ e ^ { q } = \frac { | | p - p _ { m s } | | _ { q } } { | | p | | _ { q } },$$

Figure 3: Computational grids and heterogeneous properties for Test 4 and or Test 5 (three-dimensional problem). First row: coarse grid (blue color), fine grid (green color) and fracture geometry (red color). Second row: Heterogeneous coefficient k 1 ,s ( x ) (left) and k 2 ,s ( x ) for Test 5 ( k m,s ( x ) for Test 4 ) (right)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-9c352149-c564-43ec-9a07-ae6d20f7f8b5.jpg)

where glyph[negationslash]

$$| | v _ { \alpha } | | _ { L _ { 2 } } ^ { 2 } = ( v _ { \alpha }, v _ { \alpha } ), \ \ | | v _ { \alpha } | | _ { a } ^ { 2 } = a _ { \alpha } ( v _ { \alpha }, v _ { \alpha } ), \ \ | | v | | _ { q } ^ { 2 } = \sum _ { \alpha } \left ( a _ { \alpha } ( v _ { \alpha }, v _ { \alpha } ) + \sum _ { \beta \neq \alpha } q _ { \alpha } ( v _ { \alpha } - v _ { \beta }, v _ { \alpha } ) \right ),$$

and p ms and p are the multiscale and reference (fine-scale) solutions. We use GMSH software [38] to construct computational domains and grids. The implementation is based on the open-source library FEniCS [39].

## 5.1 Two-dimensional model problems

We consider the unsaturated fluid flow in two-dimensional computational domain Ω = [0 1] , 2 for Tests 1,2 and 3 (Figure 2). Unstructured fine grid contains approximately 31 409 triangle cells and 14 376 vertices. We use 10 × 10 structured coarse grid with rectangular cells to find multiscale solution with 100 quadratic cells and 121 vertices .

In Test 1 , we consider problem for p = ( p m , p f ) with homogeneous porous matrix and fracture properties. The fine-scale solution and multiscale solution using 8 basis function are presented on Figure 4 for four time steps. Relative L 2 and energy errors are presented for different number of multiscale basis functions in Table

Figure 4: Numerical results for Test 1 . Solution p m and p ms,m for different times, t 25 , t 50 , t 100 and t 200 (from left to right). First row: fine scale solution DOF f = 14376. Second row: multiscale solution using 8 basis functions DOF c = 968. Third row: multiscale solution using simplified basis functions DOF c = 295

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-93858b10-bfa8-4d46-80fd-9eebeacd22dc.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-6e069658-792f-4b78-a95a-bfbf8bcde95e.jpg)

1. We presented error comparison for fixed number of multiscale basis functions in each local domain ω i ( M i = M for i = 1 , ..., N Ω c ) and investigate the adaptive approach that based on the choosing of the number of basis functions that related to the smallest eigenvalues. This adaptive approach, automatically identify number of basis functions that sufficient for obtaining a good results. Furthermore, we show the results for simplified multiscale basis functions that can be applied for some simple cases when we a priory can define main local modes of the solution. We observe that the 1 and 2 basis functions are not sufficient for obtaining a good results, but 4 multiscale basis functions can be sufficient for given geometry. Presented results demonstrates the dynamic of the multiscale solution behavior with different number of basis functions. We observe, a good convergence of the solution, when we took sufficient number of multiscale basis functions and with increasing number of basis functions.

Table 1: Numerical results for Test 1 . Relative L 2 and energy errors (%) for different number of multiscale basis functions

| Number of multiscale   | DOF c   | t 25    | t 25    | t 50    | t 50    | t 100   | t 100   | t 200   | t 200   |
|------------------------|---------|---------|---------|---------|---------|---------|---------|---------|---------|
| basis functions        |         | e L 2 m | e a m   | e L 2 m | e a m   | e L 2 m | e a m   | e L 2 m | e a m   |
| 1                      | 121     | 87.47   | 1605.09 | 87.44   | 2243.09 | 87.67   | 2064.03 | 87.87   | 2705.61 |
| 2                      | 242     | 48.06   | 1836.22 | 43.85   | 2406.58 | 40.39   | 2140.27 | 33.01   | 1989.89 |
| 4                      | 484     | 1.81    | 24.41   | 1.09    | 12.41   | 0.99    | 9.16    | 0.84    | 8.61    |
| 8                      | 968     | 0.31    | 3.16    | 0.21    | 1.85    | 0.18    | 1.25    | 0.15    | 1.21    |
| 12                     | 1452    | 0.15    | 1.35    | 0.09    | 0.69    | 0.09    | 0.62    | 0.08    | 0.56    |
| 16                     | 1936    | 0.12    | 1.05    | 0.08    | 0.54    | 0.05    | 0.26    | 0.03    | 0.16    |
| Simplified             | 295     | 1.73    | 28.23   | 1.06    | 17.72   | 1.19    | 17.67   | 1.01    | 17.11   |
| Adaptive by λ          | 290     | 2.48    | 30.95   | 1.54    | 18.74   | 2.13    | 24.47   | 1.81    | 24.35   |

Figure 5: Numerical results for Test 2 . Solution p m and p ms,m for different times, t 25 , t 50 , t 100 and t 200 (from left to right). First row: fine scale solution DOF f = 14376. Second row: multiscale solution using 8 basis functions DOF c = 968

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-a503026b-c166-46c6-9ab1-0a76406e9460.jpg)

Next, we consider Test 2 for p = ( p m , p f ) with heterogeneous matrix property k m,s ( x ) (see Figure 2). The fine-scale solution and multiscale solution using 8 basis function are presented on Figure 5. In Table 2, we present a dynamic of the L 2 and energy errors for different number of multiscale basis functions for four

Table 2: Numerical results for Test 2 . Relative L 2 and energy errors (%) for different number of multiscale basis functions

| Number of multiscale   | DOF c   | t 25    | t 25    | t 50    | t 50    | t 100   | t 100   | t 200   | t 200   |
|------------------------|---------|---------|---------|---------|---------|---------|---------|---------|---------|
| basis functions        |         | e L 2 m | e a m   | e L 2 m | e a m   | e L 2 m | e a m   | e L 2 m | e a m   |
| 1                      | 121     | 77.37   | 1204.73 | 76.57   | 1418.63 | 74.58   | 1087.39 | 67.86   | 1726.95 |
| 2                      | 242     | 13.27   | 187.27  | 19.76   | 450.39  | 17.92   | 318.65  | 10.41   | 212.56  |
| 4                      | 484     | 1.96    | 46.45   | 1.05    | 32.52   | 1.44    | 38.32   | 0.57    | 23.98   |
| 8                      | 968     | 0.37    | 3.54    | 0.23    | 1.73    | 0.36    | 3.23    | 0.11    | 0.56    |
| 12                     | 1452    | 0.23    | 1.31    | 0.14    | 0.88    | 0.21    | 1.81    | 0.05    | 0.24    |
| 16                     | 1936    | 0.13    | 0.91    | 0.07    | 0.39    | 0.11    | 0.78    | 0.02    | 0.09    |

time layers. We observe similarly good results as for previous homogeneous porous matrix property. The GMsFEM demonstrate a good convergence with at least 4 multiscale basis functions.

Table 3: Numerical results for Test 3 . Relative errors (%) for different number of multiscale basis functions

| Number of multiscale   | DOF c   | p 1     | p 1     | p 2     | p 2   | p q     |
|------------------------|---------|---------|---------|---------|-------|---------|
| basis functions        |         | e L 2 1 | e a 1   | e L 2 2 | e a 2 | e       |
| 1                      | 242     | 47.51   | 1228.31 | 0.46    | 2.71  | 1301.74 |
| 2                      | 484     | 9.12    | 385.99  | 0.27    | 1.98  | 63.74   |
| 4                      | 968     | 0.71    | 6.47    | 0.11    | 0.62  | 0.39    |
| 8                      | 1936    | 0.09    | 0.57    | 0.02    | 0.04  | 0.02    |
| 12                     | 2904    | 0.03    | 0.11    | 0.01    | 0.01  | 0.01    |
| 16                     | 3872    | 0.13    | 2.85    | 0.01    | 0.02  | 0.04    |

Next, we consider Test 3 for p = ( p , p 1 2 , p f ) with heterogeneous properties k 1 ,s ( x ) and k 2 ,s ( x ) with σ 12 = σ 21 = k 2 ,s ( x ). Reference (fine-scale) and multiscale solutions are presented in Figure 6. In Table 3, we present a relative errors and observe good convergence of the presented method with coupled multiscale basis functions for multicontinua unsaturated flow problem. For example, when we take 4 multiscale basis functions we have less than one percent of L 2 errors for p 1 and p 2 with near 5 % of energy error.

For two-dimensional test problems ( Test 1, 2 and 3 ), we observe small errors, when we take at least 4 multiscale basis functions with huge reduction of the system size. For Test 1 and 2 , the fine grid system has size DOF f = 14376, but using coarse grid multiscale solver with 4 basis functions, we reduce size to DOF c = 484. For Test 3 , fine grid and coarse grid systems have size DOF f = 28752 and DOF C = 968 for 4 multiscale basis functions, respectively.

Figure 6: Numerical results for Test 3 . Solution p , p 1 2 and p ms, 1 , p ms, 2 for different times, t 25 , t 50 , t 100 and t 200 (from left to right). First row: fine scale solution. (First continuum). Second row: multiscale solution using 8 basis functions (First continuum). Third row: fine scale solution (Second continuum). Fourth row: multiscale solution using 8 basis functions (Second continuum). DOF f = 28752, DOF c = 968

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-1f0ea3b7-8b98-4659-8ca3-d2a3777ca620.jpg)

## 5.2 Three-dimensional model problems

In this section, we consider three-dimensional model problems and present a numerical results for Test 4,5 and 6 . We consider, similar test problems as for two - dimensional problems. We start with Test 4 with homogeneous porous matrix properties and fracture distribution that presented on Figure 3. Heterogeneous coefficient for Test 4 and 5 are presented in Figure 3. Fine grid has 125740 tetrahedral cells and 21609

vertices. In GMsFEM, we use 5 × 5 × 5 structured coarse grid with 125 cubic cells and 216 vertices.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-5f2929dc-bf6e-46f0-8041-dccd6629a4cf.jpg)

Figure 7: Numerical results for Test 4 . Solution p m and p ms,m for final time. Left: fine scale solution DOF f = 21609. Right: multiscale solution using 16 basis functions DOF c = 3456

Table 4: Numerical results for Test 4 . Relative errors (%) for different number of multiscale basis functions

| Number of multiscale basis functions   | DOF c   | t 25    | t 25   | t 50    | t 50   | t 100   | t 100   | t 200   | t 200   |
|----------------------------------------|---------|---------|--------|---------|--------|---------|---------|---------|---------|
| Number of multiscale basis functions   | DOF c   | e L 2 m | e a m  | e L 2 m | e a m  | e L 2 m | e a m   | e L 2 m | e a m   |
| 1                                      | 216     | 29.68   | 102.99 | 36.41   | 101.46 | 37.23   | 103.83  | 17.39   | 108.95  |
| 2                                      | 432     | 22.04   | 194.11 | 31.17   | 286.69 | 31.23   | 292.51  | 13.91   | 240.97  |
| 4                                      | 864     | 6.73    | 89.81  | 2.21    | 13.91  | 4.59    | 44.39   | 0.95    | 10.79   |
| 8                                      | 1728    | 5.07    | 61.19  | 2.21    | 13.33  | 2.09    | 18.59   | 0.74    | 7.05    |
| 12                                     | 2592    | 2.32    | 18.16  | 0.73    | 4.74   | 1.05    | 7.33    | 0.31    | 2.78    |
| 16                                     | 3456    | 1.76    | 14.35  | 0.68    | 4.26   | 0.73    | 4.19    | 0.21    | 1.77    |

Table 5: Numerical results for Test 5 . Relative errors (%) for different number of multiscale basis functions

| Number of multiscale basis functions   | DOF c   | t 25    | t 25   | t 50    | t 50   | t 100   | t 100   | t 200   | t 200   |
|----------------------------------------|---------|---------|--------|---------|--------|---------|---------|---------|---------|
| Number of multiscale basis functions   | DOF c   | e L 2 m | e a m  | e L 2 m | e a m  | e L 2 m | e a m   | e L 2 m | e a m   |
| 1                                      | 216     | 27.64   | 105.15 | 34.61   | 100.73 | 38.99   | 103.66  | 28.34   | 106.47  |
| 2                                      | 432     | 18.97   | 174.97 | 28.97   | 260.03 | 33.77   | 311.93  | 22.62   | 264.14  |
| 4                                      | 864     | 6.82    | 86.51  | 2.72    | 23.15  | 3.37    | 31.11   | 2.38    | 21.24   |
| 8                                      | 1728    | 7.61    | 116.12 | 3.17    | 28.67  | 2.16    | 18.87   | 1.23    | 10.15   |
| 12                                     | 2592    | 3.97    | 39.72  | 1.36    | 10.52  | 1.13    | 8.79    | 0.72    | 5.61    |
| 16                                     | 3456    | 4.12    | 43.41  | 1.46    | 10.41  | 0.92    | 5.99    | 0.55    | 3.83    |

The fine-scale solution and multiscale solution using 16 basis function are presented on Figures 7 and

Figure 8: Numerical results for Test 5 . Solution p m and p ms,m for final time. Left: fine scale solution DOF f = 21609. Right: multiscale solution using 16 basis functions DOF c = 3456

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-fed7bb91-33d2-451b-ab2c-af7d9d51192b.jpg)

8 for Test 4 and 5 , respectively. In Tables 4 and 5, we present a dynamic for the relative L 2 and energy errors for Test 4 and 5 . From the results we can conclude that this method has a good performance for three-dimensional problems with heterogeneous and homogeneous porous matrix properties. For example, when we take 16 multiscale basis functions, we obtain solution with 0 21% of relative . L 2 error and 1 77% . of energy error at final time for Test 4 . Fine grid system has size DOF f = 21609 and DOF c = 3456 for multiscale solver with 16 basis functions.

Table 6: Numerical results for Test 6 . Relative errors (%) for different number of multiscale basis functions

| Number of multiscale basis functions   | DOF c   | p 1     | p 1    | p 2     | p 2    | p q     |
|----------------------------------------|---------|---------|--------|---------|--------|---------|
|                                        |         | e L 2 1 | e a 1  | e L 2 2 | e a 2  | e       |
| 2                                      | 864     | 21.77   | 263.19 | 12.88   | 111.66 | 2849.94 |
| 4                                      | 1728    | 2.23    | 19.44  | 1.17    | 7.56   | 3.64    |
| 8                                      | 3456    | 1.18    | 9.89   | 0.68    | 3.74   | 0.29    |
| 12                                     | 5184    | 0.65    | 5.53   | 0.34    | 1.91   | 0.12    |
| 16                                     | 6912    | 0.49    | 3.73   | 0.24    | 1.15   | 0.08    |

Fine grid and multiscale solutions for Test 6 with p = ( p , p 1 2 , p f ) are presented in Figure 9 for final time. Relative errors are presented in Table 6 for final time. Coefficient σ 12 = σ 21 = 10 3 · k 2 ,s ( x ). Fine grid system has size DOF f = 43218 and DOF c = 6912 for multiscale solver with 16 basis functions. We observe a good convergence of the presented method with coupled multiscale basis functions for unsaturated multicontinua flow problems in three - dimensional formulations.

Figure 9: Numerical results for Test 6 . Solution p , p 1 2 and p ms, 1 , p ms, 2 for final time. Left: fine scale solution DOF f = 43218. Right: multiscale solution using 16 basis functions DOF c = 6912

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00247v1-figure-8e47a9f5-3dbf-404c-856b-65415aa631c3.jpg)

## 6 Conclusion

We presented a multiscale method for simulations of the multicontinua unsaturated flow problems in heterogeneous fractured porous media. The mathematical model is presented as general multicontinua model. The fine grid approximation are presented for coupled system of equations based on the finite element approximation and Discrete Fracture Model with unstructured grids that resolve fracture geometries explicitly. We presented the multiscale solver using GMsFEM with coupled multiscale basis function construction for multicontinua problems. Adaptive approach is investigated with simplified multiscale basis functions. To illustrate the idea of our approach, we considered unsaturated flow in fractured porous media and dual continua background model with discrete fractures networks. Numerical results are presented for two and three dimensional test problems. Presented multiscale solver accurately capture detailed interactions between multiple continua and provide a good results with huge reduction of the discrete system size. The comparison of the relative error for different number of basis functions and for adaptive approaches are presented for six test problems.

## 7 Acknowledgements

MV's and DS's works are supported by the mega-grant of the Russian Federation Government N14.Y26.31.0013 and RSF N17-71-20055. The research of Eric Chung is partially supported by the Hong Kong RGC General Research Fund (Project numbers 14304217 and 14302018) and CUHK Faculty of Science Direct Grant 2017-18.

## References

- [1] Michael A Celia, Efthimios T Bouloutas, and Rebecca L Zarba. A general mass-conservative numerical solution for the unsaturated flow equation. Water resources research , 26(7):1483-1496, 1990.
- [2] Michael A Celia and Philip Binning. A mass conservative numerical solution for two-phase flow in porous media with application to unsaturated flow. Water Resources Research , 28(10):2819-2828, 1992.
- [3] Roland Haverkamp, Michel Vauclin, Jaoudat Touma, PJ Wierenga, and Georges Vachaud. A comparison of numerical simulation models for one-dimensional infiltration 1. Soil Science Society of America Journal , 41(2):285-294, 1977.
- [4] Seong H Lee, MF Lough, and CL Jensen. Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water resources research , 37(3):443-455, 2001.
- [5] Liyong Li, Seong Hee Lee, et al. Efficient field-scale simulation for black oil in a naturally fractured reservoir via discrete fracture networks and homogenized media. In International oil &amp; gas conference and exhibition in China . Society of Petroleum Engineers, 2006.
- [6] Eric T Chung, Yalchin Efendiev, Tat Leung, and Maria Vasilyeva. Coupling of multiscale and multicontinuum approaches. GEM-International Journal on Geomathematics , 8(1):9-41, 2017.
- [7] Qiuqi Li, Yuhe Wang, and Maria Vasilyeva. Multiscale model reduction for fluid infiltration simulation through dual-continuum porous media with localized uncertainties. Journal of Computational and Applied Mathematics , 336:127-146, 2018.
- [8] Yalchin Efendiev and Thomas Y Hou. Multiscale finite element methods: theory and applications. 4, 2009.
- [9] Hadi Hajibeygi, Giuseppe Bonfigli, Marc Andre Hesse, and Patrick Jenny. Iterative multiscale finitevolume method. Journal of Computational Physics , 227(19):8604-8621, 2008.
- [10] E Weinan, Bjorn Engquist, Xiantao Li, Weiqing Ren, and Eric Vanden-Eijnden. Heterogeneous multiscale methods: a review. Commun. Comput. Phys , 2(3):367-450, 2007.
- [11] Eric Chung, Yalchin Efendiev, and Thomas Y Hou. Adaptive multiscale model reduction with generalized multiscale finite element methods. Journal of Computational Physics , 320:69-95, 2016.

| [12]   | Eric T Chung, Yalchin Efendiev, and Wing Tat Leung. Constraint energy minimizing generalized multiscale finite element method. Computer Methods in Applied Mechanics and Engineering , 339:298- 319, 2018.                              |
|--------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [13]   | Xinguang He and Li Ren. An adaptive multiscale finite element method for unsaturated flow problems in heterogeneous porous media. Journal of hydrology , 374(1-2):56-70, 2009.                                                          |
| [14]   | Xinguang He and Li Ren. A multiscale finite element linearization scheme for the unsaturated flow problems in heterogeneous porous media. Water resources research , 42(8), 2006.                                                       |
| [15]   | Victor Eralingga Ginting. Computational upscaled modeling of heterogeneous porous media flow uti- lizing finite volume method. 2005.                                                                                                    |
| [16]   | Zhiming Chen, Weibing Deng, and Huang Ye. Upscaling of a class of nonlinear parabolic equations for the flow transport in heterogeneous porous media. Communications in Mathematical Sciences , 3(4):493-515, 2005.                     |
| [17]   | IY Akkutlu, Yalchin Efendiev, and Maria Vasilyeva. Multiscale model reduction for shale gas transport in fractured media. Computational Geosciences , pages 1-21, 2015.                                                                 |
| [18]   | Matei Tene, ¸ Mohammed Saad Al Kobaisi, and Hadi Hajibeygi. Algebraic multiscale method for flow in heterogeneous porous media with embedded discrete fractures (f-ams). Journal of Computational Physics , 321:819-845, 2016.          |
| [19]   | Sebastian Bosma, Hadi Hajibeygi, Matei Tene, and Hamdi A Tchelepi. Multiscale finite volume method for discrete fracture modeling on unstructured grids (ms-dfm). Journal of Computational Physics , 2017.                              |
| [20]   | Matei Tene, ¸ Sebastian BM Bosma, Mohammed Saad Al Kobaisi, and Hadi Hajibeygi. Projection-based embedded discrete fracture model (pedfm). Advances in Water Resources , 105:205-216, 2017.                                             |
| [21]   | H. Hajibeygi, D. Kavounis, and P. Jenny. A hierarchical fracture model for the iterative multiscale finite volume method. Journal of Computational Physics , 230(24):8729-8743, 2011.                                                   |
| [22]   | I Yucel Akkutlu, Yalchin Efendiev, Maria Vasilyeva, and Yuhe Wang. Multiscale model reduction for shale gas transport in poroelastic fractured media. Journal of Computational Physics , 353:356-376, 2018.                             |
| [23]   | Eric T Chung, Yalchin Efendiev, Wing Tat Leung, Maria Vasilyeva, and Yating Wang. Non-local multi-continua upscaling for flows in heterogeneous fractured media. Journal of Computational Physics , 372:22-34, 2018.                    |
| [24]   | Maria Vasilyeva, Eric T Chung, Siu Wun Cheung, Yating Wang, and Georgy Prokopev. Nonlocal multi- continua upscaling for multicontinua flow problems in fractured porous media. Journal of Computational and Applied Mathematics , 2019. |

| [25]   | Maria Vasilyeva, Eric T Chung, Yalchin Efendiev, and Jihoon Kim. Constrained energy minimization based upscaling for coupled flow and mechanics. Journal of Computational Physics , 376:660-674, 2019.                                                |
|--------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [26]   | I Yucel Akkutlu, Yalchin Efendiev, Maria Vasilyeva, and Yuhe Wang. Multiscale model reduction for shale gas transport in poroelastic fractured media. Journal of Computational Physics , 353:356-376, 2018.                                           |
| [27]   | Eric T Chung, Wing Tat Leung, Maria Vasilyeva, and Yating Wang. Multiscale model reduction for transport and flow problems in perforated domains. Journal of Computational and Applied Mathematics , 330:519-535, 2018.                               |
| [28]   | Eric T Chung, Yalchin Efendiev, Guanglian Li, and Maria Vasilyeva. Generalized multiscale finite element methods for problems in perforated heterogeneous domains. Applicable Analysis , 95(10):2254- 2279, 2016.                                     |
| [29]   | Yalchin Efendiev, Thomas Y Hou, Victor Ginting, et al. Multiscale finite element methods for nonlinear problems and their applications. Communications in Mathematical Sciences , 2(4):553-589, 2004.                                                 |
| [30]   | Yalchin Efendiev, Juan Galvis, and Thomas Y Hou. Generalized multiscale finite element methods (gmsfem). Journal of Computational Physics , 251:116-135, 2013.                                                                                        |
| [31]   | Vincent Martin, J´rˆme e o Jaffr´, e and Jean E Roberts. Modeling fractures and barriers as interfaces for flow in porous media. SIAM Journal on Scientific Computing , 26(5):1667-1691, 2005.                                                        |
| [32]   | Carlo D'Angelo and Anna Scotti. A mixed finite element method for darcy flow in fractured porous media with non-matching grids. ESAIM: Mathematical Modelling and Numerical Analysis , 46(2):465- 489, 2012.                                          |
| [33]   | Luca Formaggia, Alessio Fumagalli, Anna Scotti, and Paolo Ruffo. A reduced model for darcy's problem in networks of fractures. ESAIM: Mathematical Modelling and Numerical Analysis , 48(4):1089-1116, 2014.                                          |
| [34]   | Carlo D'angelo and Alfio Quarteroni. On the coupling of 1d and 3d diffusion-reaction equations: applica- tion to tissue perfusion problems. Mathematical Models and Methods in Applied Sciences , 18(08):1481- 1504, 2008.                            |
| [35]   | Yalchin Efendiev, Seong Lee, Guanglian Li, Jun Yao, and Na Zhang. Hierarchical multiscale modeling for flows in fractured media using generalized multiscale finite element method. GEM-International Journal on Geomathematics , 6(2):141-162, 2015. |
| [36]   | Y. Efendiev and T. Hou. Multiscale Finite Element Methods: Theory and Applications. 4, 2009.                                                                                                                                                          |
| [37]   | I Yucel Akkutlu, Yalchin Efendiev, Maria Vasilyeva, and Yuhe Wang. Multiscale model reduction for shale gas transport in a coupled discrete fracture and dual-continuum porous media. Journal of Natural Gas Science and Engineering , 2017.          |

- [38] Software gmsh. (http://geuz.org/gmsh/).
- [39] Anders Logg, Kent-Andre Mardal, and Garth Wells. Automated solution of differential equations by the finite element method: The fenics book. 84, 2012.
Denis Spiridonov Maria Vasilyeva Eric T. Chung
Generalized Multiscale Finite Element method for multicontinua unsaturated flow problems in fractured porous media
In this paper, we present a multiscale method for simulations of the
multicontinua unsaturated flow problems in heterogeneous fractured porous
media. The mathematical model is described by the system of Richards equations
for each continuum that coupled by the specific transfer term. To illustrate
the idea of our approach, we consider a dual continua background model with
discrete fractures networks that generalized as a multicontinua model for
unsaturated fluid flow in the complex heterogeneous porous media. We present
fine grid approximation based on the finite element method and Discrete
Fracture Model approach. In this model, we construct an unstructured fine grid
that take into account a complex fracture geometries for two and three
dimensional formulations. Due to construction of the unstructured grid, the
fine grid approximation leads to the very large system of equations. For
reduction of the discrete system size, we develop a multiscale method for
coarse grid approximation of the coupled problem using Generalized Multiscale
Finite Element Method. In this method, we construct a coupled multiscale basis
functions that used to construct highly accurate coarse grid approximation. The
multiscale method allowed us to capture detailed interactions between multiple
continua. We investigate accuracy of the proposed method for the several test
problems in two and three dimensional formulations. We present a comparison of
the relative error for different number of basis functions and for adaptive
approach. Numerical results illustrate that the presented method provide
accurate solution of the unsaturated multicontinua problem on the coarse grid
with huge reduction of the discrete system size.