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

: In this paper, we consider unsaturated ﬁltration in heterogeneous porous media with rough surface topography. The surface topography plays an important role in determining the ﬂow 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 inﬂuence 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 (ﬁne-grid) solutions for different numbers of multiscale basis functions. We obtain a good agreement between ﬁne-grid and coarse-grid solutions.


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

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 Ω, with the initial condition In the above equation, Θ 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 Θ are related by the following constitutive relations (Havercamp model): where k s (x) is a heterogeneous coefficient modeling the saturated hydraulic conductivity and A, S, B, Θ r , Θ s are the Haverkamp model coefficients.
We consider domains with rough boundaries. To be specific, we consider a domain Ω, illustrated in Figure 1. We consider three types of boundary conditions on the top boundary Γ D : 1. Dirichlet boundary condition: 2. Neumann boundary condition: 3. Robin boundary condition: On the other part of the boundary, denoted by Γ N , we impose the zero flux boundary condition

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: where ∂Θ(p)/∂t = C(p)∂p/∂t, C = dΘ/dp, C n = C(p n ), and k n = k(x, p n ). In (4), the number τ > 0 denotes the time step size, and the superscript n denotes the value at the time t = nτ. For instance, p n denotes the value of the solution p at the time t = nτ.
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: Let T h be the fine grid with mesh size h for the domain Ω. We write the approximate solution as follows where N f is the number of interior fine grid vertices for (5) (or number of fine grid vertices for (6) and (7)) and {ψ i } are the conforming piecewise linear basis functions. Using p h = (p h,1 , p h,2 , ..., 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: where The following summarize the definitions of n and F n for each of the considered boundary conditions: 1. Dirichlet boundary condition:

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 T H = ∪ i K i be the coarse grid for computational domain Ω with size H and K i be the i-th coarse grid cell. We define a local domain ω i as the union of coarse grid cells around one coarse grid vertex, i = 1, ..., 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 ω i , i = 1, ..., 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 ψ i k 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 ω i snap for each local domain ω i . The snapshot space is constructed by solving the following local problems: where δ 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, ..., J i , J i is number of fine grid nodes on boundary ∂ω i /Γ D ). Then, we can define the snapshot space and the projection matrix on the snapshot space in local domain ω i as follows: Note that we do not calculate snapshots for boundary Γ 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 ω i :ÃΨ We compute the solution of the spectral problem by Ψ i j = R snap,iΨ i snap,j . We note that for computation of the multiscale basis function, we use only the linear part k s (x) of coefficient k(x, p) from (3).
Next, we select the smallest M i eigenvalues and use them for multiscale basis function construction (Ψ i j , j = 1, ..., M i ). We obtain the required multiscale basis functions after multiplication by the linear partition of unity function where χ i is the standard coarse grid nodal basis function for the coarse node i in local domain ω i . To handle rough boundary Γ D with non-homogeneous boundary conditions, we calculate an additional basis function in local domains ω i : ∂ω i ∩ Γ D = ∅. We solve the following local problem:

Robin boundary condition
To compute additional basis functions, we multiply the solution of the local problem to the linear partition of unity function The solutions of problems (12)- (14) are presented in Figure 3. Finally, we construct the multiscale space where S n c = RS n R T , A n c = RA n R T , F n c = RF n and the fine grid solution can be reconstructed, p ms = R T p c .
In Figures 4-6, we depicted computational domains and heterogeneous coefficients k s (x) for these test cases. In these figures, we depict the fine grid with green color and coarse grid with blue color.   We use DOF f (Degrees of Freedom) to denote fine grid system size and DOF 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: A = 1.511 · 10 6 , B = 3.96, Θ s = 0.287, Θ r = 0.075, and 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 %: 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 Figures 7-9 for Tests 1.DBC, 1.NBC, and 1.RBC, respectively. In Tables 1-3, 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 ω 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 Figures 10-12, respectively. In Tables 4-6, 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.       Tables 7-9, 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].

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.