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

TL;DR

This paper introduces a multiscale finite element method for efficient and accurate simulation of complex multicontinua unsaturated flow in fractured porous media, significantly reducing computational complexity.

Contribution

The paper develops a generalized multiscale finite element approach tailored for multicontinua unsaturated flow in fractured media, enabling coarse grid simulations with high accuracy.

Findings

High accuracy coarse grid approximation achieved

Significant reduction in system size demonstrated

Method effective for 2D and 3D complex fracture networks

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 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.

Tables

Table 1. Table 1: Numerical results for Test 1 . Relative L 2 subscript 𝐿 2 L_{2} and energy errors ( % percent \% ) for different number of multiscale basis functions

Number of multiscale	$DO{F}_c$	$t_{25}$		$t_{50}$		$t_{100}$		$t_{200}$
basis functions		$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$
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 $\lambda$	290	2.48	30.95	1.54	18.74	2.13	24.47	1.81	24.35

Table 2. Table 2: Numerical results for Test 2 . Relative L 2 subscript 𝐿 2 L_{2} and energy errors ( % percent \% ) for different number of multiscale basis functions

Number of multiscale	$DO{F}_c$	$t_{25}$		$t_{50}$		$t_{100}$		$t_{200}$
basis functions		$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$
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

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

Number of multiscale	$DO{F}_c$	$p_1$		$p_2$		p
basis functions		$e_1^{L_2}$	$e_1^a$	$e_2^{L_2}$	$e_2^a$	$e^q$
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

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

Number of multiscale	$DO{F}_c$	$t_{25}$		$t_{50}$		$t_{100}$		$t_{200}$
basis functions		$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$
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. Table 5: Numerical results for Test 5 . Relative errors ( % percent \% ) for different number of multiscale basis functions

Number of multiscale	$DO{F}_c$	$t_{25}$		$t_{50}$		$t_{100}$		$t_{200}$
basis functions		$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$	$e_m^{L_2}$	$e_m^a$
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

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

Number of multiscale	$DO{F}_c$	$p_1$		$p_2$		p
basis functions		$e_1^{L_2}$	$e_1^a$	$e_2^{L_2}$	$e_2^a$	$e^q$
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

Equations

$$\begin{array}{rl}& \frac{\partial {\Theta }_m}{\partial t}+\nabla \cdot q_m+L_{mf}=f_m,x\in \Omega ,\\ & q_m=-k_m(x,p_m)\nabla (p_m+z),x\in \Omega ,\\ & \frac{\partial {\Theta }_f}{\partial t}+\nabla \cdot q_f-L_{fm}=f_f,x\in \gamma ,\\ & q_f=-k_f(x,p_f)\nabla (p_f+z),x\in \gamma ,\end{array}$$

$$\begin{array}{rl}& \frac{\partial {\Theta }_m}{\partial t}-\nabla \cdot (k_m(x,p_m)\nabla (p_m+z))+L_{mf}=f_m,x\in \Omega ,\\ & \frac{\partial {\Theta }_f}{\partial t}-\nabla \cdot (k_f(x,p_f)\nabla (p_f+z))-L_{fm}=f_f,x\in \gamma .\end{array}$$

$$\begin{array}{rl}& \frac{\partial {\Theta }_1}{\partial t}-\nabla \cdot (k_1(x,p_1)\nabla (p_1+z))+L_{12}+L_{1f}=f_1,x\in \Omega ,\\ & \frac{\partial {\Theta }_2}{\partial t}-\nabla \cdot (k_2(x,p_2)\nabla (p_2+z))-L_{21}+L_{2f}=f_2,x\in \Omega ,\\ & \frac{\partial {\Theta }_f}{\partial t}-\nabla \cdot (k_f(x,p_f)\nabla (p_f+z))-L_{f1}-L_{f2}=f_f,x\in \gamma .\end{array}$$

$$p_1=p_2=p_f=p^0,$$

$$\begin{array}{rl}& p_1=p_2=p_f=g,x\in {\Gamma }_D,\\ & q_1\cdot n=q_2\cdot n=q_f\cdot n=0,x\in {\Gamma }_N,\end{array}$$

$$\begin{array}{rl}& \frac{\partial {\Theta }_{\alpha }}{\partial t}-\nabla \cdot (k_{\alpha }(x,p_{\alpha })\nabla (p_{\alpha }+z))+\sum _{\beta \ne \alpha }{L}_{\alpha \beta }=f_{\alpha },\end{array}$$

$${\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 }}},$$

$$\begin{array}{r}{m}_1\left(\frac{\partial {\Theta }_1}{\partial t},v_1\right)+a_1(p_1,v_1)+q_{12}(p_1-p_2,v_1)+q_{1f}(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_{21}(p_1-p_2,v_2)+q_{2f}(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_{1f}(p_1-p_f,v_f)-q_{2f}(p_2-p_f,v_f)=l(v_f),\end{array}$$

$$V_1=V_2=\{v\in H^1(\Omega ):v=g,x\in {\Gamma }_D\},V_f=\{v\in H^1(\gamma ):v=g,x\in {\Gamma }_D\},$$

$${\hat{V}}_1=V_2=\{v\in H^1(\Omega ):v=0,x\in {\Gamma }_D\},{\hat{V}}_f=\{v\in H^1(\gamma ):v=0,x\in {\Gamma }_D\}.$$

$$m_{\alpha }\left(\frac{\partial {\Theta }_{\alpha }}{\partial t},v_{\alpha }\right)={\int }_{\Omega }\frac{\partial {\Theta }_{\alpha }}{\partial t}{v}_{\alpha }dx,m_f\left(\frac{\partial {\Theta }_f}{\partial t},v_f\right)={\int }_{\gamma }\frac{\partial {\Theta }_f}{\partial t}{v}_{f}ds,$$

$$a_{\alpha }(p_{\alpha },v_{\alpha })={\int }_{\Omega }{k}_{\alpha }(x,p_{\alpha })\nabla p_{\alpha }\cdot \nabla v_{\alpha }dx,a_f(p_f,v_f)={\int }_{\gamma }{k}_f(x,p_f)\nabla p_f\cdot \nabla v_{f}ds,$$

$$l(v_{\alpha })={\int }_{\Omega }{f}_{\alpha }{v}_{\alpha }dx+{\int }_{\Omega }\frac{\partial k_{\alpha }(x,p_{\alpha })}{\partial z}{v}_{\alpha }dx,l(v_f)={\int }_{\gamma }{f}_f{v}_{f}ds+{\int }_{\gamma }\frac{\partial k_f(x,p_f)}{\partial z}{v}_{f}ds,$$

$$q_{12}(p_1-p_2,v_1)=q_{21}(p_1-p_2,v_2)={\int }_{\Omega }{\sigma }_{12}(p_1-p_2)v_{2}dx,$$

$$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}ds$$

$$\begin{array}{r}{m}_{\alpha }\left(\frac{\partial {\Theta }_{\alpha }}{\partial t},v_{\alpha }\right)+a_{\alpha }(p_{\alpha },v_{\alpha })+\sum _{\beta \ne \alpha }{q}_{\alpha \beta }(p_{\alpha }-p_{\beta },v_{\alpha })=l(v_{\alpha }),\alpha =1,...,L,\end{array}$$

$$\begin{array}{rl}{\Theta }_{\alpha }^{n+1,m+1}& \approx {\Theta }_{\alpha }^{n+1,m}+{\frac{d{\Theta }_{\alpha }}{d{p}_{\alpha }}|}^{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}),\end{array}$$

$$\begin{array}{rl}{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 \ne \alpha }{q}_{\alpha \beta }^{n+1,m}(p_{\alpha }^{n+1,m+1}-p_{\beta }^{n+1,m+1},v_{\alpha })=l^{n+1,m}(v_{\alpha })-m_{\alpha }\left(\frac{{\Theta }_{\alpha }^{n+1,m}-{\Theta }_{\alpha }^n}{\tau },v_{\alpha }\right),\end{array}$$

$$c_{\alpha }^{n+1,m}(p_{\alpha },v_{\alpha })={\int }_{{\Omega }_{\alpha }}{C}^{n+1,m}{p}_{\alpha }{v}_{\alpha }dx,l^{n+1,m}(v_{\alpha })={\int }_{{\Omega }_{\alpha }}{f}_{\alpha }{v}_{\alpha }dx+{\int }_{{\Omega }_{\alpha }}\frac{\partial k_{\alpha }(x,p_{\alpha }^{n+1,m})}{\partial z}{v}_{\alpha }dx,$$

$$a_{\alpha }^{n+1,m}(p_{\alpha },v_{\alpha })={\int }_{{\Omega }_{\alpha }}{k}_{\alpha }(x,p_{\alpha }^{n+1,m})\nabla p_{\alpha }\cdot \nabla v_{\alpha }dx,q_{\alpha \beta }^{n+1,m}(p_{\alpha }-p_{\beta },v_{\alpha })={\int }_{{\Omega }_{\alpha }}{\sigma }_{\alpha \beta }^{n+1,m}(p_{\alpha }-p_{\beta })v_{\alpha }dx.$$

$$p_1=\sum _{i=1}^{N_f^{\Omega }}{p}_{1,i}{\varphi }_i,p_2=\sum _{i=1}^{N_f^{\Omega }}{p}_{2,i}{\varphi }_i,p_f=\sum _{i=1}^{N_f^{\gamma }}{p}_{f,i}{\psi }_i,$$

$$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},$$

$$C^{n+1,m}=\left(\begin{array}{ccc}{C}_1^{n+1,m}& 0& 0\\ 0& C_2^{n+1,m}& 0\\ 0& 0& C_f^{n+1,m}\end{array}\right),F^{n+1,m}=\left(\begin{array}{c}{F}_1^{n+1,m}\\ F_2^{n+1,m}\\ F_f^{n+1,m}\end{array}\right),$$

$$A^{n+1,m}=\left(\begin{array}{ccc}{A}_1^{n+1,m}+Q_{12}^{n+1,m}+Q_{1f}^{n+1,m}& -Q_{12}^{n+1,m}& -Q_{1f}^{n+1,m}\\ -Q_{12}^{n+1,m}& A_2^{n+1,m}+Q_{12}^{n+1,m}+Q_{2f}^{n+1,m}& -Q_{2f}^{n+1,m}\\ -Q_{1f}^{n+1,m}& -Q_{2f}^{n+1,m}& A_f^{n+1,m}+Q_{1f}^{n+1,m}+Q_{2f}^{n+1,m}\end{array}\right).$$

$$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},$$

$$C^{n+1,m}=\left(\begin{array}{cc}{C}_1^{n+1,m}+C_f^{n+1,m}& 0\\ 0& C_2^{n+1,m}\end{array}\right),F^{n+1,m}=\left(\begin{array}{c}{F}_1^{n+1,m}+F_f^{n+1,m}\\ F_2^{n+1,m}\end{array}\right),$$

$$A^{n+1,m}=\left(\begin{array}{cc}{A}_1^{n+1,m}+A_f^{n+1,m}+Q_{12}^{n+1,m}& -Q_{12}^{n+1,m}\\ -Q_{12}^{n+1,m}& A_2^{n+1,m}+Q_{12}^{n+1,m}\end{array}\right).$$

$$k_{\alpha }(x,p_{\alpha })=k_{\alpha ,r}(p_{\alpha })k_{\alpha ,s}(x).$$

$$\begin{array}{r}{a}_{\alpha }^s(p_{\alpha },v_{\alpha })+\sum _{\beta \ne \alpha }{q}_{\alpha \beta }^s(p_{\alpha }-p_{\beta },v_{\alpha })=0,\alpha =1,...,L,\end{array}$$

$$a_{\alpha }^s(p_{\alpha },v_{\alpha })={\int }_{{\Omega }_{\alpha }}{k}_{\alpha ,s}(x)\nabla p_{\alpha }\cdot \nabla v_{\alpha }dx,q_{\alpha \beta }^s(p_{\alpha }-p_{\beta },v_{\alpha })={\int }_{{\Omega }_{\alpha }}{\sigma }_{\alpha \beta }^s(p_{\alpha }-p_{\beta })v_{\alpha }dx$$

Full text

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 eigenvalues 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 $\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 following coupled system of equations

[TABLE]

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; $\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. Here $L_{mf}$ and $L_{fm}$ are the transfer term between matrix-fracture and fracture-matrix, $L_{mf}=\eta_{m}q_{mf}$, $L_{fm}=\eta_{f}q_{mf}$, $q_{mf}\approx\sigma_{mf}(x,p_{m},p_{f})(p_{m}-p_{f})$ and $\int_{\Omega}L_{mf}dx-\int_{\gamma}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}$:

[TABLE]

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}$:

[TABLE]

where subindices $1,2$ are related for first and second continua, and subindex $f$ is related to the lower-dimensional fracture model. System of equations is coupled by the mass exchange terms between continuum and $\int_{\Omega}L_{\alpha f}dx-\int_{\gamma}L_{f\alpha}ds=0$, $\int_{\Omega}L_{\alpha\beta}dx-\int_{\Omega}L_{\beta\alpha}dx=0$ for $\alpha,\beta=1,2$ and suppose following linear relation $L_{12}=L_{21}\approx\sigma_{12}(x,p_{1},p_{2})(p_{1}-p_{2})$.

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

[TABLE]

and boundary conditions

[TABLE]

where $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$, $\Gamma_{D}$ is the top boundary of the computational domain and $\Gamma_{N}$ denotes left, right and bottom boundaries of the $\Omega$. On the fracture boundaries that inside domain $\Omega$, we set a zero Neumann boundary conditions.

In general, we have following multicontinuum model:

[TABLE]

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

As constitutive relations, we use Haverkamp model [3]

[TABLE]

where $A_{\alpha}$ $B_{\alpha}$, $C_{\alpha}$, $D_{\alpha}$, $\Theta_{\alpha,r}$ and $\Theta_{\alpha,s}$ are the Haverkamp model coefficients, $k_{\alpha,s}(x)$ is the heterogeneous saturated hydraulic conductivity for $\alpha$ 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 $\mathcal{T}_{h}$ denote a finite element partition of the domain $\Omega\in\mathcal{R}^{d}$ and $\mathcal{T}_{h}=\bigcup_{i}K_{i}$, where $K_{i}$ is the triangular element for $d=2$ and tetrahedron for $d=3$. Let $\mathcal{E}_{h}$ the set of all the faces between the elements $\mathcal{T}_{h}$ and $\mathcal{E}_{\gamma}\subset\mathcal{E}_{h}$ be the subset of all faces that represent fractures. Furthermore, fracture facets $\mathcal{E}_{\gamma}$ 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_{1},p_{2},p_{f})\in V_{1}\times V_{2}\times V_{f}$ such that

[TABLE]

for $\forall(v_{1},v_{2},v_{f})\in\hat{V}_{1}\times\hat{V}_{2}\times\hat{V}_{f}$ and

[TABLE]

[TABLE]

For bilinear and linear forms, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

for $\alpha,\beta=1,2$.

Similarly, we can generalize model for multicontinuum case: find $p=(p_{1},...,p_{L})\in V_{1}\times...\times V_{L}$ such that

[TABLE]

for $\forall(v_{1},...,v_{L})\in\hat{V}_{1}\times...\times\hat{V}_{L}$.

To solve the coupled nonlinear system of equations (9), we apply implicit approximation by time with a Picard iteration scheme. For $\Theta_{\alpha}^{n+1,m+1}$ we use Taylor series with respect to $p_{\alpha}$

[TABLE]

where $C_{\alpha}=d\Theta_{\alpha}/dp_{\alpha}$, superscript $n$ indicate time iteration and $m$ is the nonlinear iteration number.

Therefore, we obtain following variational formulation

[TABLE]

where $\alpha=1,...,L$, $\tau$ is the given time step and

[TABLE]

[TABLE]

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

[TABLE]

where $\phi_{i}$ are the standard $d$ - dimensional linear element basis functions defined on $\mathcal{T}_{h}$, and $\psi_{i}$ are $(d-1)$ - dimensional linear basis function on the fracture mesh, $N_{f}^{\Omega}$ and $N_{f}^{\gamma}$ denotes the number of the nodes for $\mathcal{T}_{h}$ and for fracture mesh.

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

[TABLE]

where

[TABLE]

[TABLE]

This fine-scale discretization yields matrices of the size $N_{f}\times N_{f}$ with $N_{f}=N_{f}^{\Omega}+N_{f}^{\Omega}+N_{f}^{\gamma}$.

Remark 1**.**

In this paper, for simplification of the matrix construction, we use a modified DFM approach and consider the case when $\sigma_{2f}=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_{1},p_{2})^{T}$

[TABLE]

where

[TABLE]

[TABLE]

with matrices of the size $N_{f}\times N_{f}$ with $N_{f}=N_{f}^{\Omega}+N_{f}^{\Omega}$.

4 Coarse grid approximation using GMsFEM

Let $\mathcal{T}_{H}$ is the coarse grid for computational domain $\Omega$ and $\omega_{i}$ is the local domains, where $i=1,...,N_{c}^{\Omega}$ and $N_{c}^{\Omega}$ 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 $\omega_{i}$, $i=1,...,N_{c}^{\Omega}$ to identify the most important characteristics of the problem. Let $p_{C}=\sum_{i,k}p_{k}^{i}\psi_{k}^{i}(x)$, where $\psi_{k}^{i}$ are the multiscale basis functions that supported in local domain $\omega_{i}$, and the index $k$ represents the numbering of the basis functions, $k=1,..,M_{i}$.

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

• •

Offline steps:

• –

Coarse grid generation and local domains construction, $\omega_{i}$.

• –

Construction of the local multiscale basis functions, $\psi_{k}^{i}$ in local domain $\omega_{i}$ ($k=1,..,M_{i}$).

• –

Generation of the projection matrix, $R^{T}=(\psi_{1}^{1},...,\psi_{M_{1}}^{1},...,\psi_{1}^{N^{\Omega}_{c}},...,\psi_{M_{N^{\Omega}_{c}}}^{N^{\Omega}_{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_{\alpha}(x,p_{\alpha})$ is separable,

[TABLE]

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

For general multicontinuum case, let $p=(p_{1},...,p_{L})\in V_{1}\times,...,\times V_{L}$ and

[TABLE]

where

[TABLE]

and $V_{\alpha}=H^{1}(\Omega_{\alpha})$.

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 $\omega_{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 $\omega_{i}$ with all possible boundary conditions. Let $V_{\alpha}=\{v\in H^{1}(\omega_{i}):v=\delta_{l}\,\text{on}\,\partial\omega_{i}\}$ ($\alpha=1,...,L$), we find $\phi^{i,l}=(\phi^{i,l}_{1},...,\phi^{i,l}_{L})\in V_{1}\times...\times V_{L}$ such that

[TABLE]

for $\forall v_{\alpha}\in V^{0}_{\alpha}=\{v\in H^{1}(\omega_{i}):v=0\,\text{on}\,\partial\omega_{i}\}$. Here $\delta_{j}$ is the function, which takes the value one at node $x_{l}\in\partial\omega_{i}$ and zero elsewhere, $l=1,...,J^{i}$ ($J^{i}$ is the number of nodes on the boundary of $\omega_{i}$). Therefore, we define

[TABLE]

For multiscale basis functions $\Psi^{i}=(\Psi^{i}_{1},...,\Psi^{i}_{L})^{T}$ in $\omega_{i}$, we solve the local spectral problem on the snapshot space for

[TABLE]

where $\Psi^{i}_{\alpha,k}=R_{snap,i}\Psi^{snap,i}_{\alpha,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

[TABLE]

[TABLE]

and $A^{s}_{\alpha}=\{a^{s}_{\alpha,lk}\}$, $S^{s}_{\alpha}=\{s^{s}_{\alpha,lk}\}$ for $\alpha=1,...,L$ with

[TABLE]

where $\varphi_{\alpha,l}$ are the linear basis functions. Note that, the adaptive approach can be applied based on the calculated eigenvalues $\lambda_{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 $\chi_{i}$ are the linear basis functions in the local domains $\omega_{i}$ on the coarse grid. The multiscale space is defined as

[TABLE]

with following projection matrix

[TABLE]

where $\psi^{i}_{k}=\chi_{i}\Psi^{i}_{k}$ and $\chi_{i}$ is the standard linear partition of unity functions.

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

[TABLE]

where $C_{c}^{n+1,m}=RC^{n+1,m}R^{T}$, $A_{c}^{n+1,m}=RA^{n+1,m}R^{T}$, $F_{c}^{n+1,m}=RF^{n+1,m}$ and the fine grid solution can be reconstructed, $p_{ms}=R^{T}p_{c}$.

5 Numerical results

We present numerical results for several model problems:

• Test 1. Two-dimensional problem with $\Omega=[0,1]^{2}$, $p=(p_{m},p_{f})$ with $T_{max}=66\cdot 10^{-4}$, $k_{m,s}=10$ and $k_{f,s}=10^{9}$ (homogeneous porous matrix).
• Test 2. Two-dimensional problem with $\Omega=[0,1]^{2}$, $p=(p_{m},p_{f})$ with $T_{max}=11\cdot 10^{-4}$, $k_{m,s}=k_{m,s}(x)$ and $k_{f,s}=10^{9}$ (heterogeneous porous matrix).
• Test 3. Two-dimensional problem with $\Omega_{1}=\Omega_{2}=[0,1]^{2}$, $p=(p_{1},p_{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 $\Omega=[0,10]^{3}$, $p=(p_{m},p_{f})$ with $T_{max}=31\cdot 10^{-4}$, $k_{m,s}=10$ and $k_{f,s}=10^{9}$ (homogeneous porous matrix).
• Test 5. Three-dimensional problem with $\Omega=[0,10]^{3}$, $p=(p_{m},p_{f})$ with $T_{max}=1.25\cdot 10^{-4}$, $k_{m,s}=k_{m,s}(x)$ and $k_{f,s}=10^{9}$ (heterogeneous porous matrix).
• Test 6. Three-dimensional problem with $\Omega_{1}=\Omega_{2}=[0,10]^{3}$, $p=(p_{1},p_{2},p_{f})$ with $T_{max}=3\cdot 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_{\alpha}=1.511\cdot 10^{6}$, $B_{\alpha}=3.96$, $\Theta_{\alpha,s}=0.287$, $\Theta_{\alpha,r}=0.075$, $C_{\alpha}=1.175\cdot 10^{6}$ and $D_{\alpha}=4.74$ for $\alpha=m,f$ for Tests 1,2,4,5 and $\alpha=1,2,f$ for Tests 3, 6. We simulate for $T_{max}=66\cdot 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 ($\Gamma_{D}$) and zero flux on other boundaries ($\Gamma_{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 $\alpha$

[TABLE]

where

[TABLE]

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 $\Omega=[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\times 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 1. We presented error comparison for fixed number of multiscale basis functions in each local domain $\omega_{i}$ ($M_{i}=M$ for $i=1,...,N^{\Omega}_{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.

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 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.

Next, we consider Test 3 for $p=(p_{1},p_{2},p_{f})$ with heterogeneous properties $k_{1,s}(x)$ and $k_{2,s}(x)$ with $\sigma_{12}=\sigma_{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.

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\times 5\times 5$ structured coarse grid with $125$ cubic cells and $216$ vertices.

The fine-scale solution and multiscale solution using 16 basis function are presented on Figures 7 and 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.

Fine grid and multiscale solutions for Test 6 with $p=(p_{1},p_{2},p_{f})$ are presented in Figure 9 for final time. Relative errors are presented in Table 6 for final time. Coefficient $\sigma_{12}=\sigma_{21}=10^{3}\cdot 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.

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.