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

and Applied Analysis 3 Figure 1: Illustration of the macroscopic and microscopic computational domains. node (i, j), and F 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 themacroscopic flux F. In the absence of explicit knowledge of F, our problem reduces to approximate the flux F; 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 Iδ i±(1/2),j and Iδ i,j±(1/2) around every interior node (i, j) and these control volumes are centered at the points (i ± (1/2), j) and (i, j ± (1/2)), respectively (see Figure 2). We assume that the control volume Iδ i±(1/2),j is a square of size δ and discretize it into a fine M × M mesh, for which (ξ k , η l ) denotes the coordinate of node (k, l), where k, l = 1, . . . , M + 1. a = ξ k+1 − ξ k = η l+1 − η l denotes the fine mesh size. For ξ k and η l , we have ξ k = H − δ 2 + (k − 1) δ M , η l = H − δ 2 + (l − 1) δ M , k, l = 1, . . . , M + 1. (3) Similarly, (η l , ξ k ) denotes the coordinate of the control volume Iδ i,j±(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 K(x) 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 K(x) at the microscopic evolution. According to the conclusion of [14], in every control volume I, we can only solve the following reduced elliptic equation: ∇ ⋅ [K (x) ∇ψ (x)] = 0, in Iδ. (4) In every control volume I, 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 Iδ i+(1/2),j (i, j = 1, 2, . . . , N), the Dirichlet boundary condition is used in x-direction and noflow boundary condition is used in y-direction, and vice versa for Iδ i,j+(1/2) (i, j = 1, 2, . . . , N). Set the head of the basic elliptic problem with no dimensional change. Let φ1 and φ2 be the solutions of these two basic elliptic problems, respectively. In Iδ i+(1/2),j , 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 Iδ i,j+(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 differencemethodwithmultigrid over a finemesh to solve the original equation. For implementing themultigrid algorithm, we use directly a numerical simulator MGD9V [26]. In fact, according to the above two basic elliptic problems, we have φ 1 + φ 2 = 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 φ 2 = 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 I δ i+(1/2),j , Fx,β i+(1/2),j (α, β = 1, 2) denote, respectively, the basic macroscopic fluxes estimated by the solutions of two basic elliptic problems; we have


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][7][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][18][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 finescale 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.

New Scheme of FDHMM
where  is the hydraulic head,   is the specific storage coefficient, K is the hydraulic conductivity tensor,  is the source sink term, x = (, ) is the spatial coordinate,  is the time variable, Ω is the study area, and  is the time domain.

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), 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 Similarly, (  ,   ) denotes the coordinate of the control volume   ,±(1/2) .

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 K(x) 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 K(x) at the microscopic evolution.According to the conclusion of [14], in every control volume   , we can only solve the following reduced elliptic equation: 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 noflow 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  1 and  2 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 +  2 = 1.Then, we only need to solve the first basic elliptic problem and obtain the solution  2 of the second basic elliptic problem in the course of computation according to  2 = 1 −  1 .(i, j)

Estimation of
The control volumes at the coarse node (, ) and two basic microscopic elliptic models for the control volume where  ,+(1/2) is the geometric mean of  , and  ,+1 .

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 respectively, where Thus,  , , the solution of the microscopic elliptic problem in   +(1/2), , is obtained over the control volume   +(1/2), via the linear combination of  1 and  2 , which are the solutions of two corresponding basic elliptic problems, respectively, and then we have Similarly, in   ,+(1/2) , heads at the bottom and top sides are respectively, where The solution of the microscopic elliptic problem in this control volume is 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 By applying the above assumption, we derive the macroscopic flux   , 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 Similarly, for   ,+(1/2) , we have Combining ( 15), (16), and (17) yields where Interpolate heads at boundaries of local microscopic elliptic problems ( 9) and ( 12) Solve local basic microscopic elliptic problems (4) Estimate basic macroscopic fluxes (5)(6) Estimate macroscopic fluxes (16)(17) Solve the macroscopic model ( 18) Then, we solve the following equation at the coarse mesh by using MGD9V [26]: 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.
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 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].

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  2 norm and the relative maximum norm respectively, where   is the total number of nodes on the coarse mesh, Ψ  denotes the coarse solution at the coarse node x  , and Ψ(x  ) denotes the "exact" solution projected on the coarse mesh; that is, Ψ(x  ) = (x  ) at the coarse node.Here, (x) 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.

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.
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.
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,   = 10m, and  ln  = 2.0, and the relative  2 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.

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 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 2 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.
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) 2 unknowns; the degrees of freedom of these basic microscale problems are 2( + 1) 2 .Added ( + 1) 2 degrees of freedom of the macroscopic scheme, the total degrees of freedom of the new scheme are 2( + 1) 2 + ( + 1) 2 .Similarly, the total degrees of freedom of FDHMM are also 2( + 1) 2 + ( + 1) 2 , and degrees of freedom of the full fine scheme are (2 + 1) 2 .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.

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 3 /min, 0.24 m 3 /min, 0.36 m 3 /min, and 0.48 m 3 /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 3 /min, when  ln  = 0.5, 1.0, 1.5, 2.0, relative  2 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 3 /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  Relative maximum error 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.
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 −4 m −1 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 3 /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,   = 40m, and   = 10m.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 2 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.
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 2 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.

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.

Figure 1 :
Figure 1: Illustration of the macroscopic and microscopic computational domains.

Figure 3 :
Figure 3: Flow chart of the new scheme at a time step.

Figure 5 :
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.

Figure 6 :
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.

Figure 8 :
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.

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

Figure 10 :
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.

Figure 11 : 9 Q
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.

Figure 12 :
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.

Figure 15 :Figure 16 :
Figure 15: Relative (a)  2 and (b) maximum errors between exact and two coarse solutions for the transient flow problem with weak well drawdown.
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