A domain decomposed finite element method for solving Darcian velocity in heterogeneous porous media

Journal of Hydrology 554 (2017) 32–49

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

A domain decomposed finite element method for solving Darcian velocity in heterogeneous porous media Yifan Xie a, Jichun Wu a,⇑, Yuqun Xue a, Chunhong Xie b, Haifeng Ji c a

School of Earth Sciences and Engineering, Nanjing University, Nanjing, China Department of Mathematics, Nanjing University, Nanjing, China c School of Science, Nanjing University of Posts and Telecommunications, Nanjing, China b

a r t i c l e

i n f o

Article history: Received 24 April 2017 Received in revised form 7 August 2017 Accepted 22 August 2017 Available online 26 August 2017 This manuscript was handled by P. Kitanidis, Editor-in-Chief, with the assistance of Jian Luo, Associate Editor Keywords: Finite element method Domain decomposition technique Darcian velocity Refraction law

a b s t r a c t This paper proposes a domain decomposed finite element method (DDFEM) for groundwater flow velocity simulation, which can effectively and efficiently deal with arbitrarily oriented or intersected material interfaces in heterogeneous porous media. The main idea of this problem is to employ domain decompose technique to break the velocity problem down to subproblems by material interfaces, thereby achieving high accuracy at interface nodes and improving the velocity computation efficiency. By solving the global flow field by FEM before velocity computation, this method can capture the global information of the whole study region so as to ensure the subproblem solution accuracy by Yeh’s model. After one subproblem has been solved, the DDFEM employs refraction law to obtain the Dirichlet boundary velocities of the adjacent subdomains. Numerical examples have been done to demonstrate the applicability and accuracy of the proposed method. Comparison between the DDFEM and two classical methods shows that the DDFEM can indeed save much computational cost, while achieving high accuracy. Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction Most nature groundwater flow system is heterogeneous, the heterogeneity is often associated with features, such as different geological materials, fractures or inclusions. Often such a system can be divided into several domains by the physical material interfaces. At the interface of different materials, obeying the refraction law (Bear, 1972; Bear, 1979), the hydraulic heads and the velocities normal to interfaces should be continuous, and the tangent velocities should be discontinuous (Zhou et al., 2001). However, the continuous feather of normal velocity leads to difficulty for many classical techniques (Cainelli et al., 2012), such as finite volume methods (FVM) (Tompson and Gelhar, 1990; Chin and Wang, 1992) and finite element methods (FEM) (Bellin et al., 1992; Neuman et al., 1992), when modeling in heterogeneous media. Focused on the continuity of velocity, many continuous methods have been developed. The mixed finite element methods regards both head and velocity as degrees of freedom, which can ensure the velocity continuity and achieves high accuracy (Meissner, 1973; Segol et al., 1975; Chen and Hou, 2003; D’Angelo and Scotti, 2012; Ervin, 2013). And a mixed-hybrid finite element model is presented by Chavent and Roberts (Chavent and ⇑ Corresponding author. E-mail address: [email protected] (J. Wu). http://dx.doi.org/10.1016/j.jhydrol.2017.08.040 0022-1694/Ó 2017 Elsevier B.V. All rights reserved.

Roberts, 1991). Yeh develops a Galerkin finite element model for computing continuous velocity field after the heads distribution obtained by FEM(Yeh, 1981). Unlike the mixed finite element methods, Yeh’s model can ensure the positive definite feather of equation system and does not require iteration (Zhou et al., 2001; Zhang et al., 1994). In contrast to other head-calculationfirst continuous velocity methods, such as dual mesh method (Batu, 1984) and cubic spline method (Zhang et al., 1994), this method is more accurate and applicable. Up to now, Yeh’s model has been applied or developed by many researches (Boufadel et al., 1999; Jang and Aral, 2007; Park and Aral, 2007; Xie et al., 2016) to ensure the global continuity of the velocity field, showing the high effectiveness and applicability of this method. However, same as the other continuous methods, Yeh’s model does not satisfy the discontinuous part of the refraction law, so that it suffers inaccuracy in the vicinity of material interfaces (Zhou et al., 2001; Cainelli et al., 2012). To resolve this problem, some modifications of Yeh’s model are proposed to achieve better velocity field in the vicinity of material interfaces. Srivastava and Brusseau introduce a continuous finite element scheme for normal velocities and a discontinuous scheme for tangent velocities (Srivastava and Brusseau, 1995), but only consider the material interfaces parallel to a coordinate axis. Zhou et al. introduce a jump function (obtained by refraction law) into the equation system of Yeh’s model by iteration, which can ensure

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

33

Fig. 1. The flow chart of DDFEM procedure.

the accuracy of velocities at interface nodes (Zhou et al., 2001). Besides the modifications of Yeh’s model, some other methods can also deal with material interfaces well. Based on numerical manifold method, Hu et al. presented a continuous approach with jump functions and a discontinuous approach with lagrange multipliers to deal with arbitrarily oriented or complexly intersected material interfaces, which do not require conforming mesh and do not require iterations (Hu et al., 2015). Moreover, some works focus on fractures which also involve good interface velocity approximations (Hoteit and Firoozabadi, 2008; Martin et al., 2005; Ahmed et al., 2015), showing the importance of the study in dealing with material interfaces. However, since groundwater flow problems become more and more complex, the computational cost of velocity simulation grows rapidly. The computational efficiency is also very important for today’s numerical simulation. It is necessary to develop a simple and straightforward method to accurately deal with material interfaces in high efficiency. Our purpose of this paper is to propose a domain decomposed finite element method (DDFEM), which can accurately and efficiently solve velocity in heterogeneous porous media with material interfaces. Instead of solving the velocity in the whole study region, the domain decomposition technique allows DDFEM to break the velocity computation down to subproblems with the material interfaces as the boundaries. Before computing the velocity, the DDFEM requires to solve the global flow field in the whole study region to obtain continuous hydraulic heads so as to capture the global information by regarding the heads as the solution condition of subproblems, which is different from the classical domain decomposition. Unlike multiscale methods, the DDFEM velocity field is not computed in coarse scale. The subproblems are solved by Yeh’s model in one material subdomains in fine scale, so that different materials can’t affect each other. Therefore, DDFEM can obtain accuracy velocities at interface nodes, and the velocity continuity, mass balance and solution converge are ensured by Yeh’s model. After the one subproblem is solved, the DDFEM can employ the refraction law to obtain the boundary velocities of the adjacent subdomains, which can be regarded as the Dirichlet boundary condition of the corresponding adjacent subproblems. Therefore, the

DDFEM velocities satisfy the refraction law at material interface nodes. Altogether, the DDFEM can not only achieve high accuracy in dealing with material interfaces, but also have high computational efficiency. Furthermore, the DDFEM is more simple and straightforward than most existing interface methods, which is easy for practical application. In this paper, the DDFEM is compared with two classical methods: Yeh’s Galerkin model and Zhou’s S1 scheme in three examples, including vertical, oblique or intersected material interfaces and isotropic or anisotropic cases. Same as Zhou’s scheme, DDFEM achieves high accuracy in deal with different types interfaces, so that it is more accurate than Yeh’s model. Furthermore, the DDFEM employs the domain decomposition technique to transform the high order velocity equation system into serval low order equation systems, and does not require iteration. Therefore, DDFEM costs much less than the two Zhou’s scheme and Yeh’s Galerkin model. Moreover, the mass conservation of the DDFEM is satisfied to an acceptable degree. In numerical examples, the CPU time of DDFEM is about a half of that of Yeh’s model and less than 10% of that of Zhou’s scheme. Consequently, the DDFEM is more efficient than the two classical methods, while achieving high accuracy at material interface nodes. This paper is organized as follows: The next section introduces the refraction law. Then the process of DDFEM is introduced in Section 3, including an implementation for a 2D problem and a computational cost analysis. Then, three examples are proposed to compare DDFEM with Yeh’s Galerkin model and Zhou’s S1 scheme under different situations. In the end is the conclusion. 2. The refraction law At the interface between two different materials, the velocity and head satisfy the refraction law:

Vþ  n ¼ V   n

ð1Þ

Hþ ¼ H ;

ð2Þ

34

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

where V is the Darcian velocity, n is the unit normal vector of the interface, H is hydraulic head, the superscripts þ and  stand for the values on both sides of the interface. Since the heads are continuous at the interface nodes, the gradient vector J follows that

Jþ  s ¼ J  s;

ð3Þ

where s is the unit tangential vector of the interface. According to the Darcy’s Law, Jþ and J satisfy that 1

 Vþ

 1



Jþ ¼½Kþ  

J ¼½K 

V

ð4Þ

where K is hydraulic conductivity tensor. According Eqs. (1)–(4), the refraction law can be represented as:

V  n ¼ Vþ  n  1

½K 



ð5Þ þ 1

 V  s ¼ ½K 

þ

V s

Therefore, the jump form of the refraction law can be obtained

" Vjump ¼

V þx  V x V þy  V y

#

 ¼

nT sT

  



0

ð6Þ

1

ð1  ½K   ½Kþ  Þ  Vþ  s

where V x and V y are velocities in x; y direction, respectively. 3. The domain decomposed finite element method (DDFEM) The process of DDFEM is introduced in this section, with Fig. 1 as the flow chart. The implementation of DDFEM for 2D steady flow problem is stated in Section 3.1. Then a computational cost comparison between DDFEM and Yeh’s model is stated in SecT

The DDFEM employs the FEM to solve (9) to obtain the global head distribution so as to capture the global information of the study region X, including the material interface effects. And in each element Mijk in X, the hydraulic head H satisfies:

Hðx; yÞ ¼ Hðxi ; yi ÞN i ðx; yÞ þ Hðxj ; yj ÞNj ðx; yÞ þ Hðxk ; yk ÞN k ðx; yÞ; ð10Þ where Hðxi ; yi Þ; Hðxj ; yj Þ; Hðxk ; yk Þ are the heads at vertexes i; j; k, respectively; N i ; Nj ; N k are linear base functions associated to vertexes i; j; k, respectively. According to the material composition, the study region X can be divided into several domains. For simplicity, assume the study region has been divided into two one material domains, zone 1 and zone 2, separated by material interface AB (Fig. 2). Each domain only contains one material, the hydraulic conductivities of zone 1 and zone 2 are different. Different from Yeh’s Galerkin model, the DDFEM does not directly solve the velocity field in the whole study region. Instead, it breaks the velocity computation problem into subproblems in the one material domains, i.e., in zone 1 and zone 2. Then, the DDFEM directly solves the subproblems by using Yeh’s model (Yeh, 1981). The velocities in both zone 1 and zone 2 are all affected by material interface, the DDFEM employs global hydraulic heads in Yeh’s model to capture the information of interfaces. Therefore, the DDFEM can ensure the accuracy of the velocities in the vicinity of interface of each subprobelm. Assume the subproblem in zone 1 need to be solved first. The V h ; h ¼ x; y computation in zone 1 is presented as follows. Consider the Darcy’s equation in zone 1:

V h ¼ K h

T

tion 3.2. Assume the unit vectors n ¼ ½nx ; ny  ; s ¼ ½ny ; nx  , the   K x K xy ; K xy ¼ K yx . For simhydraulic conductivity tensor K ¼ K yx K y   Kx 0 . plicity, assume K xy ¼ 0, i.e., K ¼ 0 Ky According Eq. (5), the refraction law can be presented in the detail form:

nx V þx þ ny V þy ¼ nx V x þ ny V y

ð7Þ

ny þ nx þ ny nx V  V ¼  V x   V y Kx Ky K þx x K þy y

@H ; @h

h ¼ x; y

ð11Þ

where K h is the hydraulic conductivity in direction h; h ¼ x; y. Multiplying both sides of (11) by linear base function N I and integrate over zone 1, the weak formulation of Darcy’s law can be obtained:

ZZ

ZZ V h NI dxdy ¼  zone1

Kh zone1

@H NI dxdy; @h

I ¼ 1; 2; . . . ; nn ; h ¼ x; y

where nn is the total number of the nodes on the grid of zone 1.

The jump form of the refraction law can be written as

 1 K y ½V x jump ¼V þx  V x ¼ n2x þ  n2y Kx "  ! #  Ky Ky K y þ 2 þ n  V  1  n V n x y y y x K x K þx K þy !1 K ½V y jump ¼V þy  V y ¼ n2y þ x n2x Ky " ! #     Kx Kx K x þ 2 þ n n  V  1  n V x y x x y K y K þy K þx

ð8Þ

3.1. Implementation of DDFEM for two dimensional steady flow problem Consider the 2D steady flow equation:

r  Kðx; yÞrH ¼ W; ðx; yÞ 2 X

ð12Þ

ð9Þ

where K is hydraulic conductivity tensor, H is hydraulic head, W is the source sink term, X is the study region.

Fig. 2. The study region of Example 1.

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

According to Yeh’s model, the velocity can be represented by linear base functions in each element Mijk in zone 1:

V h ðx; yÞ ¼ V h ðxi ; yi ÞNi ðx; yÞ þ V h ðxj ; yj ÞNj ðx; yÞ þ V h ðxk ; yk ÞNk ðx; yÞ; h ¼ x; y

ð13Þ

where V h ðxi ; yi Þ; V h ðxj ; yj Þ; V h ðxk ; yk Þ are the velocities at vertexes i; j; k, respectively. By incorporating (13) into (12), the following coefficients matrix can be obtained:

ZZ V h N I dxdy ¼

X ZZ Mijk

zone1

X ZZ

¼

Mijk

Mijk

Mijk

I ¼ 1; 2; . . . ; nn ; h ¼ x; y

½V h ðxi ; yi ÞNi þ V h ðxj ; yj ÞNj þ V h ðxk ; yk ÞNk NI dxdy

X ZZ ¼ ½ ZZ þ

V h NI dxdy;

ZZ NI Ni dxdyV h ðxi ; yi Þ þ

Mijk

NI Nj dxdyV h ðxj ; yj Þ

NI Nk dxdyV h ðxk ; yk Þ

35

We remark that: 1. Since K xy ¼ K yx , if K xy – 0, the subproblem @H Eq. (11) need to be changed to: V h ¼ K h @H  K hh0 @h0 ; h ¼ x; y; @h h0 ¼ y; x. 2. The computation of V h only relates to K h and K hh0 , so the above process is also suit for anisotropic media. 3. Fig. 2 is a sample study region, the above process is suit for any types interfaces. 4. If a study region contains more than two zones. After each subproblem is solved, the DDFEM needs to obtain the adjacent subproblems’ Dirichlet boundary conditions and then solves the next subproblem, until all the subproblems are solved (Fig. 1). 5. The subproblems contain more information are recommended to be solved priority, such as those with study region boundary condition or obtained Dirichlet boundary condition. 6.In high heterogeneous media, the FEM suffers inaccuracy in the heads computation (Forsyth, 1990; Edwards, 2006). The DDFEM is very efficient, so it can employ a finer grid to obtain better accuracy. On the other hand, the DDFEM can also employ other more efficient method to replace FEM as the head solver.

3.2. The velocity computational cost of DDFEM

ð14Þ By incorporating (10) into the right side of (12), the right side item can be obtained: ZZ 

X ZZ @H @H N I dxdy ¼  NI dxdy; I ¼ 1; 2; . . . ; nn ; h ¼ x; y Kh @h @h zone1 Mijk   ZZ X @N i @Nj @N k M ¼ K h ijk Hðxi ; yi Þ þ Hðxj ; yj Þ þ Hðxk ; yk Þ N I dxdy @h @h @h Mijk " ZZ ! ! ZZ X @N i @N j M M K h ijk NI K h ijk NI ¼ dxdy Hðxi ; yi Þ þ dxdy Hðxj ; yj Þ @h @h Mijk Mijk ! ZZ @N k M þ K h ijk NI dxdy Hðxk ; yk Þ @h Mijk Kh

ð15Þ M

where K h ijk ; h ¼ x; y is the hydraulic conductivity of Mijk in direction h. Therefore, the matrix form can be obtained:

½A½V ¼ ½F

ð16Þ

with

X Z AIJ ¼ M ijk

FI ¼ 

Mijk

N I NJ dxdy

X X ZZ Mijk h¼i;j;k

Instead of solving the velocity in the whole study region, the DDFEM solves the velocity by subproblems. Assume that the study region X contains n nodes and needs to be divided into m elements. In Yeh’s Galerkin model, the velocity problem is directly solved in X, so it requires to solve a n  n equation system. In DDFEM, assume X can be divided into c one material subdomains X1 ; X2 ; . . . ; Xc by material interfaces. Assume the numbers of nodes need to be solved in each domain are n1 ; n2 ; . . . ; nc . Since DDFEM only needs to solve each interface node for one time, so that n1 þ n2 þ . . . þ nc ¼ n. The DDFEM only needs to solve a ns  ns equation system in each subdomain Xs to obtain all the velocities in X. Therefore, the computational cost of DDFEM is much less than that of Yeh’s Galerkin model. For example, assume Yeh’s Galerkin model and DDFEM employs the mesh of Fig. 2. Yeh’s Galerkin model needs to solve a 961  961 equation system, while DDFEM only needs to solve two low order equation systems, a 496  496 equation system and a 465  465 equation system. The 31 nodes of interface AB only need to be solved in zone 1. Since the sum order of the low order equation systems is 961, so that the cost of DDFEM is much less than Yeh’s Galerkin model.

M

Mijk

K h ijk NI

! @Nh dxdy Hðxh ; yh Þ; h ¼ x; y @h

where ½A is the global coefficient matrix, which is symmetric and positive definite; ½V is the unknown velocity vector; ½F is the right hand item vector. By solving the equation system (16), the nodal Darcian velocity V h in zone 1 can be obtained. Therefore, the nodal Darcian velocities at interface AB of the zone 1 side are known. The velocities at interface AB of the zone 2 side can be obtained by equation Eq. (8), which can be regarded as the Dirichlet boundary condition of the subproblem (Eq. (11)) in zone 2. In other words, the DDFEM can bring these zone 2 side interface velocities into the left side of equation system (16) of zone 2 to obtain some known items, and move them to right side as a part of right side item. Therefore, the DDFEM does not require to solve the interface velocities again in zone 2, the Dirichlet boundary also ensures the velocities of the two sides of AB are unique and exactly satisfy the refraction law. After the boundary condition of the zone 2 is specified, the DDFEM solves Eq. (11) in zone 2 to obtain the velocities. Therefore, the velocity field in Fig. 2 can be obtained.

4. Numerical verification In this section, the DDFEM is compared with two classical methods-Yeh’s finite element model (Yeh, 1981) and Zhou’s S1 scheme (Zhou et al., 2001)-for solving Darcian velocity in heterogeneous media. The applicability of DDFEM to deal with vertical, oblique and intersected material interfaces is demonstrated in Example 1, 2, 3, respectively. Both of the isotropic and anisotropic cases are considered, so as to show the efficiency and effectiveness of DDFEM. All programs used to implement the methods described in this section are compiled in C++ and run on the same computer, using Cholesky decomposition method as the matrix solver. For convenience, some abbreviations are used: hydraulic head (H); nodal Darcian velocity in x direction (V x ); nodal Darcian velocity in y direction (V y ); Analytical solution (AS); finite element method (FEM); finite element method with fine grid mesh (FEMF); Yeh’s finite element model (Method-Yeh); Yeh’s finite element model employs the analytical head (Method-Yeh-AS); Zhou’s S1 scheme (Method-Zhou); Method-Zhou with fine grid mesh (Method-Zhou-F); the domain decomposed finite element method (DDFEM); the domain decomposed finite element method employs the analytical head(DDFEM-AS).

36

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

In this section, some conditions are specified. The head distributions of all the methods are calculated by FEM, except for MethodZhou-F (by FEM-F). The K xy ¼ K yx ¼ 0. The prescribed precision of the iteration of Method-Zhou and Method-Zhou-F is 105 . In the computation of relative error of velocity, if the analytical V x ¼ 0 and the calculated V x < 105 at a node, the relative error of V x at this node is 0%, otherwise is 100%. 4.1. Example 1 This example is presented to demonstrate the applicability of DDFEM to deal with interfaces vertical to axis. The problem equation is Eq. (9), the study region X ¼ ½0; 1  ½0; 1. The study region is divided by material interface ABðy ¼ 0:5Þ into two domains: zone 1 and zone 2 (Fig. 2). The hydraulic conductivities of zone 1 in x and y directions are K 1x and K 1y , respectively. Those in zone 2 are K 2x and K 2y , respectively. This example has analytical solution, which specified the Dirichlet boundary condition and the source/sink term. The hydraulic heads in zone 1 and 2 are:

H ¼ sinð10xyÞ  x3 y4 þ 1;

ð17Þ

in zone 1 4

ð2y þ 9Þx 2y þ 9 H ¼ sin½   x3 ð Þ þ 1; 2 20

in zone 2

The analytical Darcy velocities V x in zone 1 and 2 are:

  V x ¼ K 1x 10ycosð10xyÞ  3x2 y4 ; in zone 1 ð18Þ    4 ð2y þ 9Þ ð2y þ 9Þx 2y þ 9 cos  3x2 g; in zone 2 V x ¼ K 2x f 2 2 20 The analytical Darcy velocities V y in zone 1 and 2 are:

V y ¼ K 1 y½10xcosð10xyÞ  4x3 y3 ; (    3 ) ð2y þ 9Þx 3 2y þ 9  0:4x V y ¼ K 2 y xcos ; 2 20

in zone 1 ð19Þ in zone 2

4.1.1. Case 1: isotropic media In case 1, the hydraulic conductivity in each zone is isotropic, i.e.,K 1 x ¼ K 1 y ¼ 1; K 2 x ¼ K 2 y ¼ 10. Case 1 is solved by DDFEM, Method-Zhou and Method-Yeh. The study region is divided into 1800 (30  30  2) triangular elements by all the three methods (Fig. 2). The hydraulic heads of all the three methods are solved by FEM, which are displayed in Fig. 3. It can be seen that the

hydraulic heads vary differently in zone 1 and 2, and the heads at interface AB are continuous. After the flow field is obtained, the three methods can solve the velocity field. Fig. 4 shows the V x of analytical solution, DDFEM, Method-Zhou and Method-Yeh. The V x fields of the analytical solution, DDFEM and Method-Zhou are nearly the same. As for Method-Yeh, the calculated velocities near interface AB have larger errors and vary oscillatory, those far from AB are close to the analytical solution. In order to compare the velocities near interface AB, the V x along section x ¼ 0:5 and interface AB ðy ¼ 0:5Þ are shown in Fig. 5(a) and (b), respectively. In Fig. 5(a), it can be seen that the results of DDFEM and Method-Zhou are very close, and they are almost coincide with the analytical solution. Same as the analytical solution, the DDFEM and Method-Zhou obtain two different velocities at node ð0:5; 0:5Þ for the two sides of interface AB. Method-Yeh achieves only one velocity at node ð0:5; 0:5Þ, which is larger than those of the other two methods at the lower hydraulic conductivity side and smaller than those at the higher hydraulic conductivity side. Thus, by averaging, Method-Yeh reduces the velocities at the adjacent nodes at the lower hydraulic conductivity side and increases those at the higher hydraulic conductivity side. Moreover, same as that observed in the work of Zhou et al. (2001), the velocities at the nodes behind the adjacent nodes are also affected. Those at the lower hydraulic conductivity side are increased and those at the higher hydraulic conductivity side are decreased. In Fig. 5(b), same as the values at node ð0:5; 0:5Þ in Fig. 5(a), DDFEM and Method-Zhou achieve two different velocities for the two sides at each node of AB, respectively. The values of the higher curves of DDFEM and Method-Zhou are 10 times of the those of the lower curves, showing the velocities satisfy the refraction law. The curve of Method-Yeh is in the middle of the two curves of analytical solution, which are continuous along interface AB. Method-Yeh does not consider the refraction law, which leads to the inaccuracy at the interface nodes. In order to show the efficiency of DDFEM, the study region is divided into 7200 (60  60  2) and 16200 (90  90  2) triangular elements. The average relative errors of V x and the CPU time with the three different meshes (with 1800, 7200, 16200 elements) are represented in Table 1. It can also be seen that the average relative errors of DDFEM are very close to those of Method-Zhou, and much smaller than those of Method-Yeh. The accuracy of Method-Zhou is a little higher than DDFEM, which is because the velocities of Method-Zhou are solved in the global study region and those of

Fig. 3. The hydraulic head distribution of Example 1,case 1.

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

37

Fig. 4. The V x fields of different methods in Example 1,case 1.

Fig. 5. The V x of different methods along section x = 0.5 (a) and interface AB(b) in Example 1,case 1.

DDFEM are solved in local subdomains. Finer grid can decrease this difference. The CPU time presented in Table 1 is the sum for heads and V x . The Method-Zhou requires to incorporate the jump funcTable 1 Average relative errors of V x and CPU time of different methods in Example 1,case 1. Methods

Number of elements

Average RE of V x (%)

CPU time

DDFEM DDFEM DDFEM Method-Zhou Method-Zhou Method-Zhou Method-Yeh Method-Yeh Method-Yeh

1800 7200 16200 1800 7200 16200 1800 7200 16200

3.25 1.17 0.67 2.80 0.75 0.35 9.14 4.01 2.58

less than 1 s 21 s 242 s 6s 349 s 3635 s 1s 37 s 420 s

tion to the velocity equation system by iteration, so that it is more expensive than the other two methods. The iteration times of Method-Zhou with the three grids are the same (15 times), showing the mesh refinement can’t effectively reduce iteration times. The increase of prescribed precision of Method-Zhou can help to reduce iteration times with the sacrifice of accuracy, but can’t make its CPU time lower than that of Method-Yeh. The order of the equation system of DDFEM is much lower than the other methods, and the DDFEM does not require iteration. Therefore, the CPU time of DDFEM is less than 7% of that of Method-Zhou and about a half of that of Method-Yeh. Meanwhile, the CPU time for heads of all the three methods is the same, which is: less than 1s (1800 elements), 16s (7200 elements) and 197s (16200 elements). The CPU time of DDFEM for solving V x is about 25% of that of Method-Yeh and less than 2% of that of Method-Zhou.

38

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

In this example, zone 2 is the first solved subdomain, and zone 1 is solved by employing the Dirichlet boundary condition from zone 2. To see whether the affection of the choice of first-solved subdomain is much, the results of ”zone 1 solved first” situation and ‘‘zone 2 solved first” situations are compared. Let the first solved subdomain of DDFEM being zone 1 and then solve zone 2 with the Dirichlet boundary condition obtained from zone 1. The V x field computed with 1800 elements is presented in Fig. 6, which is very

close to Fig. 4(b). The average relative error of DDFEM in Fig. 6 is 3.31%, and that in Fig. 4(b) is 3.25%. This result indicate that the results of ‘‘zone 1 solved first” situation and ”zone 2 solved first” situations are very close but not exactly the same. Furthermore, finer grid can further decrease the affection of the choice of the first-solved subdomain. Let the number of elements being 7200 and 16,200, the average relative errors of DDFEM are 1.21% (7200 elements), 0.70% (16200 elements), which are also very close to

Fig. 6. The V x field of DDFEM with the zone 1 as the first solved subdomain in Example 1,case 1, the element number is 1800.

Fig. 7. The V x fields of different methods in Example 1, case 2.

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

the results of ‘‘zone 2 solved first” situation (Table 1). Compare to the affection of grid, the affection of the choice of the first-solved subdomain is very little, which can’t determine the accuracy of DDFEM. 4.1.2. Case 2: anisotropic media In case 2, the hydraulic conductivity in each zone is anisotropic, i.e., K 1x ¼ ð1 þ xÞð1 þ yÞ;K 1y ¼ 10ð1 þ xÞð1 þ yÞ; K 2x ¼ 10ð1 þ xÞð1 þ yÞ; K 2y ¼ 100ð1 þ xÞð1 þ yÞ. The study region is divided into 7200

39

(60  60  2) triangular elements. The hydraulic head field is solved by FEM, and the velocity fields are solved by DDFEM, Method-Zhou and Method-Yeh. The V x fields of the three methods and analytical solution are shown in Fig. 7. The subfigures of Fig. 7 have the similar shape as those of Fig. 4, but the ranges of values are different.It can be seen that the accuracy from high to low is DDFEM, Method-Zhou and Method-Yeh. The results of Method-Yeh of the both sides of the interface AB vary oscillatory, and those of the DDFEM and

Fig. 8. The V y fields of different methods in Example 1, case 2.

Fig. 9. The V x of different methods along section x = 0.5 (a) and interface AB (b) in Example 1, case 2.

Fig. 10. The study region of Example 2.

40

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

Method-Zhou are nearly the same as that of analytical solution. The V y fields are displayed in Fig. 8. Following the refraction law, the results of V y of DDFEM and Method-Zhou along AB are continuous, so the results of Method-Zhou and Method-Yeh are exactly the same.

Fig. 9(a) and (b) show the V x along section x ¼ 0:5 and interface AB ðy ¼ 0:5Þ, respectively.In Fig. 9(a), the curves of the DDFEM and Method-Zhou are very close to the analytical solution, the anisotropic media does not affect the ability of the two method of dealing with the material interface. Since the hydraulic conductivity is

Fig. 11. The hydraulic head distribution of Example 2, case 1.

Fig. 12. The V x fields (a–d) and V y fields (e–h) of different methods in Example 2, case 1.

Fig. 13. V x of different methods along section x ¼ 80 m (a) and y ¼ 80 m (b) in Example 2, case 1.

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

larger than that in case 1, the DDFEM and Method-Zhou have a larger jump at node ð0:5; 0:5Þ. Method-Yeh achieves continuous result at node ð0:5; 0:5Þ, which is in the middle of the two values of analytical solution. By averaging, the velocities of Method-Yeh at the adjacent nodes are also affected by interface AB, which have larger errors. In Fig. 9 (b), it can be clearly seen that the DDFEM and MethodZhou ensure their accuracy along the material interface, their curves are coincide to that of the analytical solution. Both of the two methods have two different values for the two sides of AB, showing that the refraction law can indeed make DDFEM and Method-Zhou achieve accurate velocities along material interface. In this case, Method-Zhou requires 353s for solving heads and V x (15 times iteration), that of Method-Yeh is 37s. The DDFEM achieves nearly the same result as Method-Zhou, while only using 17s for heads and 5s for V x , i.e, 22s CPU time. In conclusion, the results in case 1 and case 2 show that the DDFEM is much more efficient than the other two methods, while ensuring the accuracy at interface nodes. 4.2. Example 2 Since most urban centers are found on alluvial plains, this example is considered in such situation. The ability of DDFEM to deal with oblique material interfaces is tested, and compared with Method-Zhou and Method-Yeh. The problem equation is Eq. (9), the study region X ¼ ½0; 120 m  ½0; 120 m, source/sink term W ¼ 0. The hydraulic conductivities of zone 1 in x and y directions are K 1x and K 1y , respectively. Those in zone 2 are K 2x and K 2y , respectively. Those in zone 3 are K 3x and K 3y , respectively. The hydraulic heads are 1m on the right boundary and 0 m on the left boundary. The top boundary has a specified inflow flux (q), and the bottom boundary is impermeable. The study region is divided by material interfaces AB (y ¼ 20 m) and CD (x þ y ¼ 160 m) into three domains: zone 1, zone 2 and zone 3 (Fig. 10). The interface AB is vertical to the y axis, and CD does not be vertical to x or y axis. 4.2.1. Case 1: isotropic media In case 1, the hydraulic conductivity in each zone is isotropic, i.e.,K 1x ¼ K 1y ¼ K 3x ¼ K 3y ¼ 100m=d; K 2x ¼ K 2y ¼ 10m=d. The top boundary is impermeable, i.e, q ¼ 0m=d. Case 1 does not have analytical solution, so the solution of Method-Zhou-F is regarded as the standard solution. Method-Zhou-F divides the study region into 18432 (96  96  2) triangular elements. The heads of Method-Zhou-F are solved by FEM-F. The DDFEM, Method-Zhou and Method-Yeh divide the study region into 1152 (24  24  2) triangular elements. The head distribution of DDFEM is shown in Fig. 11, which is obtained by FEM. The hydraulic heads vary differently in the three different zones, the material interfaces AB and CD can be seen clearly.

41

The V x fields and V y fields are displayed in Fig. 12. It can be seen that the curve in the subfigure of Method-Zhou-F is more smoother than those of other methods. The results of DDFEM and MethodZhou are nearly the same, and very close to that of MethodZhou-F. The velocities of Method-Yeh near material interfaces have larger errors, and those far from interfaces are accurate. Since AB is vertical to axis y, the V y of all methods at interface AB are continuous. The V x of Method-Zhou-F, DDFEM, Method-Zhou and MethodYeh at the nodes of section x ¼ 80m and y ¼ 80m are presented in Fig. 13. The results show that the DDFEM achieves the same accuracy as Method-Zhou, and almost coincide with the standard solution (Method-Zhou-F). At the nodes of interfaces AB and CD, the DDFEM, Method-Zhou and Method-Zhou-F achieve two values for the two sides of interfaces, discontinuous at material interfaces. On the other hand, the solution of Method-Yeh is continuous at the two interfaces, and does not satisfy the refraction law. Moreover, the velocities of Method-Yeh at the nearest 2–3 nodes of each interface node are also affected. The nearest one connecting to the material nodes makes up for the value at interfaces node, and the second nearest one makes up for the value of the nearest one. Therefore, the curve of Method-Yeh has oscillatory wiggling variation in this figure. Fig. 14 shows the V x of different methods along material interfaces AB and CD. The solutions of DDFEM and Method-Zhou are very close to that of Method-Zhou-F, showing the high accuracy of DDFEM and Method-Zhou. Method-Yeh achieves average values of the two values of Method-Zhou-F, which are continuous but not satisfies the refraction law. The detail values of V x along interface AB are shown in Table 2 the symbol ‘‘+” means the zone 1 side and symbol ‘‘” means zone 2 side. The values of Method-Zhou1 F, DDFEM and Method-Zhou at zone 2 side are 10 of those at zone 1 side, obtained from the refraction law. The detail values of V x along interface CD are shown in Table 3, the symbol ‘‘+” means the zone 3 side and symbol ‘‘” means zone 2 side. It can be seen that only the values of Method-Yeh do not satisfy the refraction law, and the Method-Zhou and DDFEM keep high accuracy. The errors of results of Method-Zhou and DDFEM are a little higher than those in Table 2, which is because the head near CD changes more severe than those near AB. The CPU time of DDFEM in this case is less than 1s, that of Method-Yeh is 1s. Method-Zhou costs 4s for heads and V x , including 15 times iterations for V x . Method-Zhou-F costs 10119s, and iterates 14 times for V x . If DDFEM employs the same mesh as Method-Zhou-F, the CPU times for heads and V x only increase to 295s and 143s, i.e., 438s (5% of that of Method-Zhou-F). Assume K 1x ¼ K 1y ¼ K 3x ¼ K 3y ¼ K, and keep the other conditions the same as above. Fig. 15 shows the V x along section x ¼ 80 m with different K. In Fig. 15(a), Fig. 13(a), and Fig. 15(b), K equals to 10m=d, 100m=d, 1000m=d, respectively. Denote the ratio of hydraulic conductivity between the two sides of interfaces

Fig. 14. V x of different methods along interface AB(a) and CD (b) in Example 2, case 1.

42

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

Table 2 V x of different methods along AB in Example 2, case 1. x

Zhou-F+

Zhou-F-

DDFEM+

DDFEM-

Zhou+

Zhou-

Yeh

5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115

0.883354 0.882330 0.880629 0.878265 0.875259 0.871634 0.867421 0.862656 0.857381 0.851642 0.845493 0.838992 0.832200 0.825185 0.818020 0.810783 0.803564 0.796464 0.789614 0.783187 0.777436 0.772730 0.769562

0.0883354 0.0882330 0.0880629 0.0878265 0.0875259 0.0871634 0.0867421 0.0862656 0.0857381 0.0851642 0.0845493 0.0838992 0.0832200 0.0825185 0.0818020 0.0810783 0.0803564 0.0796464 0.0789614 0.0783187 0.0777436 0.077273 0.0769562

0.883497 0.882437 0.880726 0.878358 0.875344 0.871709 0.867485 0.862708 0.857421 0.851670 0.845508 0.838995 0.832191 0.825166 0.817991 0.810746 0.803520 0.796414 0.789562 0.783119 0.777408 0.772416 0.770163

0.0883497 0.0882437 0.0880726 0.0878358 0.0875344 0.0871709 0.0867485 0.0862708 0.0857421 0.0851670 0.0845508 0.0838995 0.0832191 0.0825166 0.0817991 0.0810746 0.080352 0.0796414 0.0789562 0.0783119 0.0777408 0.0772416 0.0770163

0.883789 0.882333 0.880755 0.878349 0.875344 0.871707 0.867484 0.862707 0.857420 0.851669 0.845507 0.838994 0.832191 0.825166 0.817991 0.810747 0.803522 0.796418 0.789566 0.783127 0.777402 0.772434 0.770314

0.0883789 0.0882333 0.0880755 0.0878349 0.0875344 0.0871707 0.0867484 0.0862707 0.0857420 0.0851669 0.0845507 0.0838994 0.0832191 0.0825166 0.0817991 0.0810747 0.0803522 0.0796418 0.0789566 0.0783127 0.0777402 0.0772434 0.0770314

0.477941 0.484160 0.484865 0.483004 0.481436 0.479431 0.477107 0.474480 0.471572 0.468409 0.465020 0.461438 0.457696 0.453833 0.449887 0.445902 0.441929 0.438021 0.434248 0.430781 0.427164 0.425810 0.430649

Table 3 V x of different methods along CD in Example 2, case 1. x

Zhou-F+

Zhou-F-

DDFEM+

DDFEM-

Zhou+

Zhou-

Yeh

45 50 55 60 65 70 75 80 85 90 95 100 105 110 115

0.383499 0.314764 0.279260 0.255454 0.237411 0.222648 0.209869 0.198282 0.187335 0.176586 0.165613 0.153934 0.140851 0.125055 0.103055

0.173141 0.140917 0.123930 0.112450 0.103771 0.096762 0.090831 0.085614 0.080854 0.076342 0.071874 0.067209 0.062004 0.055632 0.046504

0.421947 0.320561 0.283942 0.259072 0.240229 0.225090 0.212072 0.200355 0.189357 0.178626 0.167723 0.156293 0.143179 0.129747 0.102615

0.196444 0.145393 0.127710 0.115207 0.106037 0.098744 0.092664 0.087381 0.082623 0.078175 0.073822 0.069412 0.064377 0.059528 0.049426

0.394480 0.318609 0.280882 0.257161 0.238702 0.223789 0.210895 0.199239 0.188253 0.177494 0.166529 0.154956 0.141722 0.127624 0.099883

0.177533 0.142302 0.124826 0.113319 0.104518 0.097452 0.091493 0.086269 0.081524 0.077047 0.072632 0.068085 0.062951 0.057199 0.048552

0.291284 0.230592 0.202691 0.185276 0.171606 0.160620 0.151194 0.142754 0.134888 0.127271 0.119578 0.111529 0.102303 0.092478 0.074737

Fig. 15. V x of different methods with different K along section x ¼ 80 m in Example 2, case 1.

as h ¼ KK2 . Therefore, the Fig. 15 (a) presents the results in a study region with h ¼ 1, the solutions of DDFEM, Method-Zhou and Method-Yeh are continuous at all the interfaces nodes. Compare the results of Fig. 15(a) and (b), the jump between different sides of interface increases with the increase of h. The errors of Method-Yeh in Fig. 15(b) are larger than those in Fig. 15 and Fig. 13(a), showing larger h leads to larger errors at interface nodes.

Meanwhile, Method-Zhou requires to iterate 0, 15, 18 times for the result in Fig. 15(a), Fig. 13(a), and Fig. 15 (b), showing larger h leads to more iteration times. Therefore, the Method-Zhou needs more computational cost with h increases. Furthermore, in all the three different conditions, the CPU time of DDFEM is the same, which is much less than that of Method-Zhou and Method-Yeh. The DDFEM achieves nearly the same accuracy as Method-Zhou in all

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49 Table 4 The left side values of mass balance Eq. (20) of DDFEM and Method-Zhou in Example 2, case 1. Methods

Local region X10

Local region X20

Local region X30

DDFEM (gird 1) Method-Zhou (grid 1) DDFEM (gird 2) Method-Zhou (grid 2)

0.00062 0.00138 0.00033 0.00069

0.00070 0.00132 0.00036 0.00066

0.00361 0.00188 0.00108 0.00060

the three different conditions, showing that DDFEM can be more efficient than Method-Zhou with a larger h. To see whether the DDFEM is mass balance (Yeh, 1981; Chen and Hou, 2003), the DDFEM and Method-Zhou divides the study region into 4608 triangular elements (48  48  2, grid 1), and 18432 triangular elements (96  96  2, grid 2). The hydraulic conductivity is isotropic, i.e.,K 1x ¼ K 1y ¼ K 3x ¼ K 3y ¼ 100m=d; K 2x ¼ K 2y ¼ 10m=d. After solving for the velocities, the following equation is considered:

Z

Z

@ X0

V h  nds ¼

X0

WdX0 ; 8X0 2 X

ð20Þ

where X0 is an arbitrary local region belong to X; n the unit outer normal vector; V h is the velocity in the h direction, where h ¼ x; y; and W is the source/sink term. Let us consider Eq. (20) in three local regions. All local regions are 10 m  10 m rectangular domains, with their origin points at coordinates of ð10 m; 10 mÞ (local region X10 ), ð10 m; 20 mÞ (local region X20 ), ð80 m; 80 mÞ (local region X30 ). In this case, the lefthand side values of Eq. (20) are computed by DDFEM and Method-Zhou. Since the source/sink term W ¼ 0, the right-hand side of Eq, (20) is equal to 0. So the analytical solution of lefthand side value is 0. The left-hand side values computed by DDFEM and Method-Zhou are shown in Table 4. It can be seen that the absolute values of the left-hand side values of the DDFEM are very close to 0, showing that the mass conservation of the DDFEM is satisfied to an acceptable degree. The values of DDFEM with different local regions are close to those of Method-Zhou, indicates the abilities of mass balance of the two methods are close. When the grid is finer, the differences of the values between DDFEM and MethodZhou are much smaller. The results of DDFEM and Method-Zhou with grid 2 are much smaller than those with grid 1, showing that finer grid can improve the ability of mass balance of the DDFEM and Method-Zhou. 4.2.2. Case 2: anisotropic media In case 2, the DDFEM, Method-Zhou and Method-Yeh are compared to deal with anisotropic media., i.e., K 1x ¼ K 3x ¼ 100m=d; K 1y ¼ K 3y ¼ P  100m=d; K 2x ¼ 10; K 2y ¼ P  10m=d, and P ¼ 3; 5; 10. The top boundary has a specified inflow flux q ¼ 2m=d.

43

Case 2 does not have analytical solution, so the solution of Method-Zhou-F is regarded as the standard solution. Method-Zhou-F divides the study region into 18432 (96  96  2) triangular elements. The heads of Method-Zhou-F are solved by FEM-F. The DDFEM, Method-Zhou and Method-Yeh divide the study region into 1152 (24  24  2) triangular elements. The head fields of Method-Zhou-F are solved by FEM-F, and those of the other methods are solved by FEM. The head fields of DDFEM with different P are represented in Fig. 16. It can be seen that the inflow flux of the top boundary has larger influence when P increases, which makes the heads vary more severer in the y direction than x direction. The V x fields of the four methods with different P along section x ¼ 80 m and y ¼ 80 m are displayed in Fig. 17. The curves of DDFEM and Method-Zhou are coincide with that of MethodZhou-F, showing the anisotropic media does not affect the accuracy of DDFEM and Method-Zhou. The Method-Yeh has error at interface nodes, its result not close to either interface side result of Method-Zhou-F. The V y fields of the four methods with different P along section x ¼ 80 m and y ¼ 80 m are displayed in Fig. 18. Since interface AB is vertical to y axis, the V y values of all methods at the nodes of AB are continuous. The DDFEM and Method-Zhou keep high accuracy at the nodes of CD, but Method-Yeh suffers inaccuracy. The results of Method-Yeh, Method-Zhou and Method-Zhou-F vary oscillatory at the 2–3 nearby nodes of the CD. The oscillatory amplitude of the three methods from high to low is Method-Yeh, Method-Zhou and Method-Zhou-F, showing finer grid can reduce this oscillation. Different from the other methods, the results of DDFEM at the vicinity nodes of interface CD are still close to that of Method-Zhou-F (the DDFEM nearby nodes of CD are not belong to the oscillatory area of MethodZhou-F), showing the accuracy of DDFEM. The V x and V y with different P along interface AB and CD are shown in Fig. 19 and Fig. 20, respectively. In all situations DDFEM and Method-Zhou achieve almost the same results as that of Method-Zhou-F. Method-Yeh has errors at interface nodes, except in Fig. 20(a), (c), (e) (V y vertical to AB). The jumps of V x and V y of DDFEM, Method-Zhou and Method-Zhou-F at CD nodes become smaller with P increases. The V x jumps of DDFEM, Method-Zhou and Method-Zhou-F at AB nodes become larger with P increases, except for those near node ð80; 20Þ (V x close to 0). The V y values in Fig. 20(a), (c), (e) become larger with P increase. Therefore, the jumps change differently in AB; CD, which is because the increase of P makes variation of heads become larger at AB nodes in the y direction (Fig. 16). The results in Fig. 19 and Fig. 20 also show that the values of DDFEM and Method-Zhou at interface nodes satisfy the fraction law, the accuracy of the two methods is not affected by anisotropic media. When P ¼ 3; 5; 10, the CPU time of DDFEM for solving heads and V x is: less than 1s, less than 1s, 1s, which are same as those of Method-Yeh. Meanwhile, the Method-Zhou costs 4s, 4s, 4s for 15,

Fig. 16. The head fields of different P in Example 2, case 2.

44

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

Fig. 17. The V x of different methods with different P along section x ¼ 80m (a,c,e) and y ¼ 80 m (b,d,f) in Example 2, case 2.

Fig. 18. The V y of different methods with different P at section x ¼ 80 m (a,c,e) and y ¼ 80 m (b,d,f) in Example 2, case 2.

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

16, 16 times iteration with P ¼ 3; 5; 10, respectively. If the DDFEM employs the same division of the study region of Method-Zhou-F. When P ¼ 3; 5; 10, the DDFEM only needs 293s,291s,292s for solving heads and 145s, 146s, 145s for solving V x , i.e., 438s, 437s, 437s for solving heads and V x . The CPU time of Method-Zhou-F for heads and V x is much more than that of DDFEM. Method-Zhou-F costs 9786s, 9828s, 10426s for 14, 14, 15 times iteration with P ¼ 3; 5; 10, respectively. Since the program of all the methods are run in the same condition, the above CPU time shows the DDFEM is more efficient than Method-Zhou. 4.3. Example 3 In this example, the applicability of DDFEM to deal with intersected material interfaces is demonstrated. The problem equation is Eq. (9), the study region X ¼ ½0; 1  ½0; 1. The study region is divided by material interface ABðy ¼ 0:5Þ and CDðx ¼ 0:5Þ into four zones (Fig. 21). This example is considered in isotropic media, the hydraulic conductivities of zone 1,2,3,4 are 1,10,10,100, respectively. This example has analytical solution, which specified the Dirichlet boundary condition and the source/sink term. The analytical heads of zone 1 and zone 3 are the same as those of zone 1 and zone 2 of Example 1, respectively. The analytical heads of zone 2 and zone 4 are:

   3 ð2x þ 9Þy 2x þ 9  y4 þ 1; in zone 2 ð21Þ 2 20    3 ð2x þ 9Þð2y þ 9Þ 2x þ 9 2y þ 9 4  Þ þ 1; in zone 4 ð H ¼ sin 40 20 20

H ¼ sin

45

This example is from the work of Cainelli et al. (2012), the detail velocity and source/sink term can refer the work of Cainelli. The DDFEM, the DDFEM-AS, the Method-Yeh and the MethodYeh-AS are employed to solve this example. The four methods divide each edge of the study region into 120 equal parts, i.e., 28800 (120  120  2) triangular elements. The head fields of DDFEM and Method-Yeh are solved by FEM, which is shown in Fig. 22. The heads used by DDFEM-AS and the Method-Yeh-AS are analytical solution, so as to remove the affection of head field in the velocity computations. Different from that of Example 1 and Example 2, the DDFEM does not solve the subproblems subsequently. The subproblems of zone 1 and zone 4 are mutually independent, so that these two subproblems are solved at the same time, before those of zone 2 and zone 3. Then, the subproblems of zone 2 and zone 3 are solved by using the Dirichlet boundary condition obtained by refraction law. Therefore, the DDFEM values of V x and V y at intersection node solved by subproblems of zone 1 and zone 4 can be compared. The V x and V y at intersection node ð0:5; 0:5Þ are presented in Table 5. The values of analytical solution, DDFEM and DDFEM-AS are very close, showing the accuracy of DDFEM. The V x has smaller values in zone 1,2 and larger values in zone 3, 4, the V y has smaller values in zone 1,3 and larger values in zone 2, 4. Different from the analytical solution, the smaller (larger) values of DDFEM or DDFEM-AS are very close, but not exactly the same. This is because the values of the DDFEM and DDFEM-AS of zone 1, 4 are solved by subproblem, while those of zone 2,3 are obtained by refraction law. The values of zone 2 or zone 3 can be obtained from zone 1 or zone 4, so the values of zone 2,3 in Table 5 are the average value of the

Fig. 19. The V x of different methods with different P at interface AB (a,c,e) and interface CD (b,d,f) in Example 2, case 2.

46

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

Fig. 20. The V y of different methods with different P at interface AB (a,c,e) and interface CD (b,d,f) in Example 2, case 2.

Fig. 21. The study region of Example 3.

47

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

values from zone 1 and zone 4. For illustration, the DDFEM V x of zone 2 can be obtained by the value of zone 1 and the refraction law of CD(3.90248) or obtained by the value of zone 4 and the refraction law of AB (3.87762), so the V x of zone 2 is 3.89005. It can be seen that the values obtained from zone 1 and zone 4 are not equal but very close, the difference is only about 0:6%. And finer grid can further reduce this difference. This result indicates that the choice of first solved subproblem can not affect the accuracy of DDFEM much. It can also be seen that the results DDFEM-AS (Method-Yeh-AS) are more accurate than those DDFEM (MethodYeh), showing the errors are main caused by heads. Compare to the affection of heads, the affection of the choice of the first-

solved is very little. As for Method-Yeh and the Method-Yeh-AS, the values of V x (V y ) for the four different zones are identical, which are not accurate. The V x and V y of the four methods along interface AB and CD are presented in Fig. 23 and Fig. 24, respectively. Fig. 23(a) and Fig. 24 (b) show the velocities of DDFEM and DDFEM-AS parallel to interface are discontinuous, which have two different values. Meanwhile, the velocities of Method-Yeh and Method-Yeh-AS are continuous, but not accurate. Fig. 23(b) and Fig. 24(a) show the velocities of DDFEM, DDFEM-AS vertical to the interface are continuous except at node ð0:5; 0:5Þ, which is because node ð0:5; 0:5Þ is the intersection of interface AB and CD. In Fig. 23 and Fig. 24, it also

Fig. 22. The hydraulic head distribution of Example 3.

Table 5 The values of V x and V y at intersection node ð0:5; 0:5Þ in Example 3. Zone

Direction

AS

DDFEM

DDFEM-AS

Method-Yeh

Method-Yeh-AS

1 2 3 4 1 2 3 4

x x x x y y y y

4.05259 4.05259 40.5259 40.5259 4.06822 40.6822 4.06822 40.6822

3.90248 3.89005 38.9005 38.7762 3.91641 39.0469 3.90469 38.9296

4.08224 4.06594 40.6594 40.4964 4.09720 40.8127 4.08127 40.6533

21.3105 21.3105 21.3105 21.3105 21.3909 21.3909 21.3909 21.3909

22.2000 22.2000 22.2000 22.2000 22.2849 22.2849 22.2849 22.2849

Fig. 23. The V x of different methods along interfaces in Example 3.

48

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49

Fig. 24. The V y of different methods along interfaces in Example 3.

can be seen that the results of DDFEM are very close to the analytical solution, but has small errors near node ð0:5; 0:5Þ. Meanwhile, the curve of DDFEM-AS is almost coincide with that of the analytical solution, even near node ð0:5; 0:5Þ. The Method-Yeh-AS achieves better results than Method-Yeh near node ð0:5; 0:5Þ but still not accurate. This phenomenon also indicates that the errors of DDFEM are main caused by heads, the way of DDFEM for solving velocity does not introduce much errors. The above results also show that the DDFEM can deal with the intersected material interfaces well. The CPU time of DDFEM for solving head and V x is only 1144s, which is about a half of that of Method-Yeh (2223s). Since the head computation takes up the same 1070s for the two methods, so that the CPU time for V x of DDFEM (74s) is much less than 7% of Method-Yeh (1153s). This result shows that the DDFEM is much more efficient than Method-Yeh, while achieving better results at interfaces nodes. 5. Conclusion This paper proposes a domain decomposed finite element method (Fig. 1) for the simulation of Darcian velocity field in heterogeneous porous media with material interfaces, which is modified from Yeh’s Galerkin model. Same as Yeh’s model, this method employs finite element method to solve the global hydraulic heads distribution before velocity computation, so as to capture the global information. Instead of solving the velocity problem in the whole study region, this method employs domain decomposition technique to divide the velocity problem into subproblems in serval one material domains. Therefore, this method breaks the high order equation system of the whole study region down to low order equation systems, which can save much computational cost. In the numerical verification section, the DDFEM is compared with Method-Yeh and Method-Zhou in three different examples, including vertical, oblique or intersected material interfaces situations and isotropic or anisotropic cases. The DDFEM, Method-Zhou and Method-Yeh obtain almost identical solutions at the nodes far from the material interfaces. At interface nodes, the DDFEM achieves a nearly same accuracy as that of Method-Zhou, and more accurate than Method-Yeh. The affection of the choice of the firstsolved subproblem to DDFEM is very little, which is much smaller than the affection of grid or head. Furthermore, when the hydraulic conductivity ratio between the two sides of an interface increases, the DDFEM keeps the same efficiency, which is superior than Method-Zhou. Moreover, the mass conservation of the DDFEM is satisfied to an acceptable degree, the ability of mass balance of DDFEM is close to that of Method-Zhou. On the other hand, the DDFEM solves the velocities by subproblems and does not require any iteration, so it costs much less than Method-Yeh and Method-

Zhou. The numerical results show that the CPU time of DDFEM is less than 10% CPU time of that of Method-Zhou and about a half of Method-Yeh. Consequently, the DDFEM can achieve accurate Darcian velocities in isotropic or anisotropic heterogeneous porous media with vertical, oblique or intersected material interfaces, while saving much computational cost. In this work, the DDFEM is limited by FEM in head accuracy and mesh subdivision. In the future, we plan to combine the DDFEM frame work with multiscale finite element methods so as to further improve the effectiveness in solving velocity in high heterogeneous media. Furthermore, we also wish to extend this work to 3D cases and consider more complex conditions. Acknowledgements This study is financially supported by the National Natural Science Foundation of China Xinjiang Project (No. U1503282), the National Natural Science Foundation of China (No. 41702243), the Special Funding of China Postdoctoral Science Foundation (No. 2017T100354), the China Postdoctoral Science Foundation (No. 2016M591826) and the Natural Science Foundation of Jiangsu Province (Grant No. BK20160880). References Ahmed, R., Edwards, M.G., Lamine, S., et al., 2015. Control-volume distributed multipoint flux approximation coupled with a lower-dimensional fracture model. J. Comp. Phys. 284, 462–489. Batu, V., 1984. A finite element dual mesh method to calculate Nodal Darcy velocities in nonhomogeneous and anisotropic aquifers. Water Resour. Res. 20 (11), 1705–1717. Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier Publishing Company, New York. Bear, J., 1979. Hydraulics of Groundwater. McGraw-Hill, New York. Bellin, A., Salandin, P., Rinaldo, A., 1992. Simulation of dispersion in heterogeneous porous formations: statistics, first-order theories, convergence of computations. Water Resour. Res. 28 (9), 2211–2227. Boufadel, M.C., Suidan, M.T., Venosa, A.D., 1999. A numerical model for density-andviscosity-dependent flows in two-dimensional variably saturated porous media. J. Contam. Hydrol. 37 (1), 1–20. Cainelli, O., Bellin, A., Putti, M., 2012. On the accuracy of classic numerical schemes for modeling flow in saturated heterogeneous formations. Adv. Water Resour. 47, 43–55. Chavent, G., Roberts, J.E., 1991. A unified physical presentation of mixed, mixedhybrid finite elements and standard finite difference approximations for the determination of velocities in waterflow problems. Adv. Water Resour. 14 (6), 329–348. Chen, Z., Hou, T., 2003. A mixed multiscale finite element method for elliptic problems with oscillating coefficients. Math. Comput. 72 (242), 541–576. Chin, D.A., Wang, T., 1992. An investigation of the validity of first-order stochastic dispersion theories in isotropic porous media. Water Resour. Res. 28 (6), 1531– 1542. D’Angelo, C., Scotti, A., 2012. A mixed finite element method for Darcy flow in fractured porous media with non-matching grids. Esaim-Math. Model Numer. Anal. 46 (2), 465–489. Edwards, M.G., 2006. Higher-resolution hyperbolic-coupled-elliptic flux-continuous CVD schemes on structured and unstructured grids in 2-D. Int. J. Numer. Methods Fluids 51, 1059–1077.

Y. Xie et al. / Journal of Hydrology 554 (2017) 32–49 Ervin, V.J., 2013. Approximation of axisymmetric Darcy flow using mixed finite element methods. Siam J Numer. Anal 51 (3), 1421–1442. Forsyth, P., 1990. A control-volume finite element method for local mesh refinement in thermal reservoir simulation. SPERE, 561. Hu, M., Wang, Y., Rutqvist, J., 2015. On continuous and discontinuous approaches for modeling groundwater flow in heterogeneous media using the Numerical Manifold Method: model development and comparison. Adv. Water Resour. 80, 17–29. Hoteit, H., Firoozabadi, A., 2008. An efficient numerical model for incompressible two-phase flow in fractured media. Adv. Water Resour. 31 (6), 891–905. Jang, W., Aral, M.M., 2007. Density-driven transport of volatile organic compounds and its impact on contaminated groundwater plume evolution. Transp. Porous Media 67 (3), 353–374. Martin, V., Jaffr, J., Roberts, J.E., 2005. Modeling fractures and barriers as interfaces for flow in porous media. SIAM J. Sci. Comput. 26 (5), 1667–1691. Meissner, U., 1973. A mixed finite element model for use in potential flow problems. Int. J. Numer. Methods Eng. 6 (4), 467–473. Neuman, S.P., Orr, S., Levin, O., Paleologos, E., 1992. Theory and high-resolution finite element analysis of 2d and 3d effective permeability in strongly heterogeneous porous media. In: Russel, T.F. (Ed.), Mathematical Modeling in Water Resources, vol. 2. Elsevier, New York, pp. 118–136.

49

Park, C.H., Aral, M.M., 2007. Sensitivity of the solution of the Elder problem to density, velocity and numerical perturbations. J. Contam. Hydrol. 92 (1), 33–49. Segol, G., Pinder, G.F., Gray, W.G., 1975. A Galerkin-finite element technique for calculating the transient position of the saltwater front. Water Resour. Res. 11 (2), 343–347. Srivastava, R., Brusseau, M.L., 1995. Darcy velocity computations in the finite element method for multidimensional randomly heterogeneous porous media. Adv. Water Res. 18 (4), 191–201. Tompson, A.F.B., Gelhar, L.W., 1990. Numerical simulation of solute transport in three dimensional, randomly heterogeneous porous media. Water Resour. Res. 26 (10), 2541–2562. Xie, Y.F., Wu, J.C., Xue, Y.Q., Xie, C.H., 2016. Combination of multiscale finiteelement method and Yeh’s finite-element model for solving nodal darcian velocities and fluxes in porous media. J. Hydrol. Eng. 04016048. Yeh, G.T., 1981. On the computation of Darcian velocity and mass balance in the finite element modeling of groundwater flow. Water Resour. Res. 17 (5), 1529– 1534. Zhang, Z.H., Xue, Y.Q., Wu, J.C., 1994. A cubic-spline technique to calculate nodal Darcian velocities in aquifers. Water Resour. Res. 30 (4), 975–981. Zhou, Q.L., Bensabat, J., Bear, J., 2001. Accurate calculation of specific discharge in heterogeneous porous media. Water Resour. Res. 37 (12), 3057–3069.