Multiscale Model Reduction of the Unsaturated Flow Problem in Heterogeneous Porous Media with Rough Surface Topography

Abstract

In this paper, we consider unsaturated filtration in heterogeneous porous media with rough surface topography. The surface topography plays an important role in determining the flow process and includes multiscale features. The mathematical model is based on the Richards’ equation with three different types of boundary conditions on the surface: Dirichlet, Neumann, and Robin boundary conditions. For coarse-grid discretization, the Generalized Multiscale Finite Element Method (GMsFEM) is used. Multiscale basis functions that incorporate small scale heterogeneities into the basis functions are constructed. To treat rough boundaries, we construct additional basis functions to take into account the influence of boundary conditions on rough surfaces. We present numerical results for two-dimensional and three-dimensional model problems. To verify the obtained results, we calculate relative errors between the multiscale and reference (fine-grid) solutions for different numbers of multiscale basis functions. We obtain a good agreement between fine-grid and coarse-grid solutions.

1. Introduction

Prediction of flows in unsaturated media is an important problem in many areas of science and engineering. In this paper, we consider the problem in heterogeneous porous media with rough surface topography. These problems occur in many applications related to the vadose (or unsaturated) zone. Although unsaturated flow conditions occur below the surface, the surface topography plays an important role in determining the flow process. Typical surface topography includes micro-scale features, such as roughness, as well as macro-scale features, such as hills, slopes, and valleys. It has been shown by previous studies on hydrological processes (e.g., [1,2,3]) that the surface topography exerts varying control over soil hydrological processes at different spatial scales. It has been demonstrated [2,4] that the topography plays a key role in redistribution of the surface water content into runoff and infiltration, especially at the hillslope (kilometer) scales. However, thus far, studies on mathematical modeling of the infiltration process have ignored the existence of surface topography and assumed a flat top boundary for the domain. While this is easier to solve mathematically, the physical reality is not matched by this assumption. In this study, we aim to address this mismatch by considering an uneven top boundary. For unsaturated infiltration, we formulate a mathematical model that is based on the Richards’ equation [5,6,7].

To describe the flow with rough boundaries, we consider three different types of boundary conditions: Dirichlet boundary conditions (represent head values on the boundary), Neumann boundary conditions (represent the flux on the boundary), and Robin boundary conditions (represent mixed boundary conditions). In addition, we consider our problem posed in highly heterogeneous media, where the conductivity varies in space [8,9,10]. Therefore, we are faced with the problem of a large number of unknowns of the discrete system. This leads to large computational resources. For this reason, we design a model reduction technique based on homogenization and multiscale approaches. Homogenization methods give macroscopic laws and parameters that are based on local computations. These approaches are often based on a priori assumptions [11,12]. However, the multiscale methods have a two-way information exchange between micro- and macro-scales.

Many multiscale methods can be used for such types of problems; for example, the multiscale finite element method (MsFEM) [9], heterogeneous multiscale methods (HMM) [13], multiscale finite volume method (MsFVM) [14], generalized multiscale finite element method (GMsFEM) [15,16], constraint energy minimizing generalized multiscale finite element method (CEM-GMsFEM) [17], and nonlocal multi-continuum method (NLMC) [18,19,20]. Multiscale methods can be applied to unsaturated filtration problems. For example, in [21,22,23], the authors present a multiscale methods for unsaturated flow in heterogeneous media. An upscaling method for saturated flow is described in [24].

In this paper, we solve an unsaturated filtration model with rough boundaries using the Generalized Multiscale Finite Element Method (GMsFEM) [9,25,26,27]. In GMsFEM, there are online and offline stages. In the offline stage, we construct an offline space computing multiscale basis functions by solving a local spectral problem on the snapshot space in each local domain. When we construct a snapshot space, we will not take into account the rough boundary of the computational domain. The main concept of snapshot space construction is that the snapshot vectors represent essential solution properties and give a good approximation space. Snapshot space helps to better take into account heterogeneities with high contrast, as well as complex heterogeneities, such as channels and fractures. To treat rough boundaries, we calculate additional basis functions to take into account the influence of top boundary conditions. Multiscale basis functions can describe small heterogeneities on the micro-scale and provide a good approximation on the macro-scale. This work is a continuation of the following works [23,28,29], where we used a similar technique of model reduction. In [28], we consider a multi-continuum filtration problem with GMsFEM; in [23], we describe a GMsFEM for an unsaturated filtration problem in heterogeneous media; in [29], we present a GMsFEM for a multi-continuum unsaturated filtration problem in fractured media. The GMsFEM framework that is presented in this paper is based on the listed works. We are expanding our approach with additional basis functions. In [30], we used a similar approach with additional basis functions for Robin boundary conditions.

Non-homogeneous boundary conditions occur in many applications. For example, pore-scale modeling and simulations of reactive flows have many applications in many branches of science, such as biology, physics, chemists, geomechanics, and geology [31,32,33]. We have already considered the GMsFEM for unsaturated flow filtration with inhomogeneous Dirichlet boundary conditions in the following papers [23,29]. In problems with rough boundaries, it is better to construct additional basis functions to take into account the influence of inhomogeneous boundary conditions. In this paper, we are going to investigate this approach for three types of boundary conditions.

We will present numerical results to illustrate the performance of our method. In our numerical results, we consider both two-dimensional and three-dimensional examples. In all cases, the surface boundaries are taken to be heterogeneous with smaller variations compared to the coarse-grid size. All numerical results show a good accuracy.

The paper is organized as follows. In Section 2, we present a mathematical model for unsaturated filtration in heterogeneous porous media with rough boundary conditions on the surface. In Section 3, we consider an approximation on the fine grid. In Section 4, we present the GMsFEM algorithm and construction of additional basis functions for rough boundaries. In Section 5, we present numerical results for two-dimensional and three-dimensional cases. We conclude the paper with a summary and future directions.

2. Mathematical Model

For the modeling of unsaturated flow in porous media, we use a mathematical model that is described by the nonlinear Richard’s equation in the domain $\mathsf{\Omega }$,

$$\frac{\partial \mathsf{\Theta }\left(p\right)}{\partial t}-\nabla ·\left(k(x,p)\nabla p\right)-\frac{\partial k(x,p)}{\partial z}=0,x\in \mathsf{\Omega },t>0,$$

(1)

with the initial condition

$$p=p_0\left(x\right),x\in \mathsf{\Omega },t=0.$$

(2)

In the above equation, $\mathsf{\Theta }$ is the water content representing the volume fraction of the porous medium filled with fluid, k is the unsaturated hydraulic conductivity tensor, and z represents the influence of the gravity on the flow processes.

The nonlinear coefficient $k(x,p)$ and the water content $\mathsf{\Theta }$ are related by the following constitutive relations (Havercamp model):

$$\mathsf{\Theta }\left(p\right)=\frac{\mathcal{A}({\mathsf{\Theta }}_s-{\mathsf{\Theta }}_r)}{\mathcal{A}+{\left|p\right|}^{\mathcal{B}}}+{\mathsf{\Theta }}_r,k(x,p)=k_s\left(x\right)\frac{\mathcal{S}}{\mathcal{S}+{\left|p\right|}^{\mathcal{B}}}+{\mathsf{\Theta }}_r,$$

(3)

where $k_s\left(x\right)$ is a heterogeneous coefficient modeling the saturated hydraulic conductivity and $\mathcal{A}$, $\mathcal{S}$, $\mathcal{B}$, ${\mathsf{\Theta }}_r$, ${\mathsf{\Theta }}_s$ are the Haverkamp model coefficients.

We consider domains with rough boundaries. To be specific, we consider a domain $\mathsf{\Omega }$, illustrated in Figure 1. We consider three types of boundary conditions on the top boundary ${\mathsf{\Gamma }}_D$:

• Dirichlet boundary condition:

  $$p=p_1,x\in {\mathsf{\Gamma }}_D,$$
• Neumann boundary condition:

  $$-k(x,p)\frac{\partial p}{\partial n}=g,x\in {\mathsf{\Gamma }}_D,$$
• Robin boundary condition:

  $$-k(x,p)\frac{\partial p}{\partial n}=\alpha (p-p_2),x\in {\mathsf{\Gamma }}_D.$$

On the other part of the boundary, denoted by ${\mathsf{\Gamma }}_N$, we impose the zero flux boundary condition

$$-k(x,p)\frac{\partial p}{\partial n}=0,x\in {\mathsf{\Gamma }}_N,$$

where ${\mathsf{\Gamma }}_N=\partial \mathsf{\Omega }/{\mathsf{\Gamma }}_D$.

3. Fine Grid Approximation

In this section, we present the fine grid discretization of (1). For the approximation of the time derivative, we use the backward Euler method coupled with simplified approximation from the previous time step:

$$C^n\frac{p^{n+1}-p^n}{\tau }-\nabla ·\left(k^n\nabla p^{n+1}\right)-\frac{\partial k^n}{\partial z}=0,n=0,1,\dots ,$$

(4)

where $\partial \mathsf{\Theta }\left(p\right)/\partial t=C\left(p\right)\partial p/\partial t$, $C=d\mathsf{\Theta }/dp$, $C^n=C\left(p^n\right)$, and $k^n=k(x,p^n)$. In (4), the number $\tau >0$ denotes the time step size, and the superscript n denotes the value at the time $t=n\tau$. For instance, $p^n$ denotes the value of the solution p at the time $t=n\tau$.

For the approximation with respect to the spatial variables, we use the standard conforming finite element method, which is based on the following variational formulations for each type of boundary condition. Given $p^n$, we find $p^{n+1}$ using the following formulations:

• Dirichlet boundary condition. Find $p^{n+1}\in V=\{v\in H^1\left(\mathsf{\Omega }\right):v\left(x\right)=p_1,x\in {\mathsf{\Gamma }}_D\}$ such that

  $${\int }_{\mathsf{\Omega }}{C}^n\frac{p^{n+1}-p^n}{\tau }vdx+{\int }_{\mathsf{\Omega }}\left(k^n\nabla p^{n+1},\nabla v\right)dx-{\int }_{\mathsf{\Omega }}\frac{\partial k^n}{\partial z}vdx=0,\forall v\in V_0,$$

  (5)

  where $V_0=\{v\in H^1\left(\mathsf{\Omega }\right):v\left(x\right)=0,x\in {\mathsf{\Gamma }}_D\}$.
• Neumann boundary condition. Find $p^{n+1}\in V=H^1\left(\mathsf{\Omega }\right)$ such that

  $${\int }_{\mathsf{\Omega }}{C}^n\frac{p^{n+1}-p^n}{\tau }vdx+{\int }_{\mathsf{\Omega }}\left(k^n\nabla p^{n+1},\nabla v\right)dx-{\int }_{\mathsf{\Omega }}\frac{\partial k^n}{\partial z}vdx+{\int }_{{\mathsf{\Gamma }}_D}gvds=0,\forall v\in V.$$

  (6)
• Robin boundary condition. Find $p^{n+1}\in V=H^1\left(\mathsf{\Omega }\right)$ such that

  $${\int }_{\mathsf{\Omega }}{C}^n\frac{p^{n+1}-p^n}{\tau }vdx+{\int }_{\mathsf{\Omega }}\left(k^n\nabla p^{n+1},\nabla v\right)dx-{\int }_{\mathsf{\Omega }}\frac{\partial k^n}{\partial z}vdx+{\int }_{{\mathsf{\Gamma }}_D}\alpha (p-p_2)vds=0,\forall v\in V.$$

  (7)

Let ${\mathcal{T}}_h$ be the fine grid with mesh size h for the domain $\mathsf{\Omega }$. We write the approximate solution as follows

$$p_h=\sum _{i=1}^{N_f}{p}_{h,i}{\psi }_i,$$

where $N_f$ is the number of interior fine grid vertices for (5) (or number of fine grid vertices for (6) and (7)) and $\left\{{\psi }_i\right\}$ are the conforming piecewise linear basis functions. Using $p_h={(p_{h,1},p_{h,2},\dots ,p_{h,N_h})}^T$ again to denote the vector of the required unknowns, we have the following matrix form for the fully discretized system:

$$S^n\frac{p_h^{n+1}-p_h^n}{\tau }+A^n{p}_h^{n+1}=F^n,$$

(8)

where

$$S^n=\{s_{ij}^n\},s_{ij}^n={\int }_{\mathsf{\Omega }}{C}^n{\psi }_i{\psi }_{j}dx.$$

The following summarize the definitions of ${}^n$ and $F^n$ for each of the considered boundary conditions:

• Dirichlet boundary condition:

  $$A^n=\{a_{ij}^n\},a_{ij}^n={\int }_{\mathsf{\Omega }}\left(k^n\nabla {\psi }_i,\nabla {\psi }_j\right)dx,$$

  $$F^n=\{f_j^n\},f_j^n={\int }_{\mathsf{\Omega }}\frac{\partial k^n}{\partial z}{\psi }_{j}dx.$$
• Neumann boundary condition:

  $$A^n=\{a_{ij}^n\},a_{ij}^n={\int }_{\mathsf{\Omega }}\left(k^n\nabla {\psi }_i,\nabla {\psi }_j\right)dx,$$

  $$F^n=\{f_j^n\},f_j^n={\int }_{\mathsf{\Omega }}\frac{\partial k^n}{\partial z}{\psi }_{j}dx-{\int }_{{\mathsf{\Gamma }}_D}g{\psi }_{j}ds.$$
• Robin boundary condition:

  $$A^n=\{a_{ij}^n\},a_{ij}^n={\int }_{\mathsf{\Omega }}\left(k^n\nabla {\psi }_i,\nabla {\psi }_j\right)dx+{\int }_{{\mathsf{\Gamma }}_D}\alpha {\psi }_i{\psi }_{j}ds,$$

  $$F^n=\{f_j^n\},f_j^n={\int }_{\mathsf{\Omega }}\frac{\partial k^n}{\partial z}{\psi }_{j}dx+{\int }_{{\mathsf{\Gamma }}_D}\alpha p_2{\psi }_{j}ds.$$

4. Coarse Grid Approximation Using GMsFEM

To reduce the size of the system, we construct a coarse grid approximation using the Generalized Multiscale Finite Element Method (GMsFEM) [9,27]. Let ${\mathcal{T}}_H={\cup }_i{K}_i$ be the coarse grid for computational domain $\mathsf{\Omega }$ with size H and $K_i$ be the i-th coarse grid cell. We define a local domain ${\omega }_i$ as the union of coarse grid cells around one coarse grid vertex, $i=1,\dots ,N_c$, and $N_c$ is the number of coarse grid nodes. For the construction of the multiscale space for coarse grid approximation, we solve spectral problems in each local domain ${\omega }_i$, $i=1,\dots ,N_c$ to identify the most important characteristics of the problem.

The GMsFEM consists of the offline and the online stages. The construction of snapshot space in offline stage solves local problems for different choices of input parameters. This space is used to construct the offline space via a spectral decomposition of the snapshot space. In many applications, the snapshot space can avoid expensive offline space construction, for example, in problems with a finely perforated medium. The offline space is constructed by spectral decomposition on the snapshot space, which is based on the eigenvalue problem. The spectral decomposition allows the selection of high-energy elements from the offline space. It chooses eigenvectors corresponding to the smallest eigenvalues. We obtain a set of multiscale basis functions ${\psi }_k^i$ that contain dominant information about micro-scale heterogeneities in the local domains. In the online stage, we compute the coarse space for an input parameter. The online space is computed via a spectral decomposition by using the eigenvectors corresponding to the smallest eigenvalues. The online coarse space is further used to solve the original global problem. An illustration of the GMsFEM algorithm is presented in Figure 2. We note that our temporal discretization is unconditionally stable. Some theoretical analysis of multiscale methods is presented in [27,34].

For the construction of multiscale basis functions, we first generate a snapshot space $V_{\mathrm{snap}}^{{\omega }_i}$ for each local domain ${\omega }_i$. The snapshot space is constructed by solving the following local problems:

$$\begin{array}{cc}\hfill -\nabla ·(k_s\left(x\right)\nabla {\varphi }_j^i)& =0,x\in {\omega }_i,\hfill \\ \hfill {\varphi }_j^i& ={\delta }_j,\partial {\omega }_i/{\mathsf{\Gamma }}_D,\hfill \\ \hfill {\varphi }_j^i& =0,{\mathsf{\Gamma }}_D,\hfill \end{array}$$

(9)

where ${\delta }_j$ is the discrete delta function that takes the value 1 at the j-th fine grid node $x=x_j$ and zero elsewhere ($j=1,\dots ,J^i$, $J^i$ is number of fine grid nodes on boundary $\partial {\omega }_i/{\mathsf{\Gamma }}_D$). Then, we can define the snapshot space and the projection matrix on the snapshot space in local domain ${\omega }_i$ as follows:

$$V^{snap,i}=\mathrm{span}\{{\varphi }_1^i,\dots ,{\varphi }_{J_i}^i\},\mathrm{and}{R}_{snap,i}={({\varphi }_1^i,\dots ,{\varphi }_{J_i}^i)}^T.$$

Note that we do not calculate snapshots for boundary ${\mathsf{\Gamma }}_D$, and we will construct separate basis functions for that.

To obtain the multiscale basis functions, we solve a local spectral problem on the snapshot space in the local domain ${\omega }_i$:

$$\tilde{A}{\tilde{\mathsf{\Psi }}}_{snap,j}^i=\lambda \tilde{S}{\tilde{\mathsf{\Psi }}}_{snap,j}^i,$$

(10)

with

$$\tilde{A}=R_{snap,i}{A}^{{\omega }_i}{R}_{snap,i}^T,\tilde{S}=R_{snap,i}{S}^{{\omega }_i}{R}_{snap,i}^T,$$

where

$$A^{{\omega }_i}=\left\{a_{ln}\right\},a_{ln}={\int }_{{\omega }_i}(k_s\left(x\right)\nabla {\psi }_l,\nabla {\psi }_n)dx,S^{{\omega }_i}=\left\{s_{ln}\right\},s_{ln}={\int }_{{\omega }_i}{k}_s\left(x\right){\psi }_l{\psi }_{n}dx.$$

(11)

We compute the solution of the spectral problem by ${\mathsf{\Psi }}_j^i=R_{snap,i}{\tilde{\mathsf{\Psi }}}_{snap,j}^i$. We note that for computation of the multiscale basis function, we use only the linear part $k_s\left(x\right)$ of coefficient $k(x,p)$ from (3).

Next, we select the smallest $M_i$ eigenvalues and use them for multiscale basis function construction (${\mathsf{\Psi }}_j^i$, $j=1,\dots ,M_i$). We obtain the required multiscale basis functions after multiplication by the linear partition of unity function

$${\psi }_j^i={\chi }_i{\mathsf{\Psi }}_j^i,$$

where ${\chi }_i$ is the standard coarse grid nodal basis function for the coarse node i in local domain ${\omega }_i$.

To handle rough boundary ${\mathsf{\Gamma }}_D$ with non-homogeneous boundary conditions, we calculate an additional basis function in local domains ${\omega }_i:\partial {\omega }_i\cap {\mathsf{\Gamma }}_D\ne \varnothing$. We solve the following local problem:

• Dirichlet boundary condition

  $$\begin{array}{cc}\hfill -\nabla ·(k_s\left(x\right)\nabla {\mathsf{\Phi }}^i)& =0,x\in {\omega }_i,\hfill \\ \hfill {\varphi }^i& =1,x\in {\mathsf{\Gamma }}_D,\hfill \\ \hfill {\varphi }^i& =0,x\in \partial {\omega }_i/{\mathsf{\Gamma }}_D;\hfill \end{array}$$

  (12)
• Neumann boundary condition

  $$\begin{array}{cc}\hfill -\nabla ·(k_s\left(x\right)\nabla {\mathsf{\Phi }}^i)& =0,x\in {\omega }_i,\hfill \\ \hfill -k_s\left(x\right)\frac{\partial {\mathsf{\Phi }}^i}{\partial n}& =1,x\in {\mathsf{\Gamma }}_D,\hfill \\ \hfill {\mathsf{\Phi }}^i& =0,x\in \partial {\omega }_i/{\mathsf{\Gamma }}_D;\hfill \end{array}$$

  (13)
• Robin boundary condition

  $$\begin{array}{cc}\hfill -\nabla ·(k_s\left(x\right)\nabla {\mathsf{\Phi }}^i)& =0,x\in {\omega }_i,\hfill \\ \hfill -k_s\left(x\right)\frac{\partial {\mathsf{\Phi }}^i}{\partial n}& =\alpha ({\mathsf{\Phi }}^i-1),x\in {\mathsf{\Gamma }}_D,\hfill \\ \hfill {\mathsf{\Phi }}^i& =0,x\in \partial {\omega }_i/{\mathsf{\Gamma }}_D.\hfill \end{array}$$

  (14)

To compute additional basis functions, we multiply the solution of the local problem to the linear partition of unity function ${\chi }_i\in {\omega }_i$

$${\varphi }^i={\chi }_i{\mathsf{\Phi }}^i.$$

The solutions of problems (12)–(14) are presented in Figure 3.

Finally, we construct the multiscale space

$$V_{ms}=\mathrm{span}({\psi }_1^1,\dots ,{\psi }_{M_1}^1,\dots ,{\psi }_1^{N_c},\dots ,{\psi }_{M_{N_c}}^{N_c},{\varphi }^1,\dots ,{\varphi }^L)$$

and projection matrix

$$R^T=({\psi }_1^1,\dots ,{\psi }_{M_1}^1,\dots ,{\psi }_1^{N_c},\dots ,{\psi }_{M_{N_c}}^{N_c},{\varphi }^1,\dots ,{\varphi }^L),$$

where L is the number of additional basis functions, which is equal to the number of ${\omega }_i:\partial {\omega }_i\cap {\mathsf{\Gamma }}_D\ne \varnothing$.

Finally, we construct a coarse grid approximation

$$S_c^n\frac{p_c^{n+1}-p_c^n}{\tau }+A_c^n{p}_c^{n+1}=F_c^n,$$

(15)

where $S_c^n=R{S}^n{R}^T$, $A_c^n=R{A}^n{R}^T$, $F_c^n=R{F}^n$ and the fine grid solution can be reconstructed, $p_{ms}=R^T{p}_c$.

5. Numerical Results

We present numerical results for 2D and 3D model problems with three types of boundary conditions on the rough top boundary. We consider the following test cases:

• Two-dimensional problem in $\mathsf{\Omega }=[1\times 0.5]$ (Figure 4) with
  • Test 1.DBC Dirichlet boundary condition on ${\mathsf{\Gamma }}_D$ with $p_1=-20.7$ and $T_{max}={10}^{-6}$.
  • Test 1.NBC Neumann boundary condition on ${\mathsf{\Gamma }}_D$ with $g={10}^4$ and $T_{max}=2·{10}^{-6}$.
  • Test 1.RBC Robin boundary condition on ${\mathsf{\Gamma }}_D$ with $\alpha ={10}^3$, $p_2=-20.7$ and $T_{max}=2·{10}^{-6}$.
• Two-dimensional problem in $\mathsf{\Omega }=[1\times 0.5]$ (Figure 5) with
  • Test 2.DBC Dirichlet boundary condition on ${\mathsf{\Gamma }}_D$ with $p_1=-20.7$ and $T_{max}={10}^{-6}$.
  • Test 2.NBC Neumann boundary condition on ${\mathsf{\Gamma }}_D$ with $g={10}^4$ and $T_{max}=2·{10}^{-6}$.
  • Test 2.RBC Robin boundary condition on ${\mathsf{\Gamma }}_D$ with $\alpha ={10}^3$, $p_2=-20.7$ and $T_{max}=2·{10}^{-6}$.
• Three-dimensional problem in $\mathsf{\Omega }=[1\times 0.5\times 0.5]$ (Figure 6) with
  • Test 3.DBC Dirichlet boundary condition on ${\mathsf{\Gamma }}_D$ with $p_1=-20.7$ and $T_{max}=2.5·{10}^{-6}$.
  • Test 3.NBC Neumann boundary condition on ${\mathsf{\Gamma }}_D$ with $g={10}^4$ and $T_{max}=2·{10}^{-6}$.
  • Test 3.RBC Robin boundary condition on ${\mathsf{\Gamma }}_D$ with $\alpha ={10}^3$, $p_2=-20.7$ and $T_{max}=2.1·{10}^{-5}$.

In Figure 4, Figure 5 and Figure 6, we depicted computational domains and heterogeneous coefficients $k_s\left(x\right)$ for these test cases. In these figures, we depict the fine grid with green color and coarse grid with blue color.

We use following fine grids and uniform coarse grids:

• Tests 1.DBC, 1.NBC, 1.RBC. Fine grid contains 13,600 vertices and 26,739 cells. Coarse grid is $10\times 5$.
• Tests 2.DBC, 2.NBC, 2.RBC. Fine grid contains $\mathrm{13,760}$ vertices and $\mathrm{26,931}$ cells. Coarse grid is $10\times 5$.
• Tests 3.DBC, 3.NBC, 3.RBC. Fine grid contains $\mathrm{43,524}$ vertices and $\mathrm{239,100}$ cells. Coarse grid is $10\times 5\times 5$.

We use $DO{F}_f$ (Degrees of Freedom) to denote fine grid system size and $DO{F}_c$ to denote the problem size of the coarse grid system using GMsFEM.

For the Haverkamp models, we use the following values of the coefficients: $\mathcal{A}=1.511·{10}^6$, $\mathcal{B}=3.96$, ${\mathsf{\Theta }}_s=0.287$, ${\mathsf{\Theta }}_r=0.075$, and $\mathcal{S}=1.175·{10}^6$. In the numerical simulations, we solve the problem for 100 time steps.

To compare the results, we use the relative $L^2$ error in %:

$$e^{L_2}=\frac{\left|\right|p-p_{ms}{\left|\right|}_{L_2}}{{\left|\right|p\left|\right|}_{L_2}}{,\left|\right|v\left|\right|}_{L_2}^2=(v,v)$$

(16)

where $p_{ms}$ and p are the multiscale and reference solutions. As a reference solution, we use a solution by the finite element method (FEM) on the fine grid. We use GMSH software [35] to construct computational domains and grids. The implementation is based on the open-source library FEniCS [36].

The fine-scale solution and multiscale solution using 16 basis functions are presented in Figure 7, Figure 8 and Figure 9 for Tests 1.DBC, 1.NBC, and 1.RBC, respectively. In

Table 1. Numerical results for Test 1.DBC. Relative $L_2$ errors (%) for different numbers of multiscale basis functions.

Number of Multiscale Basis Functions	${\mathit{DOF}}_{\mathit{c}}$	${\mathit{t}}_{25}$	${\mathit{t}}_{50}$	${\mathit{t}}_{100}$
1	76	10.271	11.519	15.397
2	142	8.679	10.183	14.026
4	274	2.308	2.162	2.388
8	538	1.163	0.883	0.885
12	802	0.779	0.542	0.509
16	1066	0.607	0.413	0.373
24	1594	0.465	0.307	0.265
32	2122	0.375	0.247	0.208

,

Table 2. Numerical results for Test 1.NBC. Relative $L_2$ error (%) for different numbers of multiscale basis functions.

Number of Multiscale Basis Functions	${\mathit{DOF}}_{\mathit{c}}$	${\mathit{t}}_{25}$	${\mathit{t}}_{50}$	${\mathit{t}}_{100}$
1	76	6.474	22.111	103.246
2	142	5.989	19.231	92.872
4	274	1.332	1.719	3.416
8	538	0.442	0.561	1.119
12	802	0.275	0.327	0.621
16	1066	0.189	0.229	0.428
24	1594	0.126	0.154	0.287
32	2122	0.105	0.123	0.221

and

Table 3. Numerical results for Test 1.RBC. Relative $L_2$ error (%) for different numbers of multiscale basis functions.

Number of Multiscale Basis Functions	${\mathit{DOF}}_{\mathit{c}}$	${\mathit{t}}_{25}$	${\mathit{t}}_{50}$	${\mathit{t}}_{100}$
1	76	7.665	21.611	83.572
2	142	7.067	18.645	63.629
4	274	1.384	1.691	2.834
8	538	0.449	0.564	0.941
12	802	0.278	0.332	0.531
16	1066	0.194	0.233	0.363
24	1594	0.131	0.155	0.242
32	2122	0.108	0.123	0.187

, we present $L_2$ relative errors for different numbers of multiscale basis functions for three time layers. We observe that the error decreases when we increase the number of basis functions for each type of boundary conditions. Therefore, to obtain a good solution, we need to take eight basis functions in each local domain ${\omega }_i$.

Next, we consider results in the 2D domain that is presented in Figure 5, where the top boundary has a more complex form. We present fine-grid and multiscale solutions for problems Tests 2.DBC, 2.NBC, 2.RBC in Figure 10, Figure 11 and Figure 12, respectively. In

Table 4. Numerical results for Test 2.DBC. Relative $L_2$ error (%) for different numbers of multiscale basis functions.

Number of Multiscale Basis Functions	${\mathit{DOF}}_{\mathit{c}}$	${\mathit{t}}_{25}$	${\mathit{t}}_{50}$	${\mathit{t}}_{100}$
1	76	11.638	13.598	18.156
2	142	9.896	11.671	15.968
4	274	2.422	2.308	2.703
8	538	1.272	0.958	1.036
12	802	0.911	0.666	0.701
16	1066	0.711	0.507	0.496
24	1594	0.526	0.378	0.356
32	2122	0.432	0.306	0.277

,

Table 5. Numerical results for Test 2.NBC. Relative $L_2$ error (%) for different numbers of multiscale basis functions.

Number of Multiscale Basis Functions	${\mathit{DOF}}_{\mathit{c}}$	${\mathit{t}}_{25}$	${\mathit{t}}_{50}$	${\mathit{t}}_{100}$
1	76	12.599	56.362	105.404
2	142	10.662	48.594	106.478
4	274	2.935	8.343	34.101
8	538	0.883	1.886	6.951
12	802	0.568	1.171	4.132
16	1066	0.383	0.724	2.328
24	1594	0.276	0.487	1.535
32	2122	0.203	0.333	1.005

and

Table 6. Numerical results for Test 2.RBC. Relative $L_2$ error (%) for different numbers of multiscale basis functions.

Number of Multiscale Basis Functions	${\mathit{DOF}}_{\mathit{c}}$	${\mathit{t}}_{25}$	${\mathit{t}}_{50}$	${\mathit{t}}_{100}$
1	76	12.907	39.998	92.513
2	142	10.889	32.043	91.271
4	274	2.607	5.064	11.835
8	538	0.841	1.324	2.765
12	802	0.551	0.827	1.671
16	1066	0.386	0.551	1.032
24	1594	0.279	0.381	0.704
32	2122	0.211	0.271	0.479

, we present the $L_2$ relative error for different numbers of multiscale basis functions for three time layers. From the obtained results, we can see that our approach can correctly take into account boundaries with such complex forms. The errors have increased slightly, but are still at a good level.

Finally, we consider the three-dimensional model problems (Tests 3.DBC, 3.NBC, 3.RBC). The fine-scale solution and multiscale solution using 32 basis functions are presented in Figure 13, Figure 14 and Figure 15. In

Table 7. Numerical results for Test 3.DBC. Relative $L_2$ error (%) for different numbers of multiscale basis functions.

Number of Multiscale Basis Functions	${\mathit{DOF}}_{\mathit{c}}$	${\mathit{t}}_{25}$	${\mathit{t}}_{50}$	${\mathit{t}}_{100}$
1	462	9.505	10.401	12.483
2	858	8.074	9.198	11.376
4	1650	4.112	4.061	5.356
8	3234	2.609	2.374	3.068
12	4818	1.411	1.253	1.556
16	6402	0.997	0.889	0.997
24	9570	0.657	0.528	0.521
32	12,738	0.469	0.339	0.324

,

Table 8. Numerical results for Test 3.NBC. Relative $L_2$ error (%) for different numbers of multiscale basis functions.

Number of Multiscale Basis Functions	${\mathit{DOF}}_{\mathit{c}}$	${\mathit{t}}_{25}$	${\mathit{t}}_{50}$	${\mathit{t}}_{100}$
1	462	8.598	26.347	54.516
2	858	7.554	24.525	52.081
4	1650	2.566	6.072	19.783
8	3234	1.292	2.594	7.792
12	4818	0.718	1.194	3.591
16	6402	0.577	0.744	2.163
24	9570	0.379	0.421	1.118
32	12,738	0.259	0.271	0.646

and

Table 9. Numerical results for Test 3.RBC. Relative $L_2$ error (%) for different numbers of multiscale basis functions.

Number of Multiscale Basis Functions	${\mathit{DOF}}_{\mathit{c}}$	${\mathit{t}}_{25}$	${\mathit{t}}_{50}$	${\mathit{t}}_{100}$
1	462	14.344	31.851	64.631
2	858	12.372	27.955	55.481
4	1650	3.834	5.904	13.551
8	3234	1.891	2.344	4.924
12	4818	0.862	1.111	2.138
16	6402	0.583	0.738	1.342
24	9570	0.339	0.409	0.691
32	12,738	0.221	0.251	0.398

, we present relative $L_2$ errors for different numbers of basis functions. We observe a good convergence of the presented method for the 3D problem. In this case, we need to take at least 12 multiscale basis functions.

From the presented results, we observe that we obtain 2–3% of relative error for four multiscale basis functions and about 1% of relative error when we increase the number of basis functions to eight for all types of boundary conditions in Test 1. For Test 2, we have similar errors for the case of Dirichlet boundary conditions, and for other types of boundary conditions, we should take more basis functions. In Test 3, we obtain accurate results for 12 multiscale basis functions. In addition, for Robyn-type boundary conditions, we have about 2% of errors when we take 4, 8, and 12 basis functions in Tests 1, 2, and 3, respectively. We note that the presented results are confirmed with the analytical error estimation presented in [27,34], where authors showed that error depends on the number of basis functions and on eigenvalues of the local spectral problems. We also mention that an adaptive approach can be applied, where, to find an optimal number of basis functions, we can use local eigenvalues in order to capture all local characteristics of the solution [34].

From the results, we can see that we can get a solution with a small error when using a small number of degrees of freedom. For example, from Table 9, we can see that we need to use 4818 degrees of freedom, which is $11\%$ of the size of fine-grid system, to get very accurate results. We can observe a decrease in the error with an increase in the number of basis functions. Error reduction is associated with the theory of the method itself; for more details, see [27,34].

6. Conclusions

In conclusion, a multiscale method for the unsaturated flow problem in heterogeneous porous media with rough boundaries on the top is presented. On the top boundary, we consider three different types of boundary conditions: Dirichlet, Neumann, and Robyn boundary conditions. We perform a fine grid approximation using the finite element method. To reduce the size of the discrete system, we construct a coarse grid approximation using GMsFEM. To take into account rough boundaries with non-homogeneous boundary conditions, we construct an additional basis function for each type of boundary conditions. Numerical results are presented for two- and three-dimensional test problems. Numerical comparisons of the relative error for different numbers of basis functions are presented. The proposed methods provide a good accuracy for all types of boundary conditions. Note that this is the first step to being validated on the synthetic model. In our future work, we plan to proceed with more realistic cases.