Abstract
We propose a fully discrete method for the multiscale Richards’ equation of van Genuchten-Mualem model which describes the flow transport in unsaturated heterogenous porous media. Under the framework of heterogeneous multiscale method (HMM), a fully discrete scheme combined with a regularized procedure is proposed. Including the numerical integration, the discretization is given by piecewise finite element in space and an implicit scheme in time. Error estimates between the numerical solution and the solution of homogenized problem are derived under the assumption that the permeability is periodic. Numerical experiments with periodic and random permeability are carried out for the van Genuchten-Mualem model of Richards’ equation to show the efficiency and accuracy of the proposed method.
1. Physical Model
Richards’ equation is most often used to model the movement of groundwater flow in saturated-unsaturated porous media. It was formulated by Richards in 1931 [1]. It is a nonlinear partial differential parabolic equation. Depending on the saturation and pressure, three main forms of Richards’ equation are usually presented: pressure-based form, saturation-based form, or mixed form. Combining the continuity equation with Darcy’s law, the mixed form can be expressed as where the reduced saturation is defined by . (the residual fluid content) and (the fluid content at saturation) are constants relying on the porous medium. denotes the saturation and denotes the pressure. Due to the complex heterogeneity of natural medium, the permeability of the medium oscillates rapidly with large contrast. is the characteristic length presenting the small scale variability of the media. The coordinate in the direction of gravity is denoted by .
In this paper, we adopt the retention curves proposed by van Genuchten [2] and Mualem [3]. The relations of the reduced saturation , pressure , and permeability are read as where , and are the parameters of porous medium. and can be presented as the absolute and relative permeability, respectively. When , , the medium is in the saturated state; otherwise, it is in the unsaturated case.
Under the unsaturated case, Richards’ equation can also be rewritten in terms of reduced saturation by using the above relations. Consider with the moisture diffusivity defined by Here, under the relations (2), we should point out that (1) can model both saturated and unsaturated cases. In saturated region, (1) is an elliptic equation; in unsaturated region, the equation is parabolic. On the other hand, in unsaturated case, (1) can be degenerated because may approximate to zero when (in almost completely dry region). Richards’ equation (3) can only model the unsaturated case. But (3) may also be degenerated, because the diffusivity may vanish or explode (see Figure 1(a)). For (in almost completely dry region), the diffusivity vanishes, while for going to 1 (in almost full saturated area) goes to infinity. The degeneracy of (1) and (3) leads to the fact that the solutions must be understood in the sense of distribution as proposed in [4–7].
(a)
(b)
Numerous papers have been published on numerical scheme for the non-multiscale degenerate parabolic problem (including Richards’ equation). For example, in [5–9], the authors firstly regularized the degenerated problem then used the FEM or mixed FEM for spatial discretization. In [10–13], considering the time discretization aspect, the authors developed some schemes for time discretization, such as linear time discretization, adaptive time stepping, or relaxation scheme. The relaxation-iteration scheme analyzed in [13] for the fast diffusion case is a fixed point iteration which was proposed in [8] (also see [12] for the mixed finite element context).
Here, we should point out the work of [6, 7]. In [6], the authors considered a class of degenerate parabolic equations including Stefan problem, porous-medium problem, and nonstationary filtration problem. Combining with a regularization approach (when necessary), the authors developed their fully discrete scheme including FEM in spatial and Backward-Euler semi-implicit scheme in time. At the same time, the effects of numerical integration and domain change were taken into account. In [7], under the relations (2), Backward-Euler implicit scheme in time with a regularization step for non-multiscale Richards’ equation (3) was established and analyzed.
On the other hand, for this kind of multiscale problem, it is impossible to account explicitly for the spatial variability at fine scale because of the computational resource limitations in realistic situation. So, traditional numerical methods such as standard finite element methods (FEM) and finite difference methods (FDM) are generally not capable of solving this multiscale problem directly. In recent years, a number of multiscale numerical methods, such as multiscale finite element method (MsFEM) [14], heterogeneous multiscale method (HMM) [15], and upscaling method [16], have been proposed to solve the general multiscale problems based on similar ideas. Among them, the so-called heterogeneous multiscale method (HMM) has proved to be an efficient tool to assemble information from microscale problems in order to perform macroscale simulations. In [17], the author discussed the modeling and analysis of finite element methods (such as hybrid methods, coupling spectral or discontinuous Galerkin methods with FEM) for multiscale problems constructed in the framework of the HMM for multiscale PDEs (such as elliptic and parabolic equations).
To our knowledge, there are few papers about the numerical analysis of multiscale Richards’ equation of van Genuchten-Mualem model. Based on the frame of HMM-FEM, we will use the idea of [6, 7] to develop a fully discrete method for multiscale Richards’ equation (5). In order to overcome the above difficulties, we regularize the equation firstly; secondly, a fully discrete multiscale numerical method based on heterogeneous multiscale method (HMM) [15] and FEM is developed on a macroscale mesh; at the same time, we use the numerical quadrature for computing the integral over each element. In the paper, we are only interested in (3).
Noticing that is separable, we set (see Figure 1(b)). By Kirchhoff transform [4], we rewrite (3) as with . So, when , (5) may degenerate to hyperbolic equation; on the other hand, , and (5) may transfer to elliptic equation.
Notations. is a bounded polyhedral domain with boundary . Set , where is fixed. For defining a solution in weak sense, we let stand for the inner product on or the duality pairing between and , stand for the norm in , and and stand for the norm in and .
2. Formulation of the Problem
In this paper, we will consider the following general nonlinear parabolic equation which includes Richards’ equation (5):
Assumption 1. (H1) is a maximal monotone graph in and .
(H2) For all , and has the following asymptotic behavior:
(H3) almost everywhere.
(H4) The matrix is bounded and uniformly elliptic. There exist constants and such that , .
(H5) is bounded, uniformly Lipschitz continuous and satisfying
A weak solution of problem (6) is defined as follows.
Definition 2. Find such that and and, ,
Remark 3. The existence, uniqueness, and regularity results for the above problem can be found in [4]. The solution of problem (6) may denote the reduced saturation and therefore should be bounded by 0 and 1. In the setting stated above, the maximum principle holds for problem (6), so the solution remains bounded by 0 and 1.
Remark 4. According to the constitution relations of the van Genuchten-Mualem, the parameter in (H2) is determined as and (). In Richards’ equation, the convective function denotes the permeability . And it is easy to verify that and when approaches to zero. Also, it is easy to verify that condition (8) is valid when the second variable approaches to zero. The other cases are obvious.
Let be a small perturbed parameter. We approximate by, for all , So, we get Hence, is also a maximal monotone graph in and it follows that
Let be a regular triangle decomposition of , where stands for the mesh size. Define the functional space
In order to simplify the computation, we will use the numerical integration scheme as in [6]. Let be the local linear interpolant operator; then, for any piecewise uniformly continuous functions and . So, for all , and every element , the numerical integration scheme satisfies that [6] and, for any [18], stands for the piecewise linear interpolant operator and satisfies and stands for the discrete projection operator; more precisely, for any , and satisfies Let and () denote the quadrature nodes and weights in . Let be an integer and let . Our fully discrete formulation for problem (6) is read as follows: find , (), such that where and the operator is defined by the following problem: where and the cell is a square of size centered at . So far, our fully discrete HMM-FEM is settled down and the porous media do not have to be periodic.
3. Some Results on the Homogenized Problem
In this section, we will review the results of homogenization and the analysis of the discrete homogenized problem. Under the assumptions (H1)–(H5) and if the function and are periodic in with unit square . In [19], the authors have established the homogenization theory for (6). It has been shown that, in suitable topology space, the solutions of problem (6) converge to the solution of the following problem as . where and are defined as and and are the periodic solutions of following cell problems, respectively. and and inherit the property of and [16].
Remark 5. According to Alt and Luckhaus [4], we have at least and . Moreover, the maximum principle leads to . We therefore conclude that . This gives us pointwise for every .
The problem (26) can be read as follows.
Definition 6. For every time interval , find such that and, for all ,
We also give out the fully discrete scheme for the problem (26) by a regularization procedure.
Definition 7. Find such that and, for all , where .
In the rest of this section, we will give out a prior estimate for problem (31) and the error estimate of and .
Theorem 8. Under the assumptions (H1) to (H5), there exist a positive constant independent of such that
Proof. Let in (31) and sum it over from 1 to .
We will estimate each term separately. Following the idea in [6], we have
Using (H4), (H5) and Poincaré inequality, we have
Since , we have . Combining all the terms and choosing properly, we complete the first part.
Again, we let in (31) and sum it over from 1 to .
Considering that , we have
By (H4), (H5), and (32), it is easy to get
So, we finish the second part of the theorem.
Define a regular G-operator and -operator as in [6], for all , and there is Firstly, we have On the other hand, so we have equivalent to . We also introduce the following two identities:
Theorem 9. and are the solutions of problems (26) and (31) at time , respectively. Denote that , , and . Under the assumptions (H1) to (H5), there exists a positive constant independent of , , , such that
Remark 10. The proof method of Theorem 9 can be found in [6, 12]. However, there is a little difference between our result and the result in [6, 12] because we consider the two degenerate points case.
4. Main Result
In this section, the main task is to give out the error between and . Here, and are the solutions of problems (31) and (22), respectively. Before the proof of our main result, the following useful lemma will be introduced [15, 16].
Lemma 11. , are defined in (27) and (28) and , are defined in (23) and (24), respectively. Then, for all ,
Then, we have the main result of the paper.
Theorem 12. , and , are the solutions of problems (22) and (31) at time , respectively. Denote and . Then, under assumptions (H1)–(H5), there exists a positive constant independent of , , , such that
Proof. Subtract (31) from (22) and notice the notation ; we have
that is,
Let and sum ; then we have
Denote the above equality by .
For the term , noticing that is also positive and bounded, we use (15), (45), and the definition of to get
Rewrite the term as
Applying (20), (12), and (32), we get
For the terms , , we have used (19), (18), and Lemma 11.
Similar to , the term can be estimated as the following:
Noticing (also ) and (H5), it follows that
So, we have
Combining all the terms and choosing the parameter properly, we have
At last, we use the Gronwall inequality to the above formulation to finish the proof completely.
5. Numerical Example
In this subsection, we consider the Richards equation under the van Genuchten-Mualem model [2]. Consider and, here, the constitutive relations are where and , , , and are the media parameters. Throughout this subsection, we set , , , and . The heterogeneity comes from absolute permeability and .
Notice that the Richards equation under the van Genuchten-Mualem model (59) is degenerated at () and () and sharp interface would be developed between saturated and unsaturated regions. Our HMM-FEM scheme can be also used to treat this kind of problems.
In our computation, in order to avoid the difficulty of degeneration, the coefficient is replaced by a regularized one as where .
5.1. van Genuchten-Mualem Model with Periodic Coefficients
For periodic case, we set In our simulation, is chosen to be .
The numerical solution of HMM-FEM is obtained on a macroscale grid . For each local sample cell , we choose and solve the microproblems (25) on a grid with size of . We compare the HMM-FEM solution with the fine-scale solution which is obtained on a grid by FVM. The time step is also chosen to be in this subsection. We compare the HMM-FEM solution with the fine-scale solution along the cross-section of (Figure 2(a)).
(a)
(b)
5.2. van Genuchten-Mualem Model with Random Coefficients
For the random model, we only consider the isotropic heterogeneities case. We generate the random log-normal permeability fields and by the same methods used in [16]. The corresponding normal distribution of is and the corresponding normal distribution of is . The correlation lengths of both and are in and directions. The solution of HMM-FEM is obtained on a macroscale grid . For each local sample cell , we choose and solve the microproblems (25) on a grid with size of . We compare the HMM-FEM solution with the fine-scale solution which is obtained by solving (59) on a grid by FVM. We compare the HMM-FEM solution with the fine-scale solution along the cross-section (Figure 2(b)).
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The research is supported by the Fundamental Research Funds for the Central Universities (no. 2013B10114), the Natural Science Foundation of Jiangxi Province (no. 20132BAB211018), and the NSF of China (no. 11271281).