Application of the multiscale finite element method to flow in ...

WATER RESOURCES RESEARCH, VOL. 40, W09202, doi:10.1029/2003WR002914, 2004

Application of the multiscale finite element method to flow in heterogeneous porous media Shujun Ye and Yuqun Xue Department of Earth Sciences, Nanjing University, Nanjing, China

Chunhong Xie Department of Mathematics, Nanjing University, Nanjing, China Received 30 November 2003; revised 27 May 2004; accepted 24 June 2004; published 9 September 2004.

[1] The multiscale finite element method is applied in this paper to simulate

groundwater flow in heterogeneous porous media. The method can efficiently capture the large-scale behavior of the solution without resolving all the small-scale features by constructing the multiscale finite element base functions that are adaptive to the local property of the differential operator. Both multiscale finite element method and conventional finite element method are applied to five 2-D and two 3-D groundwater flow problems, including a 2-D steady flow problem with continuous coefficients, a 2-D steady flow with gradual change in coefficients, a 2-D transient flow with gradual change in coefficients, a 2-D steady flow problem with an abrupt change in coefficients, a 2-D transient flow problem with an abrupt change in coefficients, and a 3-D steady and a 3-D transient flow problem with gradual change in coefficients in the horizontal direction and with abrupt change in the vertical direction, respectively. The applications demonstrate the main advantages of the multiscale finite element method, i.e., significantly reducing computational efforts and improving the accuracy INDEX TERMS: 1829 Hydrology: Groundwater hydrology; 3210 Mathematical of the solutions. Geophysics: Modeling; 3230 Mathematical Geophysics: Numerical solutions; KEYWORDS: heterogeneous, multiscale finite element method, numerical simulation of groundwater flow, porous media Citation: Ye, S., Y. Xue, and C. Xie (2004), Application of the multiscale finite element method to flow in heterogeneous porous media, Water Resour. Res., 40, W09202, doi:10.1029/2003WR002914.

1. Introduction [2] Most of the aquifer systems in nature are heterogeneous, and the heterogeneity of the subsurface formations spans over many scales. The groundwater flow problems in heterogeneous porous media can be solved more accurately by using conventional finite element method (FEM) or finite difference method based on smaller scale, which leads to more computational efforts. It is desirable to develop a method that can not only decrease the number of element but also achieve accurate results. The multiscale finite element method (MSFEM) applied in the paper can realize the purpose. [3] For the flow problems in heterogeneous porous media with coefficients spanning many scales, usually the large-scale features of the solution and the averaged effect of small scales on large scales are of the main interest. The flow problems in porous media are partially the elliptic problems in mathematics. The MSFEM used in the paper was developed to solve elliptic problem with highly oscillatory coefficients [Hou and Wu, 1997; Hou et al., 1999; Efendiev et al., 2000]. It can efficiently capture the large-scale behavior of the solution without resolving all the small-scale features by constructing the multiscale Copyright 2004 by the American Geophysical Union. 0043-1397/04/2003WR002914

finite element base functions that are adaptive to the local property of the differential operator. [4] The key difference between the MSFEM and conventional FEM is attributed to the base function. The base functions in the conventional FEM are constructed only according to coordinate values of the nodes of elements. They are linear in elements and cannot show the heterogeneity of the media in the elements. Parameters of an element are constants. However, the base functions in MSFEM are constructed by solving the reduced elliptic problem in each element. They can show the heterogeneity of the media in the elements and are nonlinear in elements. Parameters of an element can be variable. So for the flow problems in the highly heterogeneous porous media, the conventional FEM should be used on small scale, while the MSFEM can be used on large scale. [5] In addition to the MSFEM, many other methods were developed for the purpose of accuracy and efficiency, such as methods based on the homogenization theory [Dykaar and Kitanidis, 1992; Cruz and Petera, 1995] and some upscaling methods based on simple physical and/or mathematical motivation [Durlofsky, 1992; McCarthy, 1995]. These methods have some successful application, but they all have some limitation. The methods based on the homogenization theory are limited by restrictive assumptions on the media with scale separation and periodicity. The upscaling methods are more general and have been applied

W09202

1 of 9

W09202

YE ET AL.: APPLICATION OF MSFEM TO GROUND WATER FLOW

to problems with random coefficients with partial success [Dykaar and Kitanidis, 1992; Cruz and Petera, 1995], but the design principle is strongly motivated by the homogenization theory for periodic structures. Their applications to nonperiodic structures are not always guaranteed to work. These limitations are not present in the MSFEM. Babuska and Szymczak [1982], Babuska and Osborn [1983], and Babuska et al. [1994] also applied the idea of using base functions governed by the differential equations. However, these methods are based on the special property of the harmonic average in one-dimensional elliptic problems; thus they cannot explain the resonance effect occurred in genuinely multidimensional problems [Hou and Wu, 1997; Hou et al., 1999; Efendiev et al., 2000]. [6] Motivated by the numerical simulation of flow transport of two-phase flow in highly heterogeneous media, a mixed multiscale finite element (MMSFEM) was proposed by Chen and Hou [2002]. A multiscale coarse grid algorithm was developed to solve steady flow problem involving well singularities in heterogeneous porous medium based on the multiscale finite element, where the additional well singularities of the problem were resolved locally by adding finite element base functions [Chen and Yue, 2003]. [7] Some other multiscale approaches for flow in porous media have been proposed to accommodate the fine-scale description directly. Arbogast [2000, 2002, 2003] and Arbogast and Bryant [2002] presented a mixed finite element method for two-phase flow in highly heterogeneous media, where the fine-scale effects are localized by a boundary condition assumption at the coarse element boundaries. Then the small-scale influence is coupled with the coarse-scale effects by numerical Greens functions. Jenny et al. [2003] developed a multiscale finite volume method for elliptic problems describing flow in porous media, which is based on a flux-continuous finite difference approach and strongly resembles MSFEM by Hou and Wu [1997]. The effective coarse-scale transmissibilities are calculated to account for the fine-scale effects. Some numerical examples are provided in these papers, which show many similar computational results. [8] The MSFEM is a new numerical method for elliptic problems developed by mathematicians. It has been evaluated only sparsely in some literatures. There is a need for more evaluation of the applicability of this method. The MSFEM is used to simulate groundwater flow in heterogeneous porous media under many different conditions in the paper. We find it applicable not only to elliptic problems (steady groundwater flow problems) but also to parabolic problems (transient groundwater flow problems), not only to 2-D problems but also to 3-D problems. The accuracy and efficiency of this method show that it is very useful for groundwater flow simulation in the heterogeneous porous media. [9] This paper is organized as follows. The principle of MSFEM is introduced in section 2. In section 3, we present seven examples which are solved by the conventional FEM and MSFEM, respectively. They are (1) a 2-D steady flow problem in heterogeneous porous media with continuous change in coefficients; (2) and (3) about a 2-D steady flow problem and a 2-D transient flow problem in heterogeneous porous media with gradual change in coefficients, respectively; (4) and (5) about a 2-D steady flow problem and a

W09202

2-D transient flow problem in heterogeneous porous media with abrupt change in coefficients, respectively; (6) and (7) about a 3-D steady flow problem and a 3-D transient flow problem in heterogeneous porous media with gradual change in coefficients in the horizontal direction and with abrupt change in the vertical direction, respectively. Some similar numerical examples were dealt with other methods [Arbogast, 2000, 2002; Arbogast and Bryant, 2002; Arbogast, 2003]. However, for the MSFEM itself, Hou and Wu [1997] and Hou et al. [1999] only used MSFEM to compute two-dimensional steady flow problem with continuous coefficients. The results of these examples in the paper show the main advantages of MSFEM, such as significantly reducing the compute resources and improving accuracy. Section 4 is reserved for some conclusions.

2. Principle of MSFEM [10] MSFEM is developed to solve the elliptic problems with multiple scales. Steady groundwater flow problems in the heterogeneous porous media are examples of this type. The heterogeneity results in spatial variation in hydrogeologic parameters. It is difficult to obtain a direct numerical solution of the multiple-scale problems even with modern supercomputers. Direct solutions provide quantitative information of the physical processes at all scales. However, from an engineering perspective, it is often sufficient to predict the macroscopic properties of the multiple-scale systems. The MSFEM is a method to capture the small-scale effect on the large scales and obtain solutions on the large-scale accurately and efficiently. [11] MSFEM is very useful for flow simulation in a large region with highly heterogeneous porous media. Parameters of each element are constants in the conventional FEM, thus only fine grid mesh can describe the variation of parameters. It would need many computing resources if the study area were very large. Parameters in an element can vary in MSFEM, and the variation of parameters is brought to the base functions. Thus coarse grid mesh can describe the variation of parameters, which results in much less computing resources. [12] The main idea of MSFEM is to construct finite element base functions which capture the small-scale information within each element. The characteristic difference between MSFEM and the conventional FEM is attributed to base function. 2.1. Base Functions in MSFEM [13] We consider solving the second-order elliptic equation r Kðx; yÞrH ¼ f

ðx; yÞ 2 D

ð1Þ

where r is Hamilton operator, K(x, y) is hydraulic conductivity, H is hydraulic head, f is source sink term, D is the study area. Assume the element is triangular, we introduce how to construct the base functions for an element Dijk shown in Figure 1. fi denotes the base function of node i, which satisfies

2 of 9

r Kðx; yÞrfi ¼ 0

ðx; yÞ 2 Dijk

ð2Þ

W09202


W09202

boundary condition. On side ij, we have the reduced problem: @ @f KðxÞ i ¼ 0 @x @x

ð4Þ

fi j i ¼ 1 fi j ¼ 0

ð5Þ ð6Þ

The problem can R x be solved analytically, that is Rx fi ðxÞij ¼ x j KdtðtÞ= xij KdtðtÞ. If K is constant in side ij, then fi(x) = (xj x)/(xj xi) is linear.

Figure 1. Triangular element Dijk. Similarly for fj and fk. The boundary conditions of base function fi are: fiji = 1, fijj = 0, fijk = 0, fi varies from 1 at node i to 0 at node j on edge ij, fi varies from 1 at node i to 0 at node k on edge ik, and fi is 0 on edge jk. The boundary conditions of fj and fk are similar to those of fi. fi, fj and fk satisfy fi + fj + fk = 1. In general, the base functions have no analytic solutions. The numerical solutions of base functions are obtained by numerical method such as FEM. The construction of the base function is fully decoupled from element to element; thus the method is perfectly parallel. It will take some computing resources to calculate the base functions. If the size of the element for the conventional finite element method is the same as that of subcell for multiscale base functions, the total cost of MSFEM will be much less than that of the conventional finite element method [Hou and Wu, 1997]. The accuracy of the final results is relatively insensitive to the accuracy of the base functions [Hou et al., 1999]. Head at any point in an element satisfies

2.3. Oversampling Method of the MSFEM [16] To avoid the resonance between the mesh scale and the physical scale [Hou and Wu, 1997] and to improve the accuracy and convergence rate for the MSFEM, Hou and Wu [1997] proposed the oversampling method to obtain the base functions from the sampling element which is larger than the origin element. [17] Consider the element in Figure 2, Dijk is the origin element, and Dabc is the sampling element. Firstly, the temporary base functions of the sampling element is constructed by the way presented in sections 2.1 and 2.2. Denote the temporary base function at node a, b, c of Dabc as y1, y2 and y3, respectively. Denote the actual base function at node i, j and k of Dijk as f1, f2, and f3, respectively. The actual base functions are constructed from the linear combination of the temporary base functions as follows: fi ðx; yÞ ¼

3 X

Ci;j yj ðx; yÞ

ði ¼ 1; 2; 3Þ

ð7Þ

j¼1

where Ci,j are the constants determined by the condition fi (xj, yj) = dij(dij = 1, if i = j, dij = 0, if i 6¼ j). For example, to obtain f1, C1,1, C1,2 and C1,3 are the solutions of the following equation set: f1 ðxi ; yi Þ ¼ C1;1 y1 ðxi ; yi Þ þ C1;2 y2 ðxi ; yi Þ þ C1;3 y3 ðxi ; yi Þ ¼ 1 ð8Þ

f1 xj ; yj ¼ C1;1 y1 xj ; yj þ C1;2 y2 xj ; yj þ C1;3 y3 xj ; yj ¼ 0

Hðx; yÞ ¼ Hðxi ; yi Þfi ðx; yÞ þ H xj ; yj fj ðx; yÞ þ Hðxk ; yk Þfk ðx; yÞ

ð9Þ

ð3Þ

f1 ðxk ; yk Þ ¼ C1;1 y1 ðxk ; yk Þ þ C1;2 y2 ðxk ; yk Þ þ C1;3 y3 ðxk ; yk Þ ¼ 0 ð10Þ

which is similar to that in the conventional FEM. 2.2. Boundary Condition of Base Functions [14] Although the final results of flow problems are relatively insensitive to the accuracy of the base functions, they are sensitive to the boundary condition of base functions. There are two kinds of boundary conditions. One choice is linear boundary condition, which means base functions are linear along boundary, such as fi vary linearly along side ij from fiji = 1 at node i to fijj = 0 at node j. The linear boundary condition is the same as that in the conventional FEM. Another choice is oscillatory boundary condition. Solve the reduced elliptic problems on each side of boundary with values 1 and 0 at the two end points, and use the resulting solution as the boundary condition for the base function. The oscillatory boundary condition can significantly improve the accuracy of the MSFEM. [15] Take the side ij of the element in Figure 1 as an example, we introduce how to obtain the oscillatory 3 of 9

Figure 2. Sampling element.

W09202


The rationale behind the MSFEM was described in detail by Hou and Wu [1997] and Hou et al. [1999]. [18] We demonstrate the construction of base function for 2-D steady flow problem in the section. Although the MSFEM is originally developed to deal with 2-D problems, it can generalize to 3-D problems [Farmer, 2002]. The construction of base function for 3-D steady flow problem is similar to that for 2-D problem. Since for the transient problem the main difficulty is the same as that for the steady state problem, i.e., the multiple scales in the solution. The MSFEM can be easily extended to solve the transient problems [Hou and Wu, 1997]. We use the same base function for transient problem as that for steady problem in the paper. [ 19 ] To facilitate the comparison among different schemes, we use the following abbreviations: LFEM stands for the conventional finite element method, LFEM-F stands for the conventional finite element method with fine grid mesh, MSFEM-L stands for the multiscale finite element method with linear boundary condition for base functions, MSFEM-O stands for the multiscale finite element method with oscillatory boundary condition for base functions, and MSFEM-OS-O stands for the multiscale finite element method with oversampling method and oscillatory boundary condition for base functions.

3. Application of MSFEM to Flow in Heterogeneous Porous Media [20] All porous media in nature is heterogeneous. The heterogeneity is represented by the spatial variation of hydrogeologic parameters. The spatial variation of parameters of media formed under different conditions are different. Parameters of some media vary continuously. Parameters of some media vary discontinuously, including gradual change and abrupt change. MSFEM can solve flow problems in heterogeneous porous media very well. In this section, we will use seven examples to show the main advantages of MSFEM. 3.1. Two-Dimensional Steady Flow in Heterogeneous Porous Media With Continuous Change in Coefficients [21] Continuous change in coefficients here is specified as the coefficients change according to the smooth function. Consider a 2-D steady flow problem with continuous coefficient as follows: [22] The study area is a rectangle. Within the rectangle, the results satisfy: @ @H @ @H K þ K ¼0 @x @x @y @y

ð11Þ

Kðx; yÞ ¼ x2

ð12Þ

where

There will be many solutions for equation (11) if without boundary condition. We assume the analytic solution is H ¼ x2 3y2

ð13Þ

We assume the boundary condition of the study problem is first kind of boundary condition, which is determined by equation (13).

W09202

Figure 3. The errors of various numerical solutions in section y = 12.5 m for study area D5 under the condition of 2-D steady flow with continuous coefficient. [23] Given seven study areas, corresponding x coordinate ranges are: [0.001, 0.1], [0.1, 0.5], [0.5, 1], [1, 5], [5, 20], [20, 50] and [50, 150], respectively. We denote the seven study areas as D1 – D7 in turn. Seven x coordinate ranges correspond to seven ranges of hydraulic conductivity. Hydraulic head at any point in the study area can be calculated by equation (13). In addition, the problem is solved using LFEM, LFEM-F, MSFEM-L, MSFEM-O and MSFEM-OS-O, respectively. The study areas are subdivided into 1152 triangular elements for LFEM-F, and 72 triangular elements for all other methods. There are 16 triangular subcells in each element for the MSFEM-L, MSFEM-O and MSFEM-OS-O. The size of subcell in MSFEM is the same as that of element in LFEM-F. [24] The solutions of groundwater flow in the study areas D5 and D7 are used as examples to compare accuracy of different numerical method. Figure 3 shows errors of numerical solutions versus x coordinate at each point in section y = 12.5 m in D5 with x coordinate range [5,20]. Figure 4 shows the same relationship as Figure 3 in section y = 100 m in D7 with x range [50,150]. We observe that the errors of solutions of LFEM are the largest; those of MSFEM-L are only less than those of LFEM; results of MSFEM-O and LFEM-F have about the same accuracy; MSFEM-OS-O solutions are the most accurate. The results of different method for other study area are similar to those illustrated in Figures 3 and 4. The results of MSFEM-O at large scale are almost equal to those of LFEM-F at small scale, which indicates that the MSFEM can not only save computing resources but also obtain accurate results. An interesting phenomenon is that the solutions of MSFEMOS-O can be more accurate than those of LFEM-F. One would think that the resolution of the direct solution on a fine grid should be higher than that of the MSFEM on a coarser grid. This phenomenon was also mentioned by Hou and Wu [1997]. [25] The size of sampling element in the MSFEM-OS-O is 1.1 times the size of the origin element in Figures 3 and 4.

4 of 9

W09202


Figure 4. The errors of various numerical solutions in section y = 100 m for study area D7 under the condition of 2-D steady flow with continuous coefficient. The size of sampling element changes, and then the accuracy of solution changes. For the study area D5, we assume the size of sampling element is 1, 1.001, 1.01 and 1.1 times the size of origin element, respectively. The results in section y = 12.5 m calculated by MSFEM-OS-O with different size of sampling element are shown in Figure 5. In Figure 5 the larger the size of sampling element is, the more accurate the results are. The results in other study areas are similar. However, oversampling technique increases the computing resources requirement, and the larger the size of sampling element is, the more computing resources are needed. For flow simulation in a large region, the mesh scale is always much larger than the physical small

Figure 5. The errors of solutions using MSFEM-OS-O with different size of sampling domain.

W09202

Figure 6. The relative errors of flow rates between analytic method and various numerical methods. scale of the heterogeneous media, and no resonance exists. For efficiency and accuracy reasons, MSFEM-O fits for flow problem in practical large region relative well. [26] In addition, we calculate flow rate using LFEM, LFEM-F and MSFEM-O and compare results of these methods. For the 2-D steady flow we designed the flow direction to be about from the bottom up. Here we compare the flow rate through the middle section of each study area calculated by different method (see Figure 6). The middle section of study area D5 is section y = 12.5 m, that of study area D6 is section y = 35 m and so on. Figure 6 shows the absolute relative error between analytic solutions and solutions of various numerical methods versus y coordinate value of the middle section of each study area. We find that the flow rate obtained by MSFEM-O is more accurate than that obtained by LFEM. 3.2. Two-Dimensional Steady Flow Problem in Heterogeneous Porous Media With Gradual Change in Coefficients [27] The study area is a rectangle covering 10 km 10 km with point (0,0) as the origin of coordinate. Hydraulic conductivities from left to right in the study area vary gradually from 1 m d1 to 250 m d1. The value of hydraulic conductivity increases 1 m d1 every 40 m from left to right. The porous media with gradual change in coefficients exists in the alluvial plain. The left and right boundaries are first kind of boundaries. Head on the left side is 10 m, and head on the right side is 0 m. Top and bottom sides are impermeable boundaries. The problem has no analytic solution, so the fine-scale solution obtained by LFEM-F is regarded as the exact solution. In all numerical examples below, the LFEM-F solution is treated as the exact solution. The study area is divided into 125,000 and 1250 triangular elements for LFEM-F and other all methods, respectively. Heads in section y = 5200 m obtained from different methods are compared in Figure 7. Figure 7 shows head versus x coordinate for each point in the section. We observe the results of MSFEM-O are closest to those of

5 of 9

W09202


Figure 7. The results of different methods under the condition of 2-D steady flow with gradual change in coefficient. LFEM-F, the results of MSFEM-L are less accurate than those of MSFEM-O, and the results of LFEM are the worst. 3.3. Two-Dimensional Transient Flow Problem in Heterogeneous Porous Media With Gradual Change in Coefficients [28] The transient problem has the same shape, range, parameters distribution, and boundary condition of the study area as the problem in section 3.2. In addition, we assume the specific storages from left to right in the study area vary gradually from 1 105 m1 to 1 106 m1. The aquifer is 10 m thick. There is a pumping well at point (5200, 5200). The well has the constant flow rate 1000 m3 d1, and is pumped 5 days in the problem. The initial condition is heads varying linearly from 10 m on the left boundary to 0 m on the right boundary. The fine-scale solution obtained by LFEM-F is still regarded as the exact solutions. The spatial discretization for every numerical method is the same as that in section 3.2. The time step size is 1 day for every method. Figure 8 shows the results in section y = 5200 m of LFEM, LFEM-F, MSFEM-L and MSFEM-O after pumping 5 days. It shows head versus x coordinate for each point in the section. We also find that the results of MSFEM-O are closest to those of LFEM-F, the results of MSFEM-L are less accurate than those of MSFEM-O, and the results of LFEM are the worst. There are larger errors of the results of MSFEM-O near point (5200,5200), which are caused by the pumping well at that point. The heads near the well vary logarithmically with distance to the well, which cannot be well described by the present method. Chen and Yue [2003] proposed a modified MSFEM to resolve the well singularities. Better results can be obtained if there are finer-grid meshes near wells, as done by Arbogast [2002]. 3.4. Two-Dimensional Steady Flow Problem in Heterogeneous Porous Media With Abrupt Change in Coefficients [29] The study area is also a rectangle covering a range of 10 km 10 km with point (0,0) as the origin of coordinate.

W09202

The left and right sides of boundary are first kind of boundaries. Head on the left side is 10 m, and head on the right side is 0 m. Top and bottom sides are impermeable boundaries. There is an abrupt interface in the study area at x = 2480 m, and hydraulic conductivity in the left side of the interface is 2 m d1 and that in the right side of the interface is 1000 m d1. The porous media with abrupt interface is general in the field. The study area is divided into 125000 triangular elements for LFEM-F, and 200 triangular elements for all other methods. Heads on interface y = 5000 m calculated using various methods are compared in Figure 9. Figure 9 shows absolute relative error between solutions of LFEM-F and those of other methods versus x coordinate for each point in the interface. In Figure 9 the results of MSFEM-O are the best, the results of MSFEM-L are close to those of LFEM, and both are bad. [30] The reason that the MSFEM-O can obtain more accurate solutions is that the base functions can capture the small-scale information. The base functions vary nonlinearly with parameters’ variation. Take the element Dijk, which includes the abrupt interface, for example; the base functions of MSFEM and LFEM are shown in Figures 10a and 10b. The coordinates of node i, j, k of the element Dijk in Figure 10 are (2000, 0), (3000, 0) and (2000, 1000), respectively. Figure 10a shows the base function at node i of MSFEM-O, which is nonlinear in the element because of the abrupt change in coefficient. Figure 10b shows the base function at node i of LFEM, which is linear in the element. 3.5. Two-Dimensional Transient Flow Problem in Heterogeneous Porous Media With Abrupt Change in Coefficients [31] The transient problem has the same shape, range, parameters distribution, and boundary condition of the study area as the problem in section 3.4. There are 3 pumping wells at point (5000, 0), (5000, 5000), and (5000, 10000), respectively. The three wells have the constant flow rate 10,000 m3 d1, and all are pumped 5 days in the

Figure 8. The results of different methods under the condition of 2-D transient flow with gradual change in coefficient.

6 of 9

W09202


W09202

Head on the left side is 10 m, and head on the right side is 1 m. Other sides are impermeable boundaries. The study area is divided into 14 layers in the vertical direction for LFEM-F, which means the aquifers are separated from the aquitards. Every horizontal layer is divided into 1600 hexahedral elements. So the study area is divided into 22,400 hexahedral elements and 25215 nodes for LFEM-F. The study area is divided into two layers in the vertical direction for LFEM, MSFEM-L and MSFEM-O, that means

Figure 9. The relative errors between the results of LFEM-F and those of other methods under the condition of 2-D steady flow with abrupt change in coefficient. problem. Hydraulic conductivity is 2 m d1 in the left side of the interface and 1000 m d1 in the right side. Specific storage is 0.000002 m1 in the left side and 0.001 m1 in the right side. The aquifer is 10 m thick. The initial condition is heads varying linearly from 10 m on the left boundary to 0 m on the right boundary. LFEM, LFEM-F and MSFEM-O methods are used to solve the problem. The finite element spatial discretization for every numerical method is the same as that in section 3.4. The time step size is 1 day for every method. Figure 11 shows the results in interface y = 5000 m of LFEM, LFEM-F and MSFEM-O after pumping 5 days. It shows head versus x coordinate for each point in the interface. We find that the results of the MSFEM-O are relative better for the transient problem with abrupt changes of coefficients. There are larger errors of the results of MSFEM-O near the well at point (5000, 5000), which are caused by the same reason as that mentioned in section 3.3. 3.6. A 3-D Steady Flow Problem in Heterogeneous Porous Media With Gradual Change in Coefficients in the Horizontal Direction and With Abrupt Change in the Vertical Direction [32] The study area is a rectangular domain covering 10 km 10 km 140 m with point (0, 0, 0) as the origin of coordinate. The simulated aquifer system is composed of seven confined aquifers and seven aquitards. The aquifers and aquitards distribute by turn in the vertical direction. All the aquifers and aquitards are 10 m thick. They are all heterogeneous. Hydraulic conductivities from left to right in all aquifers vary gradually from 1 m d1 to 196 m d1, that means they increase 5 m d1 every 250 m from left to right. Hydraulic conductivities from left to right in all aquitards vary gradually from 0.005 m d1 to 0.98 m d1, that means they increase 0.025 m d1 every 250 m from left to right. Thus the hydraulic conductivity changes abruptly in the vertical direction from the aquifer to the aquitard. The left and right sides of boundary are first kind of boundaries.

Figure 10. Base functions of MSFEM-O and LFEM. (a) Base function at node i of an element in MSFEM-O. (b) Base function at node i of an element in LFEM.

7 of 9

W09202


Figure 11. The results of different methods under the condition of 2-D transient flow with abrupt change in coefficients.

four aquifers and three aquitards in the first layer and four aquitards and three aquifers in the second layer. Every horizontal layer is divided into 400 hexahedral elements. The study area is divided into 800 hexahedral elements and 1323 nodes for these methods. Every element in MSFEM-L and MSFEM-O is subdivided into seven layers in the vertical direction, which means the aquifers are separated from the aquitards. Every horizontal layer in an element is divided into 4 hexahedral elements. Then, the element is subdivided into 28 hexahedral elements for the two methods. Heads on section y = 5000 m and z = 70 m calculated using various methods are compared in Figure 12. Figure 12 shows errors between solutions of LFEM-F and those of other methods versus x coordinate for each point in the interface. In Figure 12 the errors of results of LFEM are the largest, those of MSFEM-L are less, and those of MSFEM-O are the least. The results of MSFEM-O are very close to those of LFEM-F. 3.7. A 3-D Transient Flow Problem in Heterogeneous Porous Media With Gradual Change in Coefficients in the Horizontal Direction and With Abrupt Change in the Vertical Direction [33] The transient problem has the same shape, range, parameters distribution, and boundary condition of the study area as the problem in section 3.6. There is 1 pumping well at point (5000,5000,70). The well has the constant flow rate 3000 m3 d1, and it is pumped 50 days in the problem. The specific storages from left to right in all aquifers and aquitards vary gradually from 1 105 m1 to 2.2 106 m1, which means they decrease 2 107 m1 every 250 m from left to right. The initial condition is heads varying linearly from 10 m on the left boundary to 1m on the right boundary. LFEM, LFEM-F, MSFEM-L and MSFEM-O are used to solve the problem. The finite element spatial discretization for every numerical method is the same as that in section 3.6. The time step size is 10 days for every method. Heads on section y = 5000 m and

W09202

Figure 12. The errors between the results of LFEM-F and those of other methods under the condition of 3-D steady flow. z = 70 m calculated using various methods are compared in Figure 13. Figure 13 shows errors between solutions of LFEM-F and those of other methods versus x coordinate for each point in the interface. Figure 13 also indicates that the results of MSFEM-O are the best. However, there are larger errors of the results of MSFEM-O near the well at point (5000, 5000, 70), which are caused by the same reason as that mentioned in section 3.3.

4. Conclusions [34] We use different numerical method to solve seven groundwater flow problems characterized by heterogeneous

Figure 13. The errors between the results of LFEM-F and those of other methods under the condition of 3-D transient flow.

8 of 9

W09202


coefficients. Our numerical experiments give convincing evidence that the MSFEM is effective to solve these problems accurately. It is capable of capturing the largescale solution without resolving the small-scale details, so that it saves computing efforts. It is effective to deal with not only the problems with continuous coefficients but also those with gradual change in coefficients and abrupt change in coefficients. It is effective to deal with not only the elliptic problems (the steady flow problems) but also the parabolic problems (the transient flow problems), not only the 2-D problems but also the 3-D problems. Different boundary condition of the base function significantly influences the accuracy of the MSFEM. The results using oscillatory boundary condition for base functions are much more accurate than those using linear boundary condition, so the former boundary condition is preferred in application. The results using oversampling technique are more accurate than those not using the technique. The larger the size of sampling domain is, the more accurate the results are. However, larger size of sampling domain leads to more computing resources requirement. In reality, the mesh scale is much coarser than the physical small scales, thus no resonance exists. For the goal of efficiency and accuracy, the multiscale finite element with oscillatory boundary and without oversampling technique is very useful for the very large practical problems with heterogeneous properties. The characteristic difference between MSFEM and the conventional FEM is attributed to the base function. The base function of MSFEM can indicate the variation of parameters in an element, so that MSFEM has advantages in dealing with the heterogeneity. However MSFEM is not able to decrease the errors caused by other factors. For example, the error of solutions of MSFEM is larger near the pumping well as mentioned in sections 3.3, 3.5 and 3.7. The well is a singular point. The base functions of elements near the well need to be modified, which will be considered in a separate paper. The MSFEM is developed for the elliptic problem, but we find it also can be relative well applied to solve the parabolic problem. We use the reduced elliptic problem to obtain the base functions for parabolic problem in the paper, however, we will study on calculating the base functions using reduced parabolic problem next step. We have demonstrated the application of MSFEM to some numerical examples. We will apply MSFEM to solve the 3-D land subsidence model in Shanghai in another paper, which covers 6400 km2 and has 12 highly heterogeneous aquifers and aquitards. The parameter values are known and assigned to the fine grid in these numerical examples in the paper. However, parameter values of large practical groundwater flow problem like that in Shanghai are unknown. They will be estimated by field tests and model calibration. [35] Acknowledgments. We thank the anonymous reviewers for many constructive comments. This paper is financially supported by the

W09202

National Natural Science Foundation of China grants 40172082 and 40335045.

References Arbogast, T. (2000), Numerical subgrid upscaling of two-phase flow in porous media, in Numerical Treatment of Multiphase Flows in Porous Media, Lect. Notes Phys., vol. 552, edited by Z. Chen, R. E. Ewing, and Z. C. Shi, pp. 35 – 49, Springer-Verlag, New York. Arbogast, T. (2002), Implementation of a locally conservative numerical subgrid upscaling scheme for two-phase Darcy flow, Comput. Geosci., 6, 453 – 481. Arbogast, T. (2003), An overview of subgrid upscaling for elliptic problems in mixed form, in Current Trends in Scientific Computing, edited by Z. Chen, R. Glowinski, and K. Li, pp. 21 – 32, Am. Math. Soc., Providence. Arbogast, T., and S. L. Bryant (2002), A two-scale numerical subgrid technique for waterflood simulations, SPE J., 446 – 457. Babuska, I., and E. Osborn (1983), Generalized finite element methods: Their performance and their relation to mixed methods, SIAM J. Numer. Anal., 20, 510 – 536. Babuska, I., and W. G. Szymczak (1982), An error analysis for the finite element method applied to convection-diffusion problems, Comput. Methods Appl. Math. Eng., 31, 19 – 42. Babuska, I., G. Caloz, and E. Osborn (1994), Special finite element methods for a class of second order elliptic problems with rough coefficients, SIAM J. Numer. Anal., 31, 945 – 951. Chen, Z., and T. Y. Hou (2002), A mixed multiscale finite element method for elliptic problems with oscillating coefficients, Math. Comput., 72, 541 – 576. Chen, Z., and X. Y. Yue (2003), Numerical homogenization of well singularities in the flow transport through heterogeneous porous media, Multiscale Model. Simul., 1(2), 260 – 303. Cruz, M. E., and A. Petera (1995), A parallel Monte-Carlo finite-element procedure for the analysis of multicomponent random media, Int. J. Numer. Methods Eng., 38, 1087 – 1121. Durlofsky, L. J. (1992), Representation of grid block permeability in coarse scale models of randomly heterogeneous porous media, Water Resour. Res., 28(7), 1791 – 1800. Dykaar, B. B., and P. K. Kitanidis (1992), Determination of the effective hydraulic conductivity for heterogeneous porous media using a numerical spectral approach: 1. Method, Water Resour. Res., 28(4), 1155 – 1166. Efendiev, Y. R., T. Y. Hou, and X. H. Wu (2000), Convergence of a nonconforming multiscale finite element method, SIAM J. Numer. Anal., 37, 888 – 910. Farmer, C. L. (2002), Upscaling: A review, Int. J. Numer. Methods Fluids, 40, 63 – 78. Hou, T. Y., and X. H. Wu (1997), A multiscale finite element method for elliptic problems in composite materials and porous media, J. Comput. Phys., 134, 169 – 189. Hou, T. Y., X. H. Wu, and Z. Cai (1999), Convergence of a multiscale finite element method for elliptic problems with rapidly oscillating coefficients, Math. Comput., 68, 913 – 943. Jenny, P. S., H. Lee, and H. A. Tchelepi (2003), Multi-scale finite-volume method for elliptic problems in subsurface flow simulation, J. Comput. Phys., 187, 47 – 67. McCarthy, J. F. (1995), Comparison of fast algorithms for estimating largescale permeabilities of heterogeneous media, Transp. Porous Media, 19, 123 – 137.

C. Xie, Department of Mathematics, Nanjing University, Nanjing 210093, China. Y. Xue and S. Ye, Department of Earth Sciences, Nanjing University, Nanjing 210093, China. ([email protected])

9 of 9