New Scheme of Finite Difference Heterogeneous Multiscale Method to Solve Saturated Flow in Porous Media

Abstract

A new finite difference scheme, the development of the finite difference heterogeneous multiscale method (FDHMM), is constructed for simulating saturated water flow in random porous media. In the discretization framework of FDHMM, we follow some ideas from the multiscale finite element method and construct basic microscopic elliptic models. Tests on a variety of numerical experiments show that, in the case that only about a half of the information of the whole microstructure is used, the constructed scheme gives better accuracy at a much lower computational time than FDHMM for the problem of aquifer response to sudden change in reservoir level and gives comparable accuracy at a much lower computational time than FDHMM for the weak drawdown problem.

1. Introduction

Natural porous media exhibit a significant spatial variability in most attributes of hydrogeological interest. For instance, it is quite typical for hydraulic conductivity to vary orders of magnitude over distances [1]. The groundwater flow problems in heterogeneous porous media can be accurately solved by using conventional finite element method or finite difference method based on smaller scale, which leads to more computational cost. Discrete schemes obtained in this way are often by far too expensive to be solved directly. For the sake of the accuracy and efficiency, several different but related multiscale methods, such as the multiscale finite element method (MsFEM) [2, 3], the heterogeneous multiscale method (HMM) [4], and the numerical homogenization method [5], for problems with oscillating coefficients have been proposed to accommodate the fine-scale description directly. Here, we should also mention the work of Babuška in the 70s [6–8], which motivated these multiscale methods in an extent.

Multiscale solution methods are currently under active investigation for the simulation of subsurface flow in heterogeneous formations [9]. Ye et al. [10] applied MsFEM to simulate two- and three-dimensional saturated flow problems. Chen and Hou [11] proposed a mixed multiscale finite element method for elliptic problems with oscillating coefficients; they demonstrated the efficiency and accuracy of the proposed method for flow transport in a porous medium with a random log-normal relative permeability. He and Ren [12] presented the finite volume MsFEM for solving saturated flow in heterogeneous porous media. E et al. [13] took a systemic analysis of HMM for elliptic homogenization problems, where the error between the numerical solution of HMM and the solution of homogenized equation is estimated, and how to construct better approximation of the exact solution from the HMM solution is discussed. Ming and Zhang [14] applied HMM to the linear parabolic homogenization problem and discriminated different types of microscopic models. Ming and Yue [15] discussed the numerical performance of HMM including comparison with other methods. Yue and E [16] developed HMM for linear and nonlinear transport equations with multiscale velocity fields in heterogeneous porous media and focused on the problem where advection is dominant at the small scale.

Most of existing multiscale methods have been limited to the finite element method [17–19]. There are also widely used finite difference flow and solute transport models in both the groundwater and oil industries. To handle multiscale problems with finite difference methods, based on the framework of HMM [4], Abdulle and E [20] proposed the finite difference heterogeneous multiscale method (FDHMM) for solving multiscale parabolic problems. This method includes a “heterogeneous” discretization which cares about the fine scale only on small representative region of the spatial domain. FDHMM has two components: a macroscopic scheme evolved on a coarse grid (the grid of interest) with unknown data recovered from the solutions of the microscopic model and a microscopic scheme, in which the original equation is solved on a sparse (heterogeneous) spatial domain. The similar idea can be found in [21]. Chen and Ren [22] applied FDHMM to Richards’ equation; they found that FDHMM could effectively simulate the transient unsaturated flow in the specific soils. In the saturated flow case, we may precompute the macroscopic flux at a preprocessing step to save the computational time. In addition, Abdulle and E [20] studied the multiscale parabolic equation without the source sink term and considered the examples where the coefficients change according to the smooth function. For transient flow problem in heterogeneous porous media, the coefficients generally change in a random form; thus there is a need for more evaluation of the applicability of FDHMM.

Here, we propose a new scheme of FDHMM for simulating not only the steady saturated flow problem but also the transient saturated flow problem in geostatistical random porous media. The constructed scheme employs an idea presented by Ming and Zhang [14]; that is, the microscale parabolic model can be reduced to the microscopic elliptic model for the problem without oscillation propagation in time. Motivated by the construction of the multiscale finite element base functions [2, 3], in every control volume, we divide the microscopic elliptic model into two basic microscopic elliptic models and estimate the basic macroscopic flux by the solutions of these two basic microscopic elliptic models. The small scale information is then brought to the large scale through the approximation of basic macroscopic fluxes. These basic macroscopic fluxes are just calculated once at the preprocessing step and will be used in the subsequent computations. In general, governing coarse-grid equation and coupling the approach of the new scheme are the same as those of FDHMM by Abdulle and E. The main difference between the two methods is the microscopic scheme, in contrast to FDHMM by Abdulle and E, where the numerical fluxes are computed on the fly using localized and more resolved computations which means that FDHMM by Abdulle and E needs the macroscopic and microscopic evolution at every time step and the new scheme adopts the idea of MsFEM of Hou et al., by numerically precomputing a finite difference analogue of a multiscale shape function, which provides a fixed expression for the numerical basic flux in terms of the coarse variables. It means that the fine-scale information is coupled into the coarse scale by this finite difference analogue of a multiscale shape function. In addition, the new scheme incorporates ideas from Ming and Zhang [14] to transform the microscopic parabolic model to a microscopic elliptic model, which allows MsFEM ideas to be adopted in the computational scheme.

Our method is also analogous to the classical upscaling method, where the upscaled hydraulic conductivities are precomputed [23, 24]. Different from the classical upscaling method, the present method only precomputes the basic macroscopic fluxes. The estimation of the macroscopic fluxes, which contain both microscopic information of the medium property and useful information about the gradients of the solutions of microscopic elliptic models, is coupled into the course of solving the coarse equation, and it makes the constructed scheme put more emphasis on the interaction between the macro- and microscale behavior. On the other hand, in the new scheme, the fine-scale global flow solution is decomposed into a series of local microscopic problems; the computations of these basic microscale problems can be carried out sequentially; this obviously saves the computational time and the memory requirement, which may provide the proposed method with a possibility to solve large flow problems under restricted computational capabilities.

This paper is organized as follows. We firstly describe the flow problem and introduce the principle and the algorithm of the constructed scheme in detail. Numerical examples to illustrate the performance of the constructed scheme are arranged in Section 3. Some conclusions are given in Section 4.

2. New Scheme of FDHMM

2.1. Flow Problem

The transient saturated flow through a heterogeneous porous medium is governed by the parabolic partial differential equation 𝑆_𝑠(𝐱)𝜕𝜓(𝐱,𝑡)𝐊(𝐱)((),(),𝜕𝑡=∇⋅∇𝜓𝐱,𝑡)+𝑅𝐱,𝑡inΩ×0,𝑇 (1) where 𝜓 is the hydraulic head, 𝑆_𝑠 is the specific storage coefficient, 𝐊 is the hydraulic conductivity tensor, 𝑅 is the source sink term, 𝐱=(𝑥,𝑦) is the spatial coordinate, 𝑡 is the time variable, Ω is the study area, and 𝑇 is the time domain.

2.2. Principle and Algorithm

The discretization in this study is the mesh-centered finite difference. To simplify the presentation of the constructed scheme, we assume that the solution domain Ω is a square and uniformly discretize it with a coarse 𝑁×𝑁 mesh. Let the Cartesian coordinates of this coarse mesh be represented by (𝑥_𝑖,𝑦_𝑗), 𝑖,𝑗=1,…,𝑁+1. 𝐻=𝑥_𝑖+1−𝑥_𝑖=𝑦_𝑗+1−𝑦_𝑗 denotes the coarse mesh size.

Notice that a macroscopic model is known to exist according to the homogenization theory, and the idea of the constructed scheme is to evolve a macroscopic model for the flux form of (1), 𝑆_𝑠(𝐱)𝜕Ψ(𝐱,𝑡)()()(),𝜕𝑡=−∇⋅𝐹𝐱,𝑡+𝑅𝐱,𝑡,inΩ×0,𝑇 (2) on a coarse mesh, where Ψ is the macroscopic state variable corresponding to 𝜓, that is, we have Ψ=𝜓 at the coarse node (𝑖,𝑗), and 𝐹 is the macroscopic flux. In fact, based on the principle of the flux balance, we can also deduce (2).

To solve (2), we firstly need to determine the macroscopic flux 𝐹. In the absence of explicit knowledge of 𝐹, our problem reduces to approximate the flux 𝐹; this will be done by locally solving a set of microscopic models. In Figure 1, we center a control volume at the midpoint of the line between any two coarse nodes except for two exterior nodes. Thus, there are four control volumes 𝐼^𝛿_{𝑖±(1/2),𝑗} and 𝐼^𝛿_{𝑖,𝑗±(1/2)} around every interior node (𝑖,𝑗) and these control volumes are centered at the points (𝑖±(1/2),𝑗) and (𝑖,𝑗±(1/2)), respectively (see Figure 2). We assume that the control volume 𝐼^𝛿_{𝑖±(1/2),𝑗} is a square of size 𝛿 and discretize it into a fine 𝑀×𝑀 mesh, for which (𝜉_𝑘,𝜂_𝑙) denotes the coordinate of node (𝑘,𝑙), where 𝑘,𝑙=1,…,𝑀+1. 𝑎=𝜉_𝑘+1−𝜉_𝑘=𝜂_𝑙+1−𝜂_𝑙 denotes the fine mesh size. For 𝜉_𝑘 and 𝜂_𝑙, we have 𝜉_𝑘=𝐻−𝛿2+()𝛿𝑘−1𝑀,𝜂_𝑙𝛿=𝐻−2+()𝛿𝑙−1𝑀,𝑘,𝑙=1,…,𝑀+1. (3) Similarly, (𝜂_𝑙,𝜉_𝑘) denotes the coordinate of the control volume 𝐼^𝛿_{𝑖,𝑗±(1/2)}.

Figure 1 

Illustration of the macroscopic and microscopic computational domains.

Figure 2 

The control volumes at the coarse node (𝑖,𝑗) and two basic microscopic elliptic models for the control volume 𝐼𝛿𝑖+(1/2),𝑗.

2.2.1. Basic Microscopic Elliptic Problems

To estimate the macroscopic flux, we need to solve a set of local microscale problems in the control volumes. Actually, the saturated hydraulic conductivity tensor 𝐊(𝐱) in this study is only a function of the spatial position and does not oscillate in the temporal direction, and we only need the spatial homogenization of 𝐊(𝐱) at the microscopic evolution. According to the conclusion of [14], in every control volume 𝐼^𝛿, we can only solve the following reduced elliptic equation: ∇⋅𝐊(𝐱)∇𝜓(𝐱)=0,in𝐼^𝛿. (4)

In every control volume 𝐼^𝛿, similar to the construction of the multiscale finite element base functions developed by Hou and Wu [2] and Hou et al. [3], we will solve two basic elliptic problems with the Dirichlet-Neumann boundary condition in which the Dirichlet boundary condition is used in one direction and no-flow boundary condition is used in the other direction. For 𝐼^𝛿_{𝑖+(1/2),𝑗} (𝑖,𝑗=1,2,…,𝑁), the Dirichlet boundary condition is used in 𝑥-direction and no-flow boundary condition is used in 𝑦-direction, and vice versa for 𝐼^𝛿_{𝑖,𝑗+(1/2)} (𝑖,𝑗=1,2,…,𝑁).

Set the head of the basic elliptic problem with no dimensional change. Let 𝜙¹ and 𝜙² be the solutions of these two basic elliptic problems, respectively. In 𝐼^𝛿_{𝑖+(1/2),𝑗}, as shown in Figure 2, for the first basic elliptic problem, the head on the left side is 1, and that on the right side is 0, and vice versa for the second basic elliptic problem. Similarly, in 𝐼^𝛿_{𝑖,𝑗+(1/2)}, for the first basic elliptic problem, the head on the bottom side is 1, and that on the top side is 0, and vice versa for the second basic elliptic problem. The cell problems are computed in parallel, and the number of the processors is reduced [16, 25]. To solve the basic elliptic problem, we considered employing the conventional finite difference method with multigrid over a fine mesh to solve the original equation. For implementing the multigrid algorithm, we use directly a numerical simulator MGD9V [26]. In fact, according to the above two basic elliptic problems, we have 𝜙¹+𝜙²=1. Then, we only need to solve the first basic elliptic problem and obtain the solution 𝜙² of the second basic elliptic problem in the course of computation according to 𝜙²=1−𝜙¹.

2.2.2. Estimation of Basic Macroscopic Fluxes

After solving the basic elliptic problems, we estimate basic macroscopic fluxes based on the solutions of the basic elliptic problems. In 𝐼^𝛿_{𝑖+(1/2),𝑗}, 𝐹^{𝑥𝛼,𝛽}_{𝑖+(1/2),𝑗}(𝛼,𝛽=1,2) denote, respectively, the basic macroscopic fluxes estimated by the solutions of two basic elliptic problems; we have 𝐹^𝑥𝛼,1_{𝑖+(1/2),𝑗}1=−𝐼^𝛿_{()𝑖+1/2,𝑗}^𝐼𝛿_{𝑖+(1/2),𝑗}𝐾(𝐱)𝜕𝜙^𝛼𝑎𝜕𝑥𝑑𝑥𝑑𝑦≈−²𝛿^2𝑀𝑀_𝑘=1𝑙=1𝐾_{𝑘+(1/2),𝑙}𝜙^𝛼_𝑘+1,𝑙−𝜙^𝛼_𝑘,𝑙𝑎,𝐹^𝑥𝛼,2_{𝑖+(1/2),𝑗}1=−𝐼^𝛿_{()𝑖+1/2,𝑗}^𝐼𝛿_{𝑖+(1/2),𝑗}𝑦𝐻𝐾(𝐱)𝜕𝜙^𝛼𝑎𝜕𝑥𝑑𝑥𝑑𝑦≈−²𝛿^2𝑀𝑀_𝑘=1𝑙=1𝜂_𝑙𝐻𝐾_{𝑘+(1/2),𝑙}𝜙^𝛼_𝑘+1,𝑙−𝜙^𝛼_𝑘,𝑙𝑎,𝛼=1,2, (5) where 𝐾_{𝑘+(1/2),𝑙} is the geometric mean of 𝐾_𝑘,𝑙 and 𝐾_𝑘+1,𝑙. Similarly, in 𝐼^𝛿_{𝑖,𝑗+(1/2)}, 𝐹^𝑦𝛼,1_{𝑖,𝑗+(1/2)}1=−𝐼^𝛿_{𝑖,𝑗+(1/2)}^𝐼𝛿_{𝑖,𝑗+(1/2)}𝐾(𝐱)𝜕𝜙^𝛼𝑎𝜕𝑦𝑑𝑥𝑑𝑦≈−²𝛿^2𝑀𝑀_𝑘=1𝑙=1𝐾_{𝑘,𝑙+(1/2)}𝜙^𝛼_𝑘,𝑙+1−𝜙^𝛼_𝑘,𝑙𝑎,𝐹^𝑦𝛼,2_{𝑖,𝑗+(1/2)}1=−𝐼^𝛿_{𝑖,𝑗+(1/2)}^𝐼𝛿_{𝑖,𝑗+(1/2)}𝑥𝐻𝐾(𝐱)𝜕𝜙^𝛼𝑎𝜕𝑦𝑑𝑥𝑑𝑦≈−²𝛿^2𝑀𝑀_𝑘=1𝑙=1𝜂_𝑘𝐻𝐾_{𝑘,𝑙+(1/2)}𝜙^𝛼_𝑘,𝑙+1−𝜙^𝛼_𝑘,𝑙𝑎,𝛼=1,2, (6) where 𝐾_{𝑘,𝑙+(1/2)} is the geometric mean of 𝐾_𝑘,𝑙 and 𝐾_𝑘,𝑙+1.

2.2.3. Estimation of Macroscopic Fluxes

Based on the estimation of the basic macroscopic fluxes, we will estimate macroscopic fluxes. Let Ψ^𝑛_𝑖,𝑗 be a coarse numerical solution of (2) at time 𝑡^𝑛, to estimate the macroscopic flux 𝐹^𝑛 at time 𝑡^𝑛; we will deal with it below.

We first solve a microscopic elliptic problem with the Dirichlet-Neumann boundary condition at every control volume 𝐼^𝛿. Head at the Dirichlet boundary of the control volume is calculated by bilinear interpolating functions.

For the control volume 𝐼^𝛿_{𝑖+(1/2),𝑗}, heads at the left and right sides are 𝜓^{𝑛,left,𝑙}_{𝑖+(1/2),𝑗}=𝑔₁𝜂_𝑙+𝜉₁𝑔₂𝜂_𝑙−𝑔₁𝜂_𝑙𝐻,𝜓^{𝑛,right,𝑙}_{𝑖+(1/2),𝑗}=𝑔₁𝜂_𝑙+𝜉_𝑀+1𝑔₂𝜂_𝑙−𝑔₁𝜂_𝑙𝐻, (7) respectively, where 𝑔₁𝜂_𝑙=𝜂_𝑙𝐻Ψ_𝑖,𝑗+𝜂1−^𝑙𝐻Ψ_{𝑖,𝑗−1}𝑀,1≤𝑙≤2,𝜂2−^𝑙𝐻Ψ_𝑖,𝑗+𝜂_𝑙𝐻−1Ψ_𝑖,𝑗+1,𝑀2𝑔<𝑙≤𝑀,₂𝜂_𝑙=𝜂_𝑙𝐻Ψ_𝑖+1,𝑗+𝜂1−^𝑙𝐻Ψ_{𝑖+1,𝑗−1}𝑀,1≤𝑙≤2,𝜂2−^𝑙𝐻Ψ_𝑖+1,𝑗+𝜂_𝑙𝐻−1Ψ_{𝑖+1,𝑗+1},𝑀2<𝑙≤𝑀. (8) Thus, 𝜓^𝑛,𝑙, the solution of the microscopic elliptic problem in 𝐼^𝛿_{𝑖+(1/2),𝑗}, is obtained over the control volume 𝐼^𝛿_{𝑖+(1/2),𝑗} via the linear combination of 𝜙¹ and 𝜙², which are the solutions of two corresponding basic elliptic problems, respectively, and then we have 𝜓^𝑛,𝑙=𝜓^{𝑛,left,𝑙}_{𝑖+(1/2),𝑗}𝜙¹+𝜓^{𝑛,right,𝑙}_{𝑖+(1/2),𝑗}𝜙²,𝑙=1,…,𝑀+1. (9) Similarly, in 𝐼^𝛿_{𝑖,𝑗+(1/2)}, heads at the bottom and top sides are 𝜓^{𝑛,bottom,𝑙}_{𝑖,𝑗+(1/2)}=𝑔₃𝜂_𝑙+𝜉₁𝑔₄𝜂_𝑙−𝑔₃𝜂_𝑙𝐻,𝜓^{𝑛,top,𝑙}_{𝑖,𝑗+(1/2)}=𝑔₃𝜂_𝑙+𝜉_𝑀+1𝑔₄𝜂_𝑙−𝑔₃𝜂_𝑙𝐻, (10) respectively, where 𝑔₃𝜂_𝑙=𝜂_𝑙𝐻Ψ_𝑖,𝑗+𝜂1−^𝑙𝐻Ψ_{𝑖−1,𝑗}𝑀,1≤𝑙≤2,𝜂2−^𝑙𝐻Ψ_𝑖,𝑗+𝜂_𝑙𝐻−1Ψ_𝑖+1,𝑗,𝑀2𝑔<𝑙≤𝑀,₄𝜂_𝑙=𝜂_𝑙𝐻Ψ_𝑖,𝑗+1+𝜂1−^𝑙𝐻Ψ_{𝑖−1,𝑗+1}𝑀,1≤𝑙≤2,𝜂2−^𝑙𝐻Ψ_𝑖,𝑗+1+𝜂_𝑙𝐻−1Ψ_{𝑖+1,𝑗+1},𝑀2<𝑙≤𝑀. (11) The solution of the microscopic elliptic problem in this control volume is 𝜓^𝑛,𝑙=𝜓^{𝑛,bottom,𝑙}_{𝑖,𝑗+(1/2)}𝜙¹+𝜓^{𝑛,top,𝑙}_{𝑖,𝑗+(1/2)}𝜙²,𝑙=1,…,𝑀+1. (12)

Like the assumption in [27], we assume that the situation investigated in this study is for locally isotropic conductivity and also assume the hydraulic conductivity tensor with principal axes oriented in the direction of the principal statistical anisotropy axes of the local parameters. It means that the conductivity tensor for a locally heterogeneous medium is 𝐾(𝐱)0(𝐱)𝐊=0𝐾. (13)

By applying the above assumption, we derive the macroscopic flux 𝐹^𝑛, 𝐹^𝑛1=−𝐼^𝛿𝐼𝛿𝐊(𝐱)∇𝜓𝐱,𝑡^𝑛1𝑑𝐱=−𝐼^𝛿𝐼𝛿𝐾(𝐱)00𝐾(𝐱)𝜕𝜓^𝑛𝜕𝑥𝜕𝜓^𝑛1𝜕𝑦𝑑𝑥𝑑𝑦=−𝐼^𝛿𝐼𝛿𝐾(𝐱)𝜕𝜓^𝑛𝐾(𝐱)𝜕𝑥𝜕𝜓^𝑛𝐹𝜕𝑦𝑑𝑥𝑑𝑦=^𝑛,𝑥𝐹^𝑛,𝑦, (14) where 𝐹^𝑛,𝑥 and 𝐹^𝑛,𝑦 are estimated macroscopic fluxes in 𝑥-direction and 𝑦-direction at time 𝑡^𝑛, respectively. For the control volume 𝐼^𝛿_{𝑖+(1/2),𝑗}, together with (3), (5), (7), (9), and (14), we have 𝐹^𝑛,𝑥_{𝑖+(1/2),𝑗}1=−𝐼^𝛿_{()𝑖+1/2,𝑗}^𝐼𝛿_{𝑖+(1/2),𝑗}𝐾(𝐱)𝜕𝜓^𝑛𝑎𝜕𝑥𝑑𝑥𝑑𝑦≈−²𝛿^2𝑀𝑀_𝑘=1𝑙=1𝐾_{𝑘+(1/2),𝑙}𝜓^𝑛,𝑙_𝑘+1,𝑙−𝜓^𝑛,𝑙_𝑘,𝑙𝑎𝑎=−²𝛿^2𝑀𝑀_𝑘=1𝑙=1𝐾_{𝑘+(1/2),𝑙}𝜓^{𝑛,left,𝑙}_{𝑖+(1/2),𝑗}𝜙¹_𝑘+1,𝑙−𝜙¹_𝑘,𝑙𝑎+𝜓^{𝑛,right,𝑙}_{𝑖+(1/2),𝑗}𝜙²_𝑘+1,𝑙−𝜙²_𝑘,𝑙𝑎≈𝐻−𝛿4𝐻2Ψ^𝑛_𝑖+1,𝑗+Ψ^𝑛_{𝑖+1,𝑗−1}−Ψ^𝑛_{𝑖+1,𝑗+1}+𝐻+𝛿4𝐻2Ψ^𝑛_𝑖,𝑗+Ψ^𝑛_{𝑖,𝑗−1}−Ψ^𝑛_𝑖,𝑗+1𝐹^𝑥1,1_{𝑖+(1/2),𝑗}+𝐻−𝛿4𝐻−Ψ^𝑛_{𝑖+1,𝑗−1}+Ψ^𝑛_{𝑖+1,𝑗+1}+𝐻+𝛿4𝐻−Ψ^𝑛_{𝑖,𝑗−1}+Ψ^𝑛_𝑖,𝑗+1𝐹^𝑥1,2_{𝑖+(1/2),𝑗}+𝐻+𝛿4𝐻2Ψ^𝑛_𝑖+1,𝑗+Ψ^𝑛_{𝑖+1,𝑗−1}−Ψ^𝑛_{𝑖+1,𝑗+1}+𝐻−𝛿4𝐻2Ψ^𝑛_𝑖,𝑗+Ψ^𝑛_{𝑖,𝑗−1}−Ψ^𝑛_𝑖,𝑗+1𝐹^𝑥2,1_{𝑖+(1/2),𝑗}+𝐻+𝛿4𝐻−Ψ^𝑛_{𝑖+1,𝑗−1}+Ψ^𝑛_{𝑖+1,𝑗+1}+𝐻−𝛿4𝐻−Ψ^𝑛_{𝑖,𝑗−1}+Ψ^𝑛_𝑖,𝑗+1𝐹^𝑥2,2_{𝑖+(1/2),𝑗}. (15) Similarly, for 𝐼^𝛿_{𝑖,𝑗+(1/2)}, we have 𝐹₍₎^𝑛,𝑦_{𝑖,𝑗+1/2}≈𝐻−𝛿4𝐻2Ψ^𝑛_𝑖,𝑗+1+Ψ^𝑛_{𝑖−1,𝑗+1}−Ψ^𝑛_{𝑖+1,𝑗+1}+𝐻+𝛿4𝐻2Ψ^𝑛_𝑖,𝑗+Ψ^𝑛_{𝑖−1,𝑗}−Ψ^𝑛_𝑖+1,𝑗𝐹^𝑦1₍₎^,1_{𝑖,𝑗+1/2}+𝐻−𝛿4𝐻−Ψ^𝑛_{𝑖−1,𝑗+1}+Ψ^𝑛_{𝑖+1,𝑗+1}+𝐻+𝛿4𝐻−Ψ^𝑛_{𝑖−1,𝑗}+Ψ^𝑛_𝑖+1,𝑗𝐹^𝑦1₍₎^,2_{𝑖,𝑗+1/2}+𝐻+𝛿4𝐻2Ψ^𝑛_𝑖,𝑗+1+Ψ^𝑛_{𝑖−1,𝑗+1}−Ψ^𝑛_{𝑖+1,𝑗+1}+𝐻−𝛿4𝐻2Ψ^𝑛_𝑖,𝑗+Ψ^𝑛_{𝑖−1,𝑗}−Ψ^𝑛_𝑖+1,𝑗𝐹^𝑦2₍₎^,1_{𝑖,𝑗+1/2}+𝐻+𝛿4𝐻−Ψ^𝑛_{𝑖−1,𝑗+1}+Ψ^𝑛_{𝑖+1,𝑗+1}+𝐻−𝛿4𝐻−Ψ^𝑛_{𝑖−1,𝑗}+Ψ^𝑛_𝑖+1,𝑗𝐹^𝑦2₍₎^,2_{𝑖,𝑗+1/2}. (16)

2.2.4. Macroscopic Evolution

Let Δ𝑡 be a time step size, where Δ𝑡=𝑡^𝑛+1−𝑡^𝑛. The macroscopic evolution on the coarse mesh is now done via the approximation of (2), and here we use the Crank-Nicolson method to construct the fully discretized version of (2), 𝑆_{𝑠𝑖,𝑗}ΨΔ𝑡^{𝑛+1𝑖,𝑗}−Ψ^𝑛_𝑖,𝑗1=−𝐹2𝐻^{𝑛+1,𝑥𝑖+(1/2),𝑗}−𝐹^𝑛+1,𝑥_{𝑖−(1/2),𝑗}+𝐹^𝑛+1,𝑦_{𝑖,𝑗+(1/2)}−𝐹^𝑛+1,𝑦_{𝑖,𝑗−(1/2)}−1𝐹2𝐻^{𝑛,𝑥𝑖+(1/2),𝑗}−𝐹^𝑛,𝑥_{𝑖−(1/2),𝑗}+𝐹^𝑛,𝑦_{𝑖,𝑗+(1/2)}−𝐹^𝑛,𝑦_{𝑖,𝑗−(1/2)}+𝑅^𝑛_𝑖,𝑗. (17) Combining (15), (16), and (17) yields where 𝑐₁1=−8𝐻²×(𝐻+𝛿)−𝐹^𝑥1,1_{𝑖−(1/2),𝑗}+𝐹^𝑥1,2_{𝑖−(1/2),𝑗}−𝐹^𝑦1,1_{𝑖,𝑗−(1/2)}+𝐹^𝑦1,2_{𝑖,𝑗−(1/2)}+(𝐻−𝛿)−𝐹^𝑥2,1_{𝑖−(1/2),𝑗}+𝐹^𝑥2,2_{𝑖−(1/2),𝑗}−𝐹^𝑦2,1_{𝑖,𝑗−(1/2)}+𝐹^𝑦2,2_{𝑖,𝑗−(1/2)},𝑐₂1=−8𝐻²(𝐻+𝛿)𝐹^𝑥1,1_{𝑖+(1/2),𝑗}−𝐹^𝑥1,2_{𝑖+(1/2),𝑗}−𝐹^𝑥2,1_{𝑖−(1/2),𝑗}+𝐹^𝑥2,2_{𝑖−(1/2),𝑗}−2𝐹^𝑦1,1_{𝑖,𝑗−(1/2)}+(𝐻−𝛿)𝐹^𝑥2₍₎^,1_{𝑖+1/2,𝑗}−𝐹^𝑥2₍₎^,2_{𝑖+1/2,𝑗}−𝐹^𝑥1₍₎^,1_{𝑖−1/2,𝑗}+𝐹^𝑥1₍₎^,2_{𝑖−1/2,𝑗}−2𝐹^𝑦2₍₎^,1_{𝑖,𝑗−1/2},𝑐₃1=−8𝐻²(𝐻+𝛿)𝐹^𝑥2₍₎^,1_{𝑖+1/2,𝑗}−𝐹^𝑥2₍₎^,2_{𝑖+1/2,𝑗}+𝐹^𝑦1₍₎^,1_{𝑖,𝑗−1/2}−𝐹^𝑦1₍₎^,2_{𝑖,𝑗−1/2}+(𝐻−𝛿)𝐹^𝑥1₍₎^,1_{𝑖+1/2,𝑗}−𝐹^𝑥1₍₎^,2_{𝑖+1/2,𝑗}+𝐹^𝑦2₍₎^,1_{𝑖,𝑗−1/2}−𝐹^𝑦2₍₎^,2_{𝑖,𝑗−1/2},𝑐₄1=−8𝐻²(𝐻+𝛿)𝐹^𝑦1₍₎^,1_{𝑖,𝑗+1/2}−𝐹^𝑦1₍₎^,2_{𝑖,𝑗+1/2}−2𝐹^𝑥1₍₎^,1_{𝑖−1/2,𝑗}−𝐹^𝑦1₍₎^,1_{𝑖,𝑗−1/2}+𝐹^𝑦2₍₎^,2_{𝑖,𝑗−1/2}+(𝐻−𝛿)𝐹^𝑦2₍₎^,1_{𝑖,𝑗+1/2}−𝐹^𝑦2₍₎^,2_{𝑖,𝑗+1/2}−2𝐹^𝑥2₍₎^,1_{𝑖−1/2,𝑗}+𝐹^𝑦1₍₎^,2_{𝑖,𝑗−1/2},𝑐₅1=−8𝐻²(𝐻+𝛿)2𝐹^𝑥1₍₎^,1_{𝑖+1/2,𝑗}+2𝐹^𝑦1₍₎^,1_{𝑖,𝑗+1/2}−2𝐹^𝑥2₍₎^,1_{𝑖−1/2,𝑗}−2𝐹^𝑦2₍₎^,1_{𝑖,𝑗−1/2}+(𝐻−𝛿)2𝐹^𝑥2₍₎^,1_{𝑖+1/2,𝑗}+2𝐹^𝑦2₍₎^,1_{𝑖,𝑗+1/2}−2𝐹^𝑥1₍₎^,1_{𝑖−1/2,𝑗}−2𝐹^𝑦1₍₎^,1_{𝑖,𝑗−1/2},𝑐₆1=−8𝐻²(𝐻+𝛿)2𝐹^𝑥2₍₎^,1_{𝑖+1/2,𝑗}−𝐹^𝑦1₍₎^,1_{𝑖,𝑗+1/2}+𝐹^𝑦1₍₎^,2_{𝑖,𝑗+1/2}+𝐹^𝑦2₍₎^,1_{𝑖,𝑗−1/2}−𝐹^𝑦2₍₎^,2_{𝑖,𝑗−1/2}+(𝐻−𝛿)2𝐹^𝑥1₍₎^,1_{𝑖+1/2,𝑗}−𝐹^𝑦2₍₎^,1_{𝑖,𝑗+1/2}+𝐹^𝑦2₍₎^,2_{𝑖,𝑗+1/2}−𝐹^𝑦1₍₎^,2_{𝑖,𝑗−1/2},𝑐₇1=−8𝐻²(𝐻+𝛿)𝐹^𝑦2₍₎^,1_{𝑖,𝑗+1/2}−𝐹^𝑦2₍₎^,2_{𝑖,𝑗+1/2}+𝐹^𝑥1₍₎^,1_{𝑖−1/2,𝑗}−𝐹^𝑥1₍₎^,2_{𝑖−1/2,𝑗}+(𝐻−𝛿)𝐹^𝑦1₍₎^,1_{𝑖,𝑗+1/2}−𝐹^𝑦1₍₎^,2_{𝑖,𝑗+1/2}+𝐹^𝑥2₍₎^,1_{𝑖−1/2,𝑗}−𝐹^𝑥2₍₎^,2_{𝑖−1/2,𝑗},𝑐₈1=−8𝐻²(𝐻+𝛿)−𝐹^𝑥1₍₎^,1_{𝑖+1/2,𝑗}+𝐹^𝑥1₍₎^,2_{𝑖+1/2,𝑗}+2𝐹^𝑦2₍₎^,1_{𝑖,𝑗+1/2}+𝐹^𝑥2₍₎^,1_{𝑖−1/2,𝑗}−𝐹^𝑥2₍₎^,2_{𝑖−1/2,𝑗}+(𝐻−𝛿)−𝐹^𝑥2₍₎^,1_{𝑖+1/2,𝑗}+𝐹^𝑥2₍₎^,2_{𝑖+1/2,𝑗}+2𝐹^𝑦1,1_{𝑖,𝑗+(1/2)}+𝐹^𝑥1₍₎^,1_{𝑖+1/2,𝑗}−𝐹^𝑥1₍₎^,2_{𝑖−1/2,𝑗},𝑐₉1=−8𝐻²(𝐻+𝛿)−𝐹^𝑥2₍₎^,1_{𝑖+1/2,𝑗}+𝐹^𝑥2₍₎^,2_{𝑖+1/2,𝑗}+𝐹^𝑦2₍₎^,1_{𝑖,𝑗+1/2}+𝐹^𝑦2₍₎^,2_{𝑖,𝑗+1/2}+(𝐻−𝛿)−𝐹^𝑥1₍₎^,1_{𝑖+1/2,𝑗}+𝐹^𝑥1₍₎^,2_{𝑖+1/2,𝑗}−𝐹^𝑦1₍₎^,1_{𝑖,𝑗+1/2}+𝐹^𝑦1₍₎^,2_{𝑖,𝑗+1/2}. (19)

Then, we solve the following equation at the coarse mesh by using MGD9V [26]: −𝑐₁Ψ^𝑛+1_{𝑖−1,𝑗−1}−𝑐₂Ψ^𝑛+1_{𝑖,𝑗−1}−𝑐₃Ψ^𝑛+1_{𝑖+1,𝑗−1}−𝑐₄Ψ^𝑛+1_{𝑖−1,𝑗}+𝑆_{𝑠𝑖,𝑗}Δ𝑡−𝑐₅Ψ^𝑛+1_𝑖,𝑗−𝑐₆Ψ^𝑛+1_𝑖+1,𝑗−𝑐₇Ψ^𝑛+1_{𝑖−1,𝑗+1}−𝑐₈Ψ^𝑛+1_𝑖,𝑗+1−𝑐₉Ψ^𝑛+1_{𝑖+1,𝑗+1}=𝑐₁Ψ^𝑛_{𝑖−1,𝑗−1}+𝑐₂Ψ^𝑛_{𝑖,𝑗−1}+𝑐₃Ψ^𝑛_{𝑖+1,𝑗−1}+𝑐₄Ψ^𝑛_{𝑖−1,𝑗}+𝑆_{𝑠𝑖,𝑗}Δ𝑡+𝑐₅Ψ^𝑛_𝑖,𝑗+𝑐₆Ψ^𝑛_𝑖+1,𝑗+𝑐₇Ψ^𝑛_{𝑖−1,𝑗+1}+𝑐₈Ψ^𝑛_𝑖,𝑗+1+𝑐₉Ψ^𝑛_{𝑖+1,𝑗+1}+𝑅^𝑛_𝑖,𝑗. (20)

Thus, the algorithm is completed. The solution procedure at a time step is illustrated in Figure 3. We also give a summary which only includes the relevant discrete equations necessary to implement the proposed method. Firstly, the basic elliptic problems (4) are solved. Then, basic macroscopic fluxes (5) and (6) are estimated, and the coefficients (19) are calculated. The above steps are finished at the preprocessing step. Finally, the macroscopic discrete equations (20) are evolved.

Figure 3 

Flow chart of the new scheme at a time step.

The algorithm described above is easy to extend to the steady flow problem in heterogeneous porous media. Under the condition of the steady flow, the left-hand side of (1) equals zero. Correspondingly, the left-hand side of macroscopic (2) equals zero. The remainder will similarly be completed; then we have 𝑐₁Ψ_{𝑖−1,𝑗−1}+𝑐₂Ψ_{𝑖,𝑗−1}+𝑐₃Ψ_{𝑖+1,𝑗−1}+𝑐₄Ψ_{𝑖−1,𝑗}+𝑐₅Ψ_𝑖,𝑗+𝑐₆Ψ_𝑖+1,𝑗+𝑐₇Ψ_{𝑖−1,𝑗+1}+𝑐₈Ψ_𝑖,𝑗+1+𝑐₉Ψ_{𝑖+1,𝑗+1}1=−2𝑅_𝑖,𝑗. (21) Although FDHMM proposed by Abdulle and E [20] works for the transient problem, it is also easy to extend to the steady problem. The given coarse- and fine-grid equations of FDHMM are elliptic equations in the steady condition. The iteration scheme of FDHMM is similar to that of the new scheme, and the solutions of both methods are the same.

The locally hydraulic conductivity (13) is assumed to be isotropic, but the macroscopic conductivity may be anisotropic or a full tensor because the final macroscopic scheme (18) is a nine-spot one. The algorithm may also be extended to general quadrilateral mesh; the method is similar to [28].

3. Evaluation of Numerical Accuracy

All porous media in nature are heterogeneous. The heterogeneity in this study comes from the hydraulic conductivity. As the standard deviation of logarithmic hydraulic conductivity increases, the heterogeneity increases. The random conductivity field is generated by the Turning Band method [29], in which the hydraulic conductivity is assumed to be locally isotropic. In this study, we used four saturated groundwater flow examples, including two steady flow examples and two transient flow examples, to show the main advantages of the constructed scheme. Also, influences of different factors are examined, such as conductivity fields with high variability as well as different correlation structures, the flow rate of the pumping well, and the size of the local microscale model, on the accuracy of the constructed scheme.

3.1. Implementation

The algorithm has been implemented in a FORTRAN code. Because it is difficult to construct interesting multiscale problem with an exact solution, people often compare the coarse scale solution obtained by the multiscale method with a computed reference solution obtained on the fine scale. We have employed the conventional finite difference method with multigrid over a fine mesh to solve the original equation and refer to this solution as the “exact” solution.

As a measure of the error, we take the relative 𝐿₂ norm and the relative maximum norm eer₂=∑^𝑁′_𝑖=1Ψ_𝑖−Ψ𝐱_𝑖²∑^𝑁′_𝑖=1Ψ𝐱_𝑖^21/2,eer_∞=max_{𝑖=1,…,𝑁′}Ψ_𝑖−Ψ𝐱_𝑖max_{𝑖=1,…,𝑁′}Ψ𝐱_𝑖, (22) respectively, where 𝑁^′ is the total number of nodes on the coarse mesh, Ψ_𝑖 denotes the coarse solution at the coarse node 𝐱_𝑖, and Ψ(𝐱_𝑖) denotes the “exact” solution projected on the coarse mesh; that is, Ψ(𝐱_𝑖)=𝜓(𝐱_𝑖) at the coarse node. Here, 𝜓(𝐱) is the “exact” solution.

In all test examples, the study domain Ω is a rectangle covering 1 km × 1 km with the point (0,0) as the origin. A uniform finite difference mesh is constructed by dividing Ω into an 𝑁×𝑁 mesh. The fine mesh is a 256×256 mesh, and the “exact” solution and the random hydraulic conductivity field are obtained on this mesh. The coarse mesh is a 16×16 mesh and the coarse solution is obtained by using the multiscale method on this mesh.

3.2. Steady Flow Problems with Isotropic and Anisotropic Microstructure

We impose the Dirichlet-Neumann boundary condition for the test steady flow problem. The left and right sides of boundary are Dirichlet boundaries. Head on the left is 20 m, and that on the right side is 10 m. The top and bottom sides are impermeable boundaries. To start the computation using the new scheme, we need to choose the size of the control volume 𝛿. In this study, 𝛿 is chosen to be equal to a half of the coarse mesh size, which means that we only use about 50% of the total data at the small scale; that is, 𝛿=(1/2)𝐻. Every control volume 𝐼^𝛿 is uniformly divided into an 8×8 mesh such that its mesh size equals the size of the fine mesh.

Four conductivity fields with isotropic correlation microstructure are first applied. Correlation lengths of these conductivity fields are 𝜆_𝑥=𝜆_𝑦=10 m, 𝜆_𝑥=𝜆_𝑦=20 m, 𝜆_𝑥=𝜆_𝑦=40 m, and 𝜆_𝑥=𝜆_𝑦=100 m, respectively. We assume that the geometric mean of hydraulic conductivity is 𝐾=0.006 m/min. Under 𝜎_ln𝐾=0.5,1.0,1.5,2.0, where 𝜎_ln𝐾 is the standard deviation of logarithmic hydraulic conductivity, these four conductivity fields vary by about one, three, five, and six orders of magnitude, respectively. A realization of the random conductivity fields with 𝜆_𝑥=𝜆_𝑦=100 m and 𝜎_ln𝐾=0.5,1.0,1.5,2.0 is plotted in Figure 4. Figure 5 plots the errors of the solutions of the constructed scheme for different correlation lengths at 𝜎_ln𝐾=0.5,1.0,1.5,2.0. It illustrates that the larger standard deviation of logarithmic hydraulic conductivity leads to the less accurate results. Furthermore, at the case with 𝜆_𝑥=𝜆_𝑦=10 m and 𝜎_ln𝐾=2.5, in which the conductivity field varies by about nine orders of magnitude, the new scheme is even not convergent, although it works for other correlation scales for the same logarithmic conductivity variance. The reason may be that conductivities of highly heterogeneous systems are highly discontinuous, which makes the direct application of the algorithm infeasible. It may also indicate that the standard deviation of logarithmic hydraulic conductivity plays an important role in the accuracy of the new scheme. At the same time, when 𝜎_ln𝐾=0.5, the results obtained under different correlation lengths have about the same accuracy and the correlation length of conductivity field shows no significant effect on the accuracy of the new scheme. It is noted that the correlation length is even larger than the size of the control volume in the case with 𝜆_𝑥=𝜆_𝑦=100 m. The heads in section 𝑦=500 m obtained from the fine-scale model on the fine mesh and the constructed scheme on the coarse mesh for the case with 𝜆_𝑥=𝜆_𝑦=100 m and 𝜎_ln𝐾=0.5,1.0,1.5,2.0 are depicted in Figure 6, and we observe that the solution obtained by the constructed scheme on a coarse mesh is able to approximate the exact solution. The above discussion indicates that the new scheme gives a reasonable accuracy for the test steady flow examples with isotropic correlation microstructure.

(a)

(b)

(c)

(d)

Figure 4 

A realization of the random saturated conductivity fields of different standard deviations under 𝜆𝑥=𝜆𝑦=100 m: (a) 𝜎ln𝐾=0.5, (b) 𝜎ln𝐾=1.0, (c) 𝜎ln𝐾=1.5, and (d) 𝜎ln𝐾=2.0.

(a)

(b)

Figure 5 

Relative (a) 𝐿2 and (b) maximum errors between the fine-scale model and the new scheme for the steady flow problem with isotropic correlation microstructure.

(a)

(b)

(c)

(d)

Figure 6 

Exact solution and the coarse solution of the new scheme in section 𝑦=500 m for the steady flow problem under isotropic correlation microstructure with different standard deviations: (a) 𝜎ln𝐾=0.5, (b) 𝜎ln𝐾=1.0, (c) 𝜎ln𝐾=1.5, and (d) 𝜎ln𝐾=2.0.

Convergence should be a necessary condition for the new scheme as a good numerical method. Here, we only consider the conductivity field with 𝜆_𝑥=𝜆_𝑦=100 m. Fixing 𝛿=(1/2)𝐻, Figure 7 plotted the relative errors for coarse meshes with 4×4, 8×8, 16×16, and 32×32 elements under four cases that 𝜎_ln𝐾=0.5,1.0,1.5,2.0, respectively. The errors monotonically decrease as the total number of coarse elements increases and tend to zero, which means that the new scheme solution converges as the coarse grid is refined.

(a)

(b)

Figure 7 

Relative (a) 𝐿2 and (b) maximum errors of the solution of the new scheme for coarse grids with 4×4, 8×8, 16×16, and 32×32 elements.

Next, we turn to consider three conductivity fields with anisotropic correlation microstructure. Fixing 𝑦-direction correlation length 𝜆_𝑦=10 m, 𝑥-direction correlation lengths of these conductivity fields are 𝜆_𝑥=20 m, 𝜆_𝑥=40 m, and 𝜆_𝑥=100 m, respectively. Assume that 𝐾=0.006 m/min. Under 𝜎_ln𝐾=0.5,1.0,1.5,2.0, these three conductivity fields vary by over one, three, five, and seven orders of magnitude, respectively. The results of the constructed scheme for different correlation lengths at 𝜎_ln𝐾=0.5,1.0,1.5,2.0 are depicted in Figure 8. As in the previous example, the standard deviation of logarithmic hydraulic conductivity shows significant effect on the accuracy of the new scheme. The errors obviously increase with 𝜎_ln𝐾 increasing. Compared to the standard deviation, the correlation length of conductivity field has relatively little influence on the accuracy of the new scheme. The maximum error in Figure 8 is attained at the case with 𝜆_𝑥=100, 𝜆_𝑦=10 m, and 𝜎_ln𝐾=2.0, and the relative 𝐿₂ and maximum errors of the solution of the constructed scheme are 2.28% and 5.88%, respectively. Similar to the isotropic case, the new scheme also gives a reasonable accuracy for the test steady flow examples with anisotropic correlation microstructure.

(a)

(b)

Figure 8 

Relative (a) 𝐿2 and (b) maximum errors between the fine-scale model and the new scheme for the steady flow problem with anisotropic correlation microstructure.

3.3. Aquifer Response to Sudden Change in Reservoir Level

We design this transient test example based on the example in [30]. Consider the confined aquifer in the study area. Initial head is equal to 20 m everywhere in the aquifer. We wish to simulate changes in head through time if, at 𝑡=0, we suddenly drop the water level in the reservoir on four sides of the study area from 20 m to 10 m. The specific storage coefficient and the thickness of the aquifer are 2.0×10⁻⁴ m⁻¹ and 10 m, respectively. To generate the random hydraulic conductivity field, we assume that 𝐾 is 0.003 m/min, the standard deviation of ln𝐾 is 1.5, and the correlation structure of the conductivity is anisotropic with 𝜆_𝑥=40 m and 𝜆_𝑦=10 m. Conductivity 𝐾 in this random field varies by over five orders of magnitude. Fix a time step size Δ𝑡=5 min and the study time 8000 min.

At first, accuracies and efficiencies of the constructed scheme and FDHMM are compared. Let the size of the control volume 𝛿=(1/2)𝐻, and the control volume 𝐼^𝛿 is uniformly divided into an 8×8 mesh. The computational results of different multiscale schemes at times =500,1000,2000,4000,5000,6000, and 8000 min are plotted in Figure 9. Figure 9 indicates that the new scheme seems to be more accurate than FDHMM. After 𝑡=1000 min, eer₂ and eer_∞ of the solution of the constructed scheme monotonically decrease from 2.06% to 0.32% and from 7.37% to 0.95%, respectively, while those of FDHMM fluctuate in the intervals 2.71%~6.60% and 6.13%~15.40%, respectively. The reason leading to the difference between the results obtained by the constructed scheme and obtained by FDHMM may be the different approaches of estimating the macroscopic flux. Compared with FDHMM in which the approximation of the macroscopic flux is determined before the coarse equation is solved, for the constructed scheme, the computation of the macroscopic flux is coupled into the course of solving the coarse equation. Thus, a quasistationary state of the computed macroscopic flux is approached in the global domain for the constructed scheme and in the local domain for FDHMM, and it may make the constructed scheme give better accuracy than FDHMM for this test problem. We plot the heads in the whole study domain at times 𝑡=1000 and 5000 min obtained from the fine-scale model on the fine mesh and two multiscale methods on the coarse mesh (Figure 10). We observe that the heads obtained by the new scheme on a coarse mesh can satisfyingly approximate the “exact” heads, and FDHMM underestimates the heads at the coarse nodes.

(a)

(b)

Figure 9 

Relative (a) 𝐿2 and (b) maximum errors between exact and two coarse solutions for aquifer response to sudden change in reservoir level.

(a)

(b)

(c)

Figure 10 

Exact solution (top), the coarse solutions of the new scheme (middle), and FDHMM (bottom) at times 𝑡=1000 min (left) and 5000 min (right) for aquifer response to sudden change in reservoir level.

The results were obtained on a computer running Windows XP with 2.66 GHz processor, 2 megabytes of cache, and 512 megabytes of RAM. For this test example, memory requirements using the conventional finite difference method, the constructed scheme, and FDHMM are about 27.7, 4.3, and 4.3 megabytes, respectively; CPU times using the three methods are about 12.1 min, 0.1 min, and 7.2 min, respectively. Compared with the computational cost of the conventional finite difference method, in our test example, the present new can save about 84.5% memory and about 99.2% CPU time, and FDHMM can save about 84.5% memory and about 40.5% CPU time. We need to solve 4𝑁(𝑁−1) basic microscopic elliptic problems in the constructed scheme. In fact, the computations of these basic microscale problems can be carried out sequentially, and, at a time, we only need to solve two basic microscale problems (one for 𝐼^𝛿_{𝑖+(1/2),𝑗} and the other for 𝐼^𝛿_{𝑖,𝑗+(1/2)}). When 𝛿=(1/2)𝐻, every basic microscale problem has (𝑀+1)² unknowns; the degrees of freedom of these basic microscale problems are 2(𝑀+1)². Added (𝑁+1)² degrees of freedom of the macroscopic scheme, the total degrees of freedom of the new scheme are 2(𝑀+1)²+(𝑁+1)². Similarly, the total degrees of freedom of FDHMM are also 2(𝑀+1)²+(𝑁+1)², and degrees of freedom of the full fine scheme are (2𝑁𝑀+1)². Thus, both the constructed scheme and FDHMM can obviously save the memory requirement. The first saving in computational time in two multiscale schemes is achieved by reducing the computation of the fine mesh on the whole domain. The fine-scale global flow solution is decomposed into a series of local microscopic problems, and this greatly saves the computational time. However, local microscopic models of the new scheme only need to be solved once at the preprocessing step, while those of FDHMM need to be solved at every time step. Thus, the constructed scheme needs much less CPU time than FDHMM.

Next, we discuss the effects of different cell sizes on the accuracy of the constructed scheme. In a coarse 16×16 mesh, the coarse mesh size 𝐻 equals 1000/16 m; we change the size of the control volume and let 𝛿=(1/2)𝐻,(3/4)𝐻,𝐻,(5/4)𝐻 in turn. To obtain the same size as the full fine mesh size, these control volumes are uniformly divided into 8×8, 12×12, 16×16, and 20×20 meshes in turn. The results obtained by the constructed scheme under different control volume sizes at times =500,1000,2000,4000,5000,6000, and 8000 min are depicted in Figure 11. We observe that the cases with 𝛿=(3/4)𝐻, 𝛿=𝐻, and 𝛿=(5/4)𝐻 have about the same accuracy, and the case with 𝛿=(1/2)𝐻 has a less accuracy. This is likely because, at three cases with 𝛿=(3/4)𝐻, 𝛿=𝐻, and 𝛿=(5/4)𝐻, the main microstructural information is efficiently captured by the control volume. It may indicate that the control volume size shows no significant effect on the accuracy of the constructed scheme when it is chosen to be near the coarse mesh size.

(a)

(b)

Figure 11 

Relative (a) 𝐿2 and (b) maximum errors between the fine-scale model and the new scheme under different cell sizes for aquifer response to sudden change in reservoir level.

3.4. Steady and Transient Flow Problems with Weak Well Drawdown

In this section, we first consider the steady flow problem with well drawdown in heterogeneous porous media. Similar to the examples discussed in [10, 31], we impose the following fixed head and no flux boundary conditions for the test example: heads on the left and right sides are 10 m and top and bottom sides are impermeable boundaries. In addition, a pumping well with the constant flow rate 𝑄 is located at the point (500 m, 500 m), and we let 𝑄=0.12 m³/min, 0.24 m³/min, 0.36 m³/min, and 0.48 m³/min, respectively. The aquifer is 10 m thick. We also choose 𝛿=(1/2)𝐻 and uniformly divide every control volume 𝐼^𝛿 into an 8×8 mesh such that its mesh size equals the size of the fine mesh.

Four conductivity fields with 𝜎_ln𝐾=0.5,1.0,1.5,2.0 are considered. Assume that the geometric mean of hydraulic conductivity is 𝐾=0.018 m/min and the anisotropic correlation microstructure with 𝜆_𝑥=40, 𝜆_𝑦=10 m. The errors of the results obtained by the constructed scheme under different well flow rates at 𝜎_ln𝐾=0.5,1.0,1.5,2.0 are plotted in Figure 12. Different from the examples of Section 3.2, the standard deviation of logarithmic hydraulic conductivity shows no significant effect on the accuracy of the new scheme. For example, given 𝑄=0.12 m³/min, when 𝜎_ln𝐾=0.5,1.0,1.5,2.0, relative 𝐿₂ errors of the solution of the new scheme are about 0.15%, 0.15%, 0.17%, and 0.20%, respectively, and relative maximum errors of the solution of the new scheme are 2.52%, 2.29%, 2.17%, and 2.00%, respectively. The important factor affecting the accuracy of the new scheme is the flow rate of the pumping well. The larger the flow rate of pumping well is, the larger the resulting errors are. Setting 𝜎_ln𝐾=1.0, Figure 13 plots the heads in section 𝑦=500 m obtained from the fine-scale model on the fine mesh and the constructed scheme on the coarse mesh for the cases 𝑄=0.12,0.24,0.36,0.48 m³/min. There are larger errors of the results of the constructed scheme near point (500 m, 500 m), which are caused by the pumping well at this point. This is likely because heads near the well vary nonlinearly with distance to the well, which cannot be well described by the constructed scheme. On the other hand, the problem of the well singularity may be related to the chosen scale. If we choose a coarse 32×32 mesh and 𝛿=(1/2)𝐻 and resolve this well drawdown problem. When 𝜎_ln𝐾=1.0, in Figure 14, we replot the curves shown in Figure 13. We observe that the accuracy of the new scheme was improved markedly.

(a)

(b)

Figure 12 

Relative (a) 𝐿2 and (b) maximum errors between the fine-scale model and the new scheme for the steady flow problem with well drawdown.

(a)

(b)

(c)

(d)

Figure 13 

Exact solution and the coarse solution of the new scheme in section 𝑦=500 m for the steady well drawdown problem under different well flow rates: (a) 𝑄=0.12, (b) 𝑄=0.24, (c) 𝑄=0.36, and (d) 𝑄=0.48  (m3(min)−1).

(a)

(b)

(c)

(d)

Figure 14 

Exact solution and the coarse solution of the new scheme at the coarse 32×32 mesh in section 𝑦=500 m under different well flow rates: (a) 𝑄=0.12, (b) 𝑄=0.24, (c) 𝑄=0.36, and (d) 𝑄=0.48  (m3(min)−1).

Next, we consider the transient well drawdown problem in heterogeneous porous media. Boundaries of the study area are Dirichlet types. Heads on four sides are all 10 m. Initial pressure head is also 10 m everywhere in the aquifer. The specific storage coefficient is 2.0×10⁻⁴ m⁻¹ and the aquifer is 10 m thick. There is a pumping well at the point (500 m, 500 m). The well has the constant flow rate of 0.24 m³/min and is pumped for 1600 min in the problem. The time step is 1 min for every method. This test example is analogous to the examples used in [10, 31]. The statistical parameters used to describe random conductivity field are 𝐾=0.018 m/min, 𝜎_ln𝐾=1.0, 𝜆_𝑥=40 m, and 𝜆_𝑦=10 m. This random conductivity field varies by over three orders of magnitude.

Similar to Section 3.3, we compare accuracies and efficiencies of the constructed scheme and FDHMM. The control volume 𝐼^𝛿 has a size of (1/2)𝐻 and is uniformly divided into an 8×8 mesh. We plot the computational results of different multiscale schemes at times =100,200,400,800,1000,1200, and 1600 min (Figure 15). The constructed scheme gives a little more accuracy than FDHMM. Over the whole simulating time, eer₂ and eer_∞ of the solution of the new scheme are less than 0.26% and 4.15%, respectively, and those of FDHMM are less than 0.38% and 5.44%, respectively. Figure 16 shows the heads at times 𝑡=200 and 1000 min in section 𝑦=500 m obtained from the fine-scale model on the fine mesh and two multiscale methods on the coarse mesh. We observe that the solutions obtained by both the constructed scheme and FDHMM on a coarse mesh are able to satisfyingly approximate the exact solution except of the well singularity. Near the well singularity, the results obtained by both multiscale methods are in rough agreement with those obtained by the fine-scale model, and the heads are overestimated to be about 0.45 m and 0.50 m by the constructed scheme and FDHMM at the well singularity, respectively, at 𝑡=200 min. This fact may imply that, although a quasibalance state of the macroscopic flux is achieved in the global domain for the constructed scheme versus in the local domain for FDHMM, this advantage of the constructed scheme is not obvious for the well drawdown problem.

(a)

(b)

Figure 15 

Relative (a) 𝐿2 and (b) maximum errors between exact and two coarse solutions for the transient flow problem with weak well drawdown.

(a)

(b)

Figure 16 

Exact solution and two coarse solutions at times (a) 𝑡=200 min and (b) 1000 min in section 𝑦=500 m for the transient flow problem with weak well drawdown.

Computational costs of the three methods in this test example are similar to those in the test example discussed in Section 3.3 and are omitted in this section. This is followed by a discussion of effects of different cell sizes on the accuracy of the constructed scheme. As described in Section 3.3, in a coarse 16×16 mesh, the control volumes are chosen to have sizes of 𝛿=(1/2)𝐻,(3/4)𝐻,𝐻,(5/4)𝐻 and are uniformly divided into 8×8, 12×12, 16×16, and 20×20 meshes in turn. Plotted in Figure 17 are the calculated results of the constructed scheme under different control volume sizes at times 𝑡=100,200,400,800,1000,1200, and 1600 min. It indicates that four cases give a reasonable accuracy in eer₂ and eer_∞. We observe that the results for 𝛿=𝐻 are the best, the results for 𝛿=(3/4)𝐻 are less accurate than those for 𝛿=𝐻, the results for 𝛿=(5/4)𝐻 are less accurate than those for 𝛿=(3/4)𝐻, and the results for 𝛿=(1/2)𝐻 are the worst. The results obtained under three cases with 𝛿=(3/4)𝐻,𝐻,(5/4)𝐻 are very similar. Thus, the control volume size may be chosen to be near the size of the coarse mesh for the sake of the accuracy of the constructed scheme. In addition, under the cases with 𝛿=(1/2)𝐻 and 𝛿=(3/4)𝐻, we only use about 50% and 75% of the information of the whole microstructure, respectively. This flexibility of choosing the size of the control volume means that the constructed scheme may be applied to the flow problem for which the microstructure cannot be completely found beforehand.

(a)

(b)

Figure 17 

Relative (a) 𝐿2 and (b) maximum errors between the fine-scale model and the new scheme under different cell sizes for the transient flow problem with weak well drawdown.

4. Conclusion

A new scheme of the finite difference heterogeneous multiscale method, which puts more emphasis on the interaction between the macro- and microscale behaviors, has been presented for solving saturated water flow problems in random porous media. The macroscopic iteration formulas of steady and transient flow problems have been explicitly deduced. By solving basic microscopic elliptic problems and estimating basic macroscopic fluxes, it is subtly brought to the large scale for microscale information of the medium property and useful information about the gradients of the solutions of basic microscopic elliptic models. For the transient flow problem, different from that FDHMM needs the macroscopic and microscopic evolution at every time step, the constructed scheme implements the microscopic evolution at the preprocessing step and only needs the macroscopic evolution at every time step, which offers substantial saving in the computational cost. The constructed scheme saves about 58.7% CPU time compared to FDHMM for aquifer response to sudden change in reservoir level on the case 𝛿=(1/2)𝐻. Different numerical examples, including two steady flow problems and two transient flow problems subject to Dirichlet-Neumann boundary type, are applied to test the efficiency and accuracy of the constructed scheme. We have considered seven correlation lengths and four standard deviations of the hydraulic conductivity field for steady flow problems with isotropic and anisotropic microstructure and considered four flow rates of the pumping well and four standard deviations of the hydraulic conductivity field for the steady flow problem with well drawdown. In every transient flow problem, we have also considered four sizes of the control volume. The numerical experiments demonstrate that the constructed scheme gives a better accuracy than FDHMM for aquifer response to sudden change in reservoir level and gives a comparable accuracy to FDHMM for the weak well drawdown problem. The numerical experiments also indicate that the constructed scheme can efficiently capture the macroscale behavior of the solution on a coarse mesh for the steady and transient flow problems without well drawdown, and the scheme can approximately handle the weak well drawdown problem. The well singularity is related to the chosen scale. We may refine the coarse mesh size to improve the accuracy of the solution to the well drawdown problem. The standard deviation of logarithmic hydraulic conductivity field plays an important role in the accuracy of the constructed scheme. The larger the standard deviation is, the less accurate the results are. The spatial correlation length of random conductivity field has relatively little influence on the accuracy of the constructed scheme. To obtain a reasonable accuracy, the size of the control volume may be chosen to be near or to be equal to the coarse mesh size or other suitable size if necessary. This flexibility of choosing the size of the control volume means that the constructed scheme can be not only applied to the flow problem for which the microstructure is completely found but may be also applied to the flow problem for which the microstructure is only partly found beforehand.

This study is limited to two-dimensional saturated flow through heterogeneous porous media. We also plan to extend this scheme to solve unsaturated water flow problems with heterogeneity which would be more difficult to simulate.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the Natural Science Foundation of China (51039007), the Natural Science Foundation of Hunan Province (13JJ3120), and the Construct Program of the Key Discipline in Hunan Province.