An Online Generalized Multiscale Finite Element Method for Unsaturated Filtration Problem in Fractured Media

: In this paper, we present a multiscale model reduction technique for unsaturated ﬁltration problem in fractured porous media using an Online Generalized Multiscale ﬁnite element method. The ﬂow problem in unsaturated soils is described by the Richards equation. To approximate fractures we use the Discrete Fracture Model (DFM). Complex geometric features of the computational domain requires the construction of a ﬁne grid that explicitly resolves the heterogeneities such as fractures. This approach leads to systems with a large number of unknowns, which require large computational costs. In order to develop a more efﬁcient numerical scheme, we propose a model reduction procedure based on the Generalized Multiscale Finite element method (GMsFEM). The GMsFEM allows solving such problems on a very coarse grid using basis functions that can capture heterogeneities. In the GMsFEM, there are ofﬂine and online stages. In the ofﬂine stage, we construct snapshot spaces and solve local spectral problems to obtain multiscale basis functions. These spectral problems are deﬁned in the snapshot space in each local domain. To improve the accuracy of the method, we add online basis functions in the online stage. The construction of the online basis functions is based on the local residuals. The use of online bases will allow us to get a signiﬁcant improvement in the accuracy of the method. We present results with different number of ofﬂine and online multisacle basis functions. We compare all results with reference solution. Our results show that the proposed method is able to achieve high accuracy with a small computational cost.


Introduction
In this paper, we consider an unsaturated filtration problem in fractured heterogeneous media. For unsaturated filtration, we formulate a mathematical model that is based on the Richards' equation [1][2][3][4][5]. Due to the high permeability of the fractures, these fractures have a significant effect on the flow processes and require a special approach in the construction of a mathematical model and computational algorithms. Such problems in fractured and heterogeneous media require very fine grids. Standard approaches use the finite element method to accurately model these processes [6][7][8]. Such approaches lead to large systems with large numbers of unknowns, and solving them is computationally expensive. Multiscale methods are designed to reduce the size of the system. A multiscale model reduction technique is based on the construction of multiscale basis functions to extract information in the micro-level. These basis functions are then used to obtain coarse scale equations.
In previous works the multiscale model reduction of Richards equation are presented, using offline multiscale basis functions [26][27][28]40]. In this paper we use an online generalized multiscale finite element method. This method consists of offline and online stages. In the offline stage, we obtain multiscale basis functions by solving spectral problems in local domains. In particular, we will solve spectral problems on the snapshot space. The snapshot space helps us to take into account some complex properties of computational domain. Using snapshot space is justified for high-contrast heterogeneous domains, and for complex domains that contains fractures and channel. Using multiscale basis function we obtain an offline space on which we will solve our problem. In the online stage we solve our problem on the coarse grid using the offline space. For complex heterogenous problems, there is a need to enrich the approximation space using online basis functions. These online basis functions are computed adaptively in the online stage using local residuals. The aim is to reduce the error significantly by using a small number of these online basis functions. The motivation of the construction of online basis functions is the optimize the error reduction locally, and mathematical theories show that the online basis functions give fast convergence [41].
We will present numerical results for some two-dimensional examples. We consider two cases of fractured media. In the first case we take domains with simple network of fractures, and in the second case we test our method in domain with complex fracture network. We expect that, online basis functions will allow to us modeling unsaturated filtration process in all types of fracture geometry.
We organize the paper as follows. In Section 2, we present a mathematical model for unsaturated filtration problem in heterogeneous fractured media. In Section 3, we consider fine-scale approximation using a Discrete Fracture Model. Next in Section 4, we describe the construction of the coarse grid approximation and describe construction of the multiscale basis functions. In this Section, we present an offline multiscale basis functions and an online multiscale basis functions. We present numerical results in Section 5 for two-dimensional model problem. Finally, we present the Conclusion.

Problem Formulation
We use Richard's equations to describe an unsaturated filtration process in fractured porous media. Let Ω is the domain of the porous matrix. We consider lower dimensional fractures due to small thickness of the fractures compared to the domain sizes. We solve the following problem of the unsaturated filtration in fractured porous media with following initial condition and boundary conditions where To define the physical properties of the domain Ω we use Havercamp model [3]. Then, the non-linear coefficient k q (x, p q ) and the water content Θ q are determined as follows where k q,s (x) is the saturated hydraulic conductivity, α q β q , γ q , A α , Θ α,r , Θ α,s are the Haverkamp model coefficients and q = m, f .

Fine Grid Approximation
To perform numerical experiments, we solve a reference solution on a fine grid. We construct unstructured triangular fine grid that explicitly resolve fractures in the level of mesh. We construct an approximation on the fine grid by the finite element method and Discrete Fracture Model (DFM) for fracture networks. We have the weak formulation For the approximation of the time derivative, we use the backward Euler method. To resolve non-linearity we use simplified approximation from the previous time step and obtain the following approximation on the fine grid where We can write approximation, in the matrix form as To approximate lower dimensional fractures we use DFM approach. We assume that p m = p f and using superposition principle [35,43] we eliminate p f from Equation (7) and obtain following system where S n = S n m + S n f , A n = A n m + A n f and F n = F n m + F n f .

Coarse Grid Approximation
We construct a coarse grid approximation using an online generalized multiscale finite element method. Let T H is the coarse grid and ω i is the local domains, where i = 1, . . . , N c and N c is the number of coarse grid nodes. A local domain ω i is obtained by the combining all the coarse cells around one vertex of the coarse grid.
The Online GMsFEM procedure consist of two parts (see Figure 1) • Offline stage. In the offline stage we define an offline space by constructing an offline multiscale basis function; • Online stage. In the online stage we construct the system on the offline space and enrich offline space by online multiscale basis functions.

Offline Stage
We start from the constructing a snapshot space V ω i snap . The snapshot space are obtained by solving next local problems with boundary condition φ i m,j = δ j , φ i f ,j = δ j , where δ j is the discrete delta function which 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 ). The main concept of constructing the snapshot space is that the snapshot vectors preserve some essential properties of the solution and provide 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. For the construction of the accurate approximation on the coarse grid for problems in fractured media, the number of dominant modes in spectral problem that corresponds to the very small eigenvalues should belong to the long fractures that cross the local domain boundaries. For the small separated fractures that do not cross the boundaries only one basis is sufficient for the approximation. When we skip construction of a snapshot space, a separate basis will be built for each small separate fracture. This approach will require more multiscale basis functions in each local domain ω i .
The snapshot space and the projection matrix on the snapshot space are defined as follows On the snapshot space we solve the next local spectral problem in each ω i to obtain an offline multiscale basis functionsÃΨ We use solution of the spectral problem Ψ i j = R snap,iΨ i snap,j in offline basis construction. We take only linear part k s (x) of coefficient k(x, p) from (4). For basis construction we choose the smallest M i eigenvalues. We obtain an offline multiscale basis functions by multiplication on the linear partition of unity function ψ i j = χ i Ψ i j , where χ i is the linear coarse grid nodal basis function that is equal to zero on the boundary of local domain ω i and one at the coarse grid node i. An example of first 4 solution of spectral problem and first 4 offline multiscale basis functions are presented in Figure 2.

Online Stage
Since our problem is non-linear, we need to take into account the change in the coefficients over time. It is possible to consider ways of updating offline basis functions at certain time steps, but this approach has large computational costs. It would be much more profitable to add one basis to each local area, which minimizes the residual. This way is based on the online basis enrichment for GMsFEM based on the residual [41]. Presented procedure will allow us to significantly reduce the number of used offline multiscale basis functions and add local online multiscale basis functions to fast error reduction based on the local residual information. With fewer basis functions, we will be able to get better accuracy, since our bases will take into account the coefficient k changes over time.
In online GMsFEM, we update projection matrix R n by adding online residual based multiscale basis functions at n-th time step. Therefore, we solve next system on the coarse grid where coarse scale matrices and vectors are constructed using current projection matrix R n Note that, when we use pre-constructed offline multiscale basis functions for given heterogeneity and fracture distribution in coarse scale system construction, we can use predefined projection matrix without updating it at n-th time step In online GMsFEM, we set R n = R at first time step (n = 0), then we construct and solve coarse scale system (12). Next, at time step n = 1, 2, . . ., we enrich multiscale space by residual based online multiscale basis functions. In order to enrich space, we solve system (12) with R n = R n−1 , then we calculate online multiscale basis functions ϑ i 1 locally in ω i using current residuals and update the projection matrix (R 1,n ) T = (ψ 1 1 , . . . , ψ 1 M 1 , . . . , ψ N c 1 , . . . , ψ N c M Nc , ϑ 1 1 , . . . , ϑ N c 1 ).
We can reiterate the process with residual calculation and add more online basis functions on the n-th time step where L is the number of online iterations. We can enrich multiscale space not for every time steps, for example, every 5-th or 10-th. Then, when we keep R n for the time steps where we do not want to update online basis functions. Next, we present construction of the online multiscale basis functions in local domain ω i in details. Construction of the local residual based online multiscale basis functions is based on the solution of the following local problem in ω i : where with zero Dirichlet boundary condition Φ i l = 0, Φ f ,i l = 0 on the ∂ω i . Finally, we obtain an online multiscale basis functions by multiplication on the linear partition of unity function ϑ i l = χ i Φ i l . We present solution of (13) and (14) for three different local domains in Figure 3 (l = 1).
where coarse scale matrices and vectors are constructed using current projection matrix R l,n S l,n c = R l,n S n (R l,n ) T , A l,n c = R l,n A n (R l,n ) T , F l,n c = R l,n F n , p l,n+1 with projection matrix where l = 1, . . . , L, R n = R L,n and L is the number of the online iterations.

Numerical Results
We present numerical results for two-dimensional problem with in fractured heterogeneous media. We consider problems in two-dimensional domains with different fractures locations. In this work, we consider the following model problems:   To compare the results with reference solution, we use the relative L 2 error (%) and relative energy error (%) (16).
where p ms and p are the multiscale and reference (fine-scale) solutions.
In Test 1 we consider domain with 14 fractures with homogeneous porous matrix. In Figure 5 we present fine grid solution, multiscale solution using 3 offline basis functions and multiscale solution using 3 offline and 1 online basis functions at final time. We perform an enrichment procedure in each fifth time layer. In Table 1, we present relative L 2 and energy errors for multiscale solution using 0, 1, and 2 online basis functions for different number of offline basis functions. In this table DOF c denotes the vector size on the coarse grid and M denotes the number of offline basis function in each local domain ω i . By this tables we can see that, when we use 3 offline and 1 online basis functions we can obtain a solution better then using 8 offline basis function. Using online basis functions, we can obtain a solution using system with smaller number of unknowns. We can see the advantage of the online approach over offline in the Figure 6. In this Figure, we depict a relative L 2 and energy error over time. The Figure 6 clearly shows the fall in error every fifth time layer, when we update the online bases. In this Figure, we observe a jump in error when we use 2 online basis functions. From the results we can see, that the method give solution with high accuracy. It is most profitable to use 1 online multiscale basis function in each local domain ω i . When using the 2-nd online basis, the fall in error is not so significant, and its addition is not worth the spent computing resources. As we can see, we get a significant reduction in the size of the original system with very little loss of accuracy.   Next, we consider Test 2 with heterogeneous matrix. In Figure 7 we present fine grid solution, multiscale solution using 3 offline basis functions and multiscale solution using 3 offline and 1 online basis functions at final time. In Table 2, we present relative L 2 and energy errors for multiscale solution using 0, 1, and 2 online basis functions for different number of offline basis functions. In this case, we also enough to use 3 offline and 1 online multiscale basis functions. We also show the advantage of the online approach over offline in the Figure 8. Here we do not observe any jumps in errors. Adding a heterogeneous matrix did't affect on the accuracy of the method. The behavior of the method is similar to the previous task (Test 1). Therefore, we can draw a similar conclusion that method provides a good solution.   Next, we consider Test 3, the case with 50 fractures. In Figure 9, we present fine grid solution, multiscale solution using 12 offline basis functions and multiscale solution using 12 offline and 1 online basis functions at final time. In this case, the enrichment procedure are performed also on each fifth time layer. In Table 3 we present relative L 2 and energy errors for multiscale solution using 0, 1, and 2 online basis functions for different number of offline basis functions and adaptive approach. From results, we can see that we need to use at least 12 offline and 1 online basis functions to obtain solution with good accuracy. Using the 2-nd online basis function improves accuracy of the method but not good enough to use it. The behavior of the obtained results is the same as in the Test 1 and Test 2, but in this case we need to use a larger number of offline bases. In Figure 10, we present a relative L 2 and energy error over time. In this Figure, we can see the error raising when we use 2 online basis functions. These results once again confirm that it is better to use 1 online multiscale basis function.

Conclusions
We presented an online generalized multiscale finite element method for unsaturated filtration problem in fractured media. We performed a multiscale modeling for domains with 14 and 50 fractures for homogeneous or heterogeneous matrix. We considered a multiscale method with different numbers of offline and online basis functions. Our experiments showed that adding 1 or 2 online multiscale basis function significantly reduces the error of the method, especially in energy norm. Moreover, we considered an adaptive online approach, namely, online basis functions are added at some selected regions with larger errors. This approach performed well in all experiments. This method allows us to obtain a solution with high accuracy with a significant decrease in the dimension of the original system. We conclude that for solving the unsaturated filtration problems in fractured media by the online generalized multiscale finite element method, it is sufficient to use one online multiscale basis function in each local domain.