(W, LIBRARY
5:2 {.2 32 «in: 4/ A/ Michigan State
V University

This is to certify that the
thesis entitled

AN IMPROVED FINITE DIFFERENCE METHOD FOR
SOLVING COMPLEX GROUNDWATER FLOW PROBLEMS
IN GENERAL 2-D ANISOTROPIC MEDIA

presented by
SOHEIL AFSHARI

has been accepted towards fulfillment
of the requirements for the

Master of degree in Department of Civil and
Science Environmental Eng‘neering

1 Looks»

. )ajor )Profejr’s Signature /

 

 

12/
I

Date

2003

 

MSU is an Affinnative Action/Equal Opportunity Institution

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c:/CIRC/DateDue.p65«p.15

 

AN IMPROVED FINITE DIFFERENCE METHOD FOR SOLVING COMPLEX
GROUNDWATER FLOW PROBLEMS IN GENERAL 2-D ANISOTROPIC MEDIA

By

Soheil Afshari

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE
Department of Civil and Environmental Engineering

2003

ABSTRACT

AN IMPROVED FINITE DIFFERENCE METHOD FOR SOLVING COMPLEX .
GROUNDWATER FLOW PROBLEMS IN GENERAL 2-D ANISOTROPIC MEDIA

By

Soheil Afshari

In this thesis, we present an improved finite difference (FD) method for solving complex
2-D groundwater flow problems in anisotropic media. In particular, we address one of the
major difficulties in solving the tensorial form of the groundwater flow equation.
Traditional FD techniques, when applied to the flow equation containing the complete
conductivity tensor, often lead to physically unrealistic results. Traditional FD methods
are particularly problematic when the anisotropy is strong and its major orientation
deviates significantly fiom the rectilinear coordinate system. Our new approach
eliminates the problem by approximating the tensorial “diffusion” terms in a rotated
coordinate system locally aligned with the major direction of anisotropy. It then
expresses and interpolates the non-nodal heads in the resulting numerical expression in
terms of global nodal heads. The result is a significantly improved scheme that is more
accurate and robust than the traditional finite difference scheme. In addition, the new
approach is monotonic and always leads to a physically meaningful solution. We
demonstrate the advantage of the improved FD approach with a number of groundwater
flow examples and a systematic comparison with exact solutions and solutions obtained

from traditional FD methods.

ACKNOWLEDGEMENTS

I gratefully acknowledge the contribution of those who have made this thesis possible.
My deepest gratitude goes to Dr. Shu-Guang Li, Associate Professor, Department of
Civil and Environmental Engineering, Michigan State University, my thesis director
and my advisor during my course work. I appreciate his fatherly advice, patience, and

continual encouragement throughout my entire graduate program.

I also wish to extend my appreciation to Dr. Roger B. Wallace, Associate Professor,
Department of Civil and Environmental Engineering; Dr. Susan Masten, Professor,
Department of Civil and Environmental Engineering, Michigan State University; and
Dr. David C Wiggert, Professor, Department of Civil and Environmental Engineering,
Michigan State University, for serving in my thesis committee and for valuable

comments in this study.

Also, I wish to thank Dr. H.S. Liao for implementing the Fortran code for improved

FD method in Interactive Groundwater Model (IGW); and Dr. Qun Liu for his

patience and restless help during my simulation study.

Research funding was obtained from National Science Foundation (NSF).

Finally, I want to thank my family, especially my parents and my uncle

Dr. Sirous Haji-Djafari for their love, care, support and understanding.

iii

TABLE OF CONTENTS

1 INTRODUCTION ..................................................................... 1
1.1 Geological Background ................................................................... 1
1.1.1 Sloping aquifers .......................................................................... 1
1.1.2 Fractured rocks ........................................................................... 3
1.1.3 Fluvial channel deposits ............................................................. 4

1.2 Darcy’s Law In General Anisotropic Media ................................... 5
1.3 Governing Equation For 2-D Anisotropic Flow ............................. 8

2 PREVIOUS RESEARCH ....................................................... 10
2.1 Traditional FD approach ................................................................ 11

2. 1. 1 Transient flow ........................................................................... 1 1
2.1.2 Steady state flow ...................................................................... 15

2.2 The MODFLOW approach ............................................................. 17

2. 2. 1 Transient flow ........................................................................... 1 7
2.2.2 Steady state flow ...................................................................... 18

3 IMPROVED FD APPROACH ................................................ 19
3.1 Transient flow ................................................................................ 19
3.2 Steady state flow ........................................................................... 23

4 ILLUSTRATIVE EXAMPLES ................................................. 25
4.1 A Simple Example .......................................................................... 25
4.2 A More Realistic Example ............................................................. 37

5 CONCLUSION ....................................................................... 43
6 REFERENCES ...................................................................... 45

iv

List of Figures

Figure 1. A sloping aquifer system [modified from M. P. Anderson and W. W.

Woessner, 1992]. ................................................................................................... 2
Figure 2. Fractured carbonate rock [modified from Freeze and Cherry, 1979] ............ 3
Figure 3. Fluvial channel deposit [modified from Freeze and Cherry, 1979]. ............. 5

Figure 4. (a) a global coordinate system and (b) a local coordinate system aligned

with the bedding (the dashed lines represents bedding orientation) [modified

from Wang, HF and Anderson, MP, 1982] ........................................................ 7
Figure 5. Finite difference grid showing index numbering convention ...................... 12
Figure 6. Nodal and non-nodal head values ................................................................ 20
Figure 7. Definition sketch for the simple example .................................................... 26
Figure 8. Exact Solution for the simple example ........................................................ 29

Figure 9. Numerical solution using improved FD method for the simple example. .. 30
Figure 10. Numerical solution using traditional FD method for the simple example.31
Figure 11. Numerical solution using the MODFLOW approach for the simple
example. ............................................................................................................... 32
Figure 12. Comparison of the weighting coefficients for (a) the improved FD method,
(b) the traditional FD method, and(c) the MODFLOW method .......................... 35

Figure 13. Conceptual sketch with assigned physical and numerical parameters for

the realistic example. ........................................................................................... 37
Figure 14. Results of realistic example applying the improved FD approach . .......... 40
Figure 15. Results of realistic example applying traditional FD approach ................. 41

List of Tables

Table 1. Input physical and numerical parameters for the simple example ................ 28
Table 2 - Comparison of the exact predicted head and velocity at three representative

points (Figure 7) with different numerical methods. ........................................... 33
Table 3. Physical parameters for three conceptual layers in the realistic example ..... 38
Table 4. Input numerical parameters for three conceptual layers in the realistic

example. ............................................................................................................... 39

vi

1 Introduction

This thesis aims at developing, documenting, and demonstrating an improved finite-
difference (FD) method for solving complex 2-D groundwater flow problems in
anisotropic media. In particular, the thesis addresses one of the major difficulties in
solving the groundwater flow equation when the aquifer is strongly anisotropic and the
major orientation of anisotropy deviates significantly from the rectilinear coordinate

system.

1.1 Geological Background

In some specific geological formations like sloping aquifers, fluvial channel deposits,
glacial deposits, or secondary permeability terrains and fractured rocks, anisotropy
expresses with orientation. In many cases, this orientation varies from point to point and
from layer to layer. Under such situations, traditional FD techniques become problematic

and lead to non-meaningfirl or physically unrealistic results.

1. 1. 1 Sloping aquifers
One of the common examples of anisotropic formations with a variable orientation is a
sloping aquifer system. Figure 1 represents a system of dipping hydrostratigraphic units

hydraulically connected with a shallow unconfined aquifer.

 

  
   
     

   

 

Overlying m

Formations

 

Confining
unit

 

 

 

 

 

 

Figure 1. A sloping aquifer system [modifiedfrom M. P. Anderson and W. W. Woessner,

1992].

While the water in the unconfined unconsolidated aquifer tends to move in the horizontal
direction, it has a natural tendency to flow at an angle along the dipping layers in the
confined aquifer system. In this case, the preferential flow is oriented in two different
directions and, obviously, it is impossible to align both simultaneously with the overall

model grid axes.

1. 1.2 Fractured rocks

Another example of an anisotropic formation with a variable orientation is a fiactured
rock system. Many rock formations have appreciable secondary permeability as a result
of fractures or openings along bedding planes. Secondary openings in rock are caused by
changes in the stress conditions and may be enlarged as a result of circulating
groundwater (e.g., in limestone or dolomite dissolution). In the case of folded rocks, the
zone of fracture concentration and enlargement are commonly associated with the crest of
anticlines and to a lesser extent with synclinal troughs [Freeze and Cherry, 1979] (Figure
2). In the fractured zone, flow tends to move much easier in the direction of natural
fracture and it has the least infiltration perpendicular to the fi'acture. In many cases, the
orientation of the fractures may vary from location to location and fiom formation to
formation. It may also likely be different from the orientation in the overlying

unconsolidated aquifers since it is created by a different mechanism.

K2 1

 

 

Figure 2. Fractured carbonate rock [modified from Freeze and Cherry, 1 979].

1. 1.3 Fluvial channel deposits

Fluvial deposits are additional example that ofien exhibits anisotropy with variable
orientation. The fluvial deposits are the materials laid down by deposition in river
channels or on flood plains. Fluvial materials occur in nearly all regions. Figure 3
illustrates the morphology and variations in deposits formed by meandering rivers.
Because of the way they are deposited and the ever-changing depositional velocities,
river deposits have characteristic textural variability that causes much heterogeneity in
the distribution of hydraulic properties [Freeze and Cherry, 197 9]. These heterogeneous
deposits are frequently stratified and oriented along the streams. When the average
properties of large volumes are considered, the heterogeneous character of fluvial
deposits imparts a strong anisotropy to the system and creates preferential flow parallel to
the stream. Obviously, the major orientation of anisotropy varies as the channel winds

around and shifis in its position [Freeze and Cherry, 1979].

  
  
   

Point bar

Channel Deposit

  

Flood basin

Figure 3. Fluvial channel deposit [modified from Freeze and Cherry, 1979].

1.2 Darcy’s Law in General Anisotropic Media
Under these situations, Darcy’s law applies only locally along the ‘principal’ directions
or the preferential flow direction and the direction normal to it. If the local principal

coordinates are denoted as x', 2' (see Figure 4), Darcy’s equation can be written as:

6CD
q x r - 'K x I 3):, (1)
6(1)
qz. "1932—, (2)
where,
qx ,,qz , the specific discharge x’ and 2' respectively [LT'l]

x ,z the local coordinate system [L]
(D the hydraulic head [L]

K x , , and K 2 , the principal hydraulic conductivities [LT"]

However, in a system wide finite—difference modeling analysis, a global rectilinear
coordinate is required. In this case, the following generalized tensorial form of Darcy’s

law must be used (Eqs. 3 and 4):

qx = _Kxx 23:2“sz g: (3)
‘12 =— 21 $492 55:9 (4)
where,

x andz the spatial coordinates in the global rectilinear system[L]
qx anqu the specific discharges along x andz respectively in the

global coordinate system [LT'I]

K ,K ,K , and K components of the hydraulic conductivity tensor which
xx xz zx 22

represent an anisotropic medium [LT'l].

 

Preferential direction of flow

 

 

 

 

   

 

Figure 4. (a) a global coordinate system and (b) a local coordinate system aligned with
the bedding (the dashed lines represents bedding orientation) [modified from Wang, HF.

and Anderson, MR, 1982].

The cross-conductivity K xz represents the x-seepage flow induced by a unit hydraulic

gradient in the z direction. For instance, the term K x2 2% in Eq. (3) reflects the fact that

a head gradient in one direction (e.g., the z direction) will induce some flow along the
beds, and hence there will be a component of flow in the other direction (e.g., the x
direction). Note that the principal conductivities and the orientation are the fundamental
properties that characterize an aquifer formation and the conductivity tensor is a derived
property and can be related to the principal conductivities and the orientation. In the case
of two dimensional (2-D) flow in the xz plane, this relationship is given by the following

equations:

K K 1 K , 0
xx xz = R — x R (5)
K 0 K ,
zx 22 z
with
_, cos 9 - sin 6
R = . (6)
51nd cos6i
where
R the rotation matrix
R" the inverse of R
6? the angle of anisotropy, which is measured from the x axes of global direction

[degree] in the counterclockwise direction.

1.3 Governing Equation For 2-D Anisotropic Flow

To systematically analyze groundwater flow using Eq. (3), and (4) we must first know the

hydraulic head distribution. An equation for determining the hydraulic head can be

obtained from the mass conservation equation:

S Q—3q;x_+.a_qL
S at 6x

— +r
62

where

SS specific storage, H

r a source/sink related term, [L/T]

t time, [T].

(7)

Substituting Eqs. (3) and (4) into Eq. (7) we obtain the following equation

sat 6x xxax 6x ”62 02 2262 62 zxax

Eq. (8) describes groundwater flow in a 2-D anisotropic and heterogeneous media. For
illustration purpose, we assume that the media is homogeneous (e. g., the hydraulic

conductivity is independent of distance).

2 Previous Research

Except for certain special situations under highly restrictive assumptions Eq. (8) must be
solved on a computer using grid-based numerical methods such as a finite difference

(FD) or a finite element (FE) method.

Historically orientation in anisotropy has been accounted for only in finite element
models using unstructured grids. For example, the commonly used finite-element code
FEMWATER [Lin and others, 2001], the newest version of SUTRA [Voss and Provost,
2002], and FEHM [Dash and others, 2002] allow mutually orthogonal anisotropy
directions to be in any orientation. In most finite difference models, the effects of the
cross-conductivity terms are ignored, implicitly assuming that horizontal anisotropy is
unimportant and the aquifer layers are always oriented approximately horizontal. For
example, the popular MODFLOW-2000 [Harbaugh and others,2000] with either the
Layer-Property Flow (LPF) Package [Harbaugh and others, 2000] or the Hydrogeologic-
Unit Flow (HUF) Package [Anderman and Hill, 2000] have been limited to grid-direction
anisotropy. That is, the simulated direction of anisotropy has been limited to being

oriented along the grid axis [USGS open file report 02-409, Anderman E. R. et al., 2002].

Most recently, some FD models such as VST2D [Friedel, 2001] have taken into account
the orientation in anisotropy. However, as we will explain in the following subsection,
these extended FD models, being unable to approximate accurately the cross-conductivity
terms, lead frequently to ill-conditioned matrices and significant errors in the magnitude

and direction of predicted flow, especially for large and complex flow systems.

10

2.1 Traditional FD approach

In this section, we review the traditional finite difference method for solving general

anisotropic flow problems. We fiirther examine the source of the numerical difficulty and

lay the foundation for the development of an improved methodology.

2. 1.1 Transient flow

Most traditional FD methods approximate Eq. (8) by invoking a standard finite difference

approximation. Specifically, the time derivative 9; at each nodal point is approximated

as follow using the following implicit finite difference:

C C
At

III

6(1) ¢n+1_(bn
32—

where,
At
n, n+1
n +1 n
(DC and (DC

(9)

the length of the time step. [T]

the superscripts are the time index, H.

the head in the center cell respectively at time (n+1 ) and (n), [L].

The spatial derivatives are approximated based on centered finite difference using the

following Eq. (10) and (l 1) (Figure 5 represents the finite difference grid with the index

numbering convention).

ll

 

 

 

 

 

 

 

 

 

 

 

 

Z
Al
r
n+1 x
«>21 [D V we /
C 0/
(Dig) ( qDIEFI
‘ n+1 . ’ x
C
/
/ Azl +__\
/ . le — .
(Du-+1 TCDIH'I mil-+1
SW 5 SE

Figure 5. Finite difi'erence grid showing index numbering convention.

For illustration purpose we assume that the media is homogeneous.

 

 

i[K 6(1) ]5K W C E
6x xx xx

 

 

10
6x (Ax)2 ( )
. ,4... WW
— K 5K (11)
6x x2 62 x2 4mm

where,

a)" +1

C the hydraulic head at the center cell at time n +1 ,[L]
(bx/+1, (1):. +1, (DZ.+ 1, (1);,+1 the hydraulic heads at the neighboring cells at time
n +1, [L]
Ax , Az grid spacing in x and z direction, [L]

Other spatial derivative terms of the general Eq. (8) can be approximated in a similar

fashion.

By substituting the finite-difference approximation Eqs. (9) through (1 1) into Eq. (8) and

reorganizing, one obtains the equation:

n+l__ n+1 n+1 n+1 n+1 n+1
(DC —aN<DN +aS<I>S +aE<I>E +aW<DW +aNW<DNW+

n+1 n+1 n+1 n
+aNE¢NE +aSW (DSW +aSE¢SE +b<I>C +s

(12)

where,

n+1 n+1 n+1 n+1
(DN ’(DS ’cDE ’(DW ’
n+1 n+1 n+1 n+1
(DNE ’¢NW’¢SE ’(DSW

the hydraulic heads in the neighboring cells at

time n+l,[L]

s a term related to general sink/source term, [L]

13

The symbols aN ,aS ,aE ,aW ’aNE ’aNW ’aSE ’aSW ,b are weighting coefficrents

givenby:
_ _ 2 .2 2
aN —aS —2AAx (Kx,srn 6+Kz,cos 6) (13)
a = =2AA22(K sin26+K c0826) (14)
E 0W 2' x'
— ~iAxAz K -K ' 20 15
“NW “’SE “ 2 ( x' 2')“ ( )
— -£AxAz(K —K )' 26 16
“NE “’sw ‘ 2 x' z' 5‘“ ( )
2 2
Ax Az
b—ZASST (17)
s =2ArAx 2A2 2 (18)
where,
A_ 0.5m

 

(Ax 2Kx,sin2 9+Ax2Kz,coszc9+Az 2Kx' cosZB+Ax 2K2, sin2 6)2At +SSAx2Az 2

(19)

Equation (12) provides a discrete representation of the governing flow equation. It is a
discrete mathematical statement of the mass balance for a cell over a time interval t" to

tn+l.

l4

2.1.2 Steady state flow

In the case of steady state flow, the numerical representation can be further simplified.

This can be obtained by simply letting S S be zero in Eq. (13) through Eq. (18). The

result is the following FD scheme:

(DC =aN<l>N +aS<DS +aE<1>E +aW<DW +aNW (DNW +aNE¢NE +aSW (DSW +aSE¢SE +3
(20)
with the coefficients given by:
_ _ 2 - 2 2
aN -aS —2AAx (Kx,sm 6+Kz,cos 0) (21)
a = =2AA22(K sin26+K c0526) (22)
E aW z' x'
- —:’i‘-AxAz(K —1< )' 29 23
“NW “'SE ’ 2 x' 2' 3m ( )
- -iAxAz K —K )' 26 24
aNE “’sw ‘ 2 ( x' z' 5‘“ ( )
s =2ArAx 21122 (25)
where
1
A - (26)

 

4(Az 2Kx,coszlSi+AxZKx,sin219-l-Ax2K2,cos26I+Az2Kz,sin2 6)

Equation (12) or (20) indicates that the central head in a cell can be expressed as the
weighted average of the heads in the surrounding cells. The weighting coefficients dictate
the performance of the numerical scheme. These coefficients are a function of orientation
and the ratio of anisotropy, as well as, time step and grid spacing.

15

A basic physical principal requires that these weighting coefficients be positive.
However, when the off-diagonal components of the hydraulic conductivity tensor are not
zero, the weighting coefficients derived fi'om traditional FD methods may not satisfy this
positivity requirement. In fact, as long as the angle of anisotropy (6) is nonzero, the

comer weighting coefficients (a N E ’aSW or a NW ,aS E ) depending on angle of

anisotropy for the central head are always negative. The coefficients become more
negative as the orientation of anisotropy and the ratio of anisotropy increases. The
magnitude of the negative coefficients are largest when the orientation is diagonal (e. g.,
close to 45°). In other words, the head in the center cell would decrease in response to an
increase in the head in the north-west (NW) and south-east (SE) cells. This is not
meaningful and it is the source of numerical problems encountered in the traditional FD

methods.

16

2.2 The MODFLOW approach

2. 2. 1 Transient flow
The popular USGS MODFLOW [MODFLOW-2000, Harbaugh and others,2000]
program circumvents the numerical problem by ignoring the orientation of anisotropy. In

particular, MODFLOW neglects the off-diagonal component of the hydraulic

conductivity tensor (assume 6 = 0’ ), and always works with the following flow equations

(the 2-D flow xz scenario):

5 ELEV ELEV 92),, (27)
3 at 6x xx 6x 62 22 62

With the application of the standard implicit finite difference over time and central finite
difference in space, one obtains the same general numerical scheme in Eq. (12) with the

following weighting coefficients:

_ _ 2
aN —aS —2AAx K2, (28)

_ _ 2
aE —aW —2AAz Kx' (29)
“NE =aNW ="SE =aSW =0 (30)

2 2

Ax Az

b—ZASST (31)
s=2ArAx2A22 (32)

17

where,

A = 1 At 1 (33)

2|:(Ax 2K2, +Az sz,)2At +SSAx2Az 2]

 

2.2.2 Steady state flow

In the case of steady state flow, the weighting coefficients can be further simplified by

assuming SS equal to zero in Eqs. (28) through (32) the result is:

 

_ _ 2

aN —aS —AKZ,Ax (34)
a = =AK A22 (35)
E 0W x'

“NE ”NW ”SE ”SW =0 (36)
s =ArAx 2Az 2 (37)
where

A = 0.5 (38)

szK ,+Az2K ,
Z x

Although the weighting coefficients associated with MODFLOW are all positive and
symmetric, the corner information is not used. Not surprisingly, the distribution of the
weighting coefficients does not properly reflect the impact of the orientation in

anisotropy.

l8

3 improved FD approach

3.1 Transient flow

In this section, we present an improved FD approach for modeling groundwater flow in
general 2-D anisotropic media. Recognizing the numerical difficulty associated with the
cross-conductivity terms, we work with the following governing equation that is free of

the K ,K terms:
xz zx

6(1) 6 6(1) 6 6(1)
Ss aT-lex'a‘vlwile'arl“ (39)

We are able to eliminate the “trouble making” cross-conductivity terms, because the
equation is written in a rotated coordinate system locally aligned with the major direction
of anisotropy. We further invoke implicit finite differencing over time (Eq. 12) as in a
standard FD method and the following central differencing for the spatial terms in the

local x ',z ' system:

 

.n+1 n+1 in+1
_6_[K My ,‘I’w 4% “DE
x x

 

 

4O

ax: an (M32 ( )
n+1 ¢in+l_2¢n+l+®m+l

31 K a") 2K N C S 41

62' z, 62' — z, (A2,)2 ( )

19

where,

 

 

 

 

 

 

 

 

 

 

Kx, K2, the principal hydraulic conductivities, [LT"].
Ax ',Az' the grid spacing in x ' and z ' direction, [L].
.n+1 m+1 m+1 m+1 -
(I) N , (I) S , (I) E , and (DW the non-nodal heads in the rotated
coordinate system, [L] as shown in Figure 6.
I
Z
"+1 n+l
(I) .1
Ni " .
d e
¢n+l
Wt
'ml
(1) W
/ .1
11+
(DSII’

 

 

 

 

 

Figure 6. Nodal and non-nodal head values.

20

Substituting Eqs. (9), (40), and (41) in Eq. (39), we obtain the following:

¢n+1=2A(KZ'AxI2¢In+1+K Axl2¢ln+l+K A212¢1n+1

 

C N z ' S x ' E
i i 42
x ' W At C
where,
A = 0.5At (43)

 

2(Ax 'sz , +Az '2Kx ,)At +SsAx '2Az '2

Note that Eq. (42) cannot be directly used yet since the central head is expressed in terms
of the non-nodal heads. They must be related to the nodal heads before the equation can
be used in an overall numerical analysis. This can be achieved in a straightforward

fashion using the following basic linear interpolations based on the nearest nodal heads

¢§G+l=tan :(¢n+Wl _¢7V+1)+(D;zv+l (44)
<D§1+1=tan622x (<1>"+1— —<D§+1)+<D§+1 (45)
(Pkg-+1=tan6%(<l>n+1—<Dn+l)+¢g,+l (46)
¢£+1=tan6% (o"+1— ¢7V+1)+<I>Z.+1 (47)

Substituting the Eqs. (44) through (47) into Eq. (42) produces the same numerical scheme
as in Eq (12) with the following different weighting coefficients. These coefficients are

for the case of 0° 5 6? < 45°:

21

aN =aS =BKZ,Ax (Ax —Az tan0)

aE =aW =BKx,Az (Az —Ax tan9)

 

 

aNW =aSE =BKZ,AxAz tan6
aNE =aSW =BKx,AxAz tan6

szAzz
b=BSS——2

Atcos (9

BrAx2A22
s = 2
cos 6
where,
1
B =
S Ax2A22

 

2K ,Az2+2K ,Ax2+ S
x Z

cos2 BAt

(48)

(49)

(50)

(51)

(52)

(53)

(54)

Similarly the coefficients could be obtained for other angles (6 ’5). Simply we have to

apply different interpolation equations for each sector of 45 degrees.

22

3.2 Steady state flow

In the case of steady state flow, Eq.(42) reduces to

n+l_ :2 ”1+1 '2 ”1+1
(DC —B(Kz,Ax <I>N +KZ,Ax q’s

 

(55)
+K ,AZ'2¢,n+l+K 'AZI2(DIn+1+rAx12AZr2)
x E x W
where,
1
B = (56)

2(Ax'2K ,+Az'2K )
z x
which is obtained from Eq. (42 ) by setting S s to zero. By substituting the same

interpolation equations in to Eq. (55), we obtain a steady state scheme or Eq. (20) with

the following coefficients:

aN =aS =BKZ,Ax (Ax -Aztan6) (57)
aE =aW =BKx,Az (A2 —Ax tana) (58)
“NW =aSE =BKZ,AxAz tan6 (59)
aNE =aSW =BKx,AxAz tang (60)
where,

B _ 0.5 (61)

 

—K ,A22+K ,sz
x Z

23

Notice that all weighting coefficients obtained from improved FD approach, for any
given anisotropy direction, are positive. Also note how the weighting coefficients vary
proportional to the natural direction of anisotropy. As expected, the center head is always
influenced more by surrounding heads along the direction of orientation than those

across.

As we will demonstrate in the next section, the improved FD method is significantly

more accurate. In addition, the new scheme is more robust than the traditional FD scheme

and always leads to a physically meaningfiil solution.

24

4 Illustrative Examples

In this section, we demonstrate the advantage of the improved F D method with two
specific examples and a systematic comparison with exact solution, when available, and

the solutions obtained using the traditional FD methods.

4.1 A Simple Example

Our first example considers a simple groundwater problem that involves steady state
uniform flow driven by a constant horizontal gradient in a 2-D homogeneous, anisotropic
media with the major direction of anisotropy oriented at 30° with the horizontal axis
(Figure 7). The top and bottom boundary conditions are imposed in such a way that flow

is uniform throughout the solution domain (Figure 7).

25

6(1)(x ,z)

 

 

 

 

 

 

 

 

 

= 0
Z 62 z = l
A 2
Z = I
Z Obs. [:0an
O
6'0 H
u Obs. polnt2 '9'
A II
N O A
o“ N“
V o
6' V
Obs. point 1 0
I
K2. K1"
2 =0 K/ 6=30°
X=O m0”) :0 X=l X
62 z = O x
Figure 7. Definition sketch for the simple example
In this case, the problem has the following exact solution:
“’1 ‘ (D0
<D(x ,z)=(I>O+ l x (62)
x
—K <D(x ,z ) — (I)
v (x,z)= xx[ 0] (63)
x n x
—K <D(x ,z )- (I)
v (x,z)= xz [ O] (64)
n x

26

where,

(”0"”1

the prescribed head on the left and right boundaries, respectively [L].

the horizontal and vertical spatial coordinate measured from the lefi
constant head boundaries and down impermeable boundary, [L]

the length of the domain in the x direction [L]

the length of the domain in the z direction [L]

porosity H

the average linear velocities in x and z direction [LT'I]

We solved the problem numerically using the traditional FD method, improved FD

method, and the MODFLOW method based on the following parameters and inputs:

27

Table I. Input physical and numerical parameters for the simple example

 

 

Parameter Value (unit)
6 30 °
(1)0 10 (m)
(D1 2 ("1)
Ix 800 (m)
[z 600 (m)
K x , 50 (tn/day)
K 2 , 5 (m/day)
n 0.3 H
Ax 5 ("1)
A2 5 (In)

 

The results are summarized graphically in Figures 8 through 11. Obtained solutions for
each method are achieved by utilizing the Interactive Groundwater Model (IGW) [Shu-

Guang Li, et a1, 2002].

28

Z
Z = 600(m)
I
Z
\t
t S
X
\ II
‘x A
x N,
\ O
\ V
‘ '9
t
\
l
h
\
\
I
\
Z = 0 ‘

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

R

Exactvelocity
v =1.29m/day
I .X
vz=0.65mlday
a<D(x,z) _0
02 2:1
7 Z
l/r/ I I’ll/z; .2 2’
12"; I I’ll// I I/ in»
Ix»; / xiv/2w; I // V, '
I
22”- z III/z) I //r,¢' x
l
12'sz I //}///z x {3901
lrr” ‘/ Izlzdzvz // ||
°° h- 1 '0 ’ A
2/w/ I’l'd'd '/ I l/N
I, n I'-
xxx; Ida/(2’ I 1122::
ii 0
xxx-z Izzfipzzz zeal/’99,
xxx; ”ir/l”” //'//e
2/(/ ’I’azzzzzz 2A”
Ill [flail/I’ll 24")
/ w/ /////// zw/z'GlObal
x = 0 5‘99”?) = 0 x =800(m)
& 2:0

 

Figure 8. Exact Solution for the simple example.

29

 

 

 

 

 

 

 

 

 

 

Exact velocity
vx=1.29m/day
v =0.65m/day
Z
62 2:12
Z=600(m) ' *
//’/ I ’II’ I I I,
2222 2 2222 2 2 22 c3\
2' 2222 2 2222 2 2 22 \” '
\ 2222 2 2222 2 2 22,”x
\ p
\\ 222022 2 2222 2 2 {’40}
\ II 2222 ‘2 22224-2122 22 ||
\ 2x co 5 In ’ A
‘\ N2222 2,212 '2; 2M22N
n e a
“ 82222 2 222 2 2 220
O
“ Iel2222 2 22/2 2 2 228
‘ I
\\ 2222 2 /222 2 2 226'
l‘ 2222 ’2’ 2222 2 2 22
‘1 222 36 2222 2 2 22
‘ 2 2222 2 2 22GlObal
z=0‘// 2 2 ‘
x
x = o 6<D(x,z) 0 x =800(m)
& z=0

Figure 9. Numerical solution using improved FD method for the simple example.

30

 

Exact velocity
vx = 1.29m lday

v2 =0.65m lday

 

 

 

 

 

 

Z & :0
2 =1
2 = 600(m)
- ' \
' . . , (90°
Z ,’ '
f I I ’ x
x . .- 2 b’
‘ O (‘1
\ q—d 2 2 2 I]
\\ II I I I A
\ A N
\ N 2 2 2 .
e , » §
\ e ’ ’ v
‘\ I I e
“ I I /
t“ 2'
‘ Global
2 = 0 \ a]
x
x = 0 x =800(m)

  

0

Figure 10. Numerical solution using traditional FD method for the simple example.

31

z
z = 600(m)
z I
\
x O
\ v—fi
\ ll
‘ A
\ N
\‘ o"
‘ V
\ '9
\
\
\
\
\
\
$
\\
z = 0

lllll

9

ollllllll

lllllllllll

1..

.2)

6¢(x

_._q

‘d

n“

‘_

—O-I

 

”d

Z

—o—o
-e-o
—o—o
‘-
—o—o
—o—o

-o—o

‘-

—o-p
0’
5d
fifl

“

fi‘

6<I>(x.z)
a:

 

 

 

 

 

Exact velocity

 

 

 

v =l.29m/day
/ x
vz=0.65m/day
=0
:1
Z
—o—on-o—o —o —o—o \
-.—on-o—o —o -o—o L003
”are—o —o —o-o I,
2’ x
”-4-.“ —. ‘fl
’ C“
n—el-o—o —o —0
2 ll
ndnpfi’T—o —o—o a
4’ ‘5' n 8
#I’1-Qd ‘ —.-.
’ 8
‘9—04-0—0 —o —o-—o e
—o—o(]—o—o -o —o—o
—e—o‘-o—o —o —o—o
_._.._._. _. _._, Global
x
=0 x=800(m)
z=0

 

Figure 11. Numerical solution using the MODFLOWapproach for the simple example.

We have also tabulated the results in Table 2 for the different methods at three

representative points (see Figure 7) to provide a more quantitative comparison.

32

Table 2 - Comparison of the exact predicted head and velocity at three representative

points (Figure 7) with diflerent numerical methods.

Observation point 1

 

 

 

 

 

 

 

 

Method x(m) z(m) <I>(m) vx (m/day) v2 (m/day)
Exact 200 63 8 l .29 0.66
Improved 200 63 8 1 .29 0.65
Traditional 200 63 6.52 1.1 0.23
MODFLOW 200 63 8 l .67 0
Observation point 2
Method x(m) z(m) (I) (m) v v

x {m/day) Z (m/day)
Exact 402 300 5.98 1.29 0.66
Improved 402 300 5.98 1.29 0.65
Traditional 402 300 5.96 1 .21 0.46
MODF LOW 402 300 5.98 1 .67 0
Observation point 3
Method x(m) z(m) (I) (m ) v v

x (m/day) Z (m/day)
Exact 670 526 3 .3 1.29 0.66
Improved 670 526 3 .3 1.29 0.65
Traditional 670 526 4.03 1.52 0.48
MODFLOW 670 526 3.3 1 .97 0

 

33

As shown by Figures 8 through 11 the solution obtained from the improved approach
agrees very well with the exact solution. In fact, they are graphically indistinguishable
from each other. By contrast, the numerical solution obtained from the traditional FD
approach is significantly different from the exact solution for both the predicted heads
and velocities. The MODFLOW approach leads to, not surprisingly, wrong velocity
direction and magnitude since it ignores the orientation of anisotropy. The additional

quantitative comparison in Table 2 leads to the same conclusions.

To further examine why the traditional FD approach, and the MODF LOW approach fail
even for such a simple flow situation, we take a close look at basic numerical
representation of the various approaches for this simple example. Consider a typical
model cell of the discretized domain as shown in Figure 12. As we discussed earlier, the
predicted head in a model cell can be expressed as the weighted average of heads in the
surrounding cells. The values of the weighting coefficients in this case are shown

graphically in the respective cell.

34

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.0262 0.0192 0.262 -0.089 0.148 0.089 0 0.045 0
0.192 0.192 0.352 0.352 0.455 0.455
0.262 0.0192 0.0262 0.089 0.148 -0.089 0 0.045 0

(a) Improved approach (b) Traditional approach (c) MODFLOW approach

Figure 12. Comparison of the weighting coeflicients for (a) the improved FD method,

(b) the traditional FD method, and(c) the MODFLOW method

As we can see, the weighting coefficients for the improved FD approach (Figure 12a) are
all positive and add up to unity. The distribution agrees with what is expected, with the
values being an order of magnitude larger in the northeast (NE) and southwest (S "0 cells

(along the orientation of anisotropy) than in the northwest (NM and southeast (SE).

In the traditional FD approach, although the weighting coefficients along the orientation
still are much larger than those across as they should be, they become negative in the
northwest (NM and southeast (SE) corner cells. This violates the basic physical principle
of groundwater flow since it suggests that an increase in the head in the
northwest/southeast (N W/SE) cell would cause a decrease in the center head. This wrong
hydraulic connection between the central cell and the surrounding cells does not express
the natural hydraulic connection between the cells. Misinterpreting the weighting

coefficients in a model cell causes the wrong solutions in the traditional FD method

(Figure 12b).
35

In the MODFLOW approach (Figure 120), all the coefficients are positive but the corner
information, which expresses the orientation of anisotropy, is missing. This absence of
comer information is what causes the wrong velocity prediction, in the direction and

magnitude.

36

4.2 A More Realistic Example

In this example, we consider a more realistic situation that involves a 2-D vertical flow
obtained from a 3-D model slice, in a three-layer aquifer system characterized by highly
variable elevations and complex sources and sinks (recharge, surface seepage). The
aquifers are assumed to be anisotropic with the bedding oriented along the layers. The
ratio of vertical anisotropy is assumed to be the same in all three aquifers. A conceptual
sketch, along with the boundary conditions, is provided in Figure 13. The input physical

and numerical parameters are provided in Tables 3 and 4.

| I l V. | I I I ' ’ I l l I;
liorltlllvitil:iv:::'lv
I I . ‘ I I | I I 1.: 1 I IIIII I I I
vvvvvv‘vvvvwvhvvvvovvvv

 

a
g a
3 N
3 . E
3 ii
i :3
0 200 400
5:
m

Figure 13. Conceptual sketch with assigned physical and numerical parameters for the

realistic example.

37

Table 3. Physical parameters for three conceptual layers in the realistic example.

 

Parameter

Value (unit)

 

0
Length of cross section
Width of cross section

Depth of cross section

K x , (layer 1)

K x , (layer 2)

K x , (layer 3)

K 2 , /K x , (for each conceptual
layer)

Recharge (throughout whole
domain)

Seepage (drain)

River leakance

Variable based on geology
2000 (m)

40 (m)

Variable around 200 (m)

0.05 (m/day)

2 (m/day)

7 (m/day)

0.1 H

2.086 E (-3) (m/day)

0. 005 (1 /day)

5 (I/day)

 

Table 4. Input numerical parameters for three conceptual layers in the realistic example.

 

 

Input numerical parameters Value (unit)

Steady state condition -

Computational layers for each 2(-)

conceptual layer

Ax 40 (M)

AZ Variable based on geology (m)

 

The problem is solved using both the improved and traditional FD approaches. The

results are shown graphically in Figures 14 and 15.

39

IIIII Illll||
'l:]:::llll‘lI!I,Oll
.. .....

it.“ .
6$$v#&6$%&

‘ 5' V l l I l |
“Macaw“.

No-flow boundary

 

0 200 400
E
m

Figure 14. Results of realistic example applying the improved FD approach .

40

No-flow boundary

Figure 15.

 

" M ~

0 200 400
E
m

Results of realistic example applying traditional FD approach.

41

A qualitative comparison of the results obtained fiom the two different methods reveals
clear differences. In fact, the solution from the traditional FD approach (Figure 15) did
not even converge and the result is wrong or does not conserve mass. We seemed to
frequently have difficulty in obtaining a converged solution with the traditional FD
approach whenever the system becomes relatively complex and the number of degree of
freedom is large. Even when the solution does converge the rate of convergence is often
very slow. These numerical difficulties can all be ultimately attributed to the “negative
weighting coefficients” and the resulting loss of diagonal dominance in the matrix
systems.

The results show that the improved FD approach leads a meaninng solution (Figure 14).
The predicted flow pattern agrees with our intuition. The flow derives from recharge and
is drained toward the lake and the wetlands. The solution also shows the impact of the
variable anisotropy and the bedrock boundary conditions. In particular, as we expected,
the flow follows, to large degree, the orientation of the anisotropy and the shape of the
bottom rock. Numerically, the improved approach is robust and always able to produce a

converged solution given the boundary conditions and hydrologic stresses.

42

5 Conclusion

The following conclusion emerge fi‘om this thesis:

Most real-world aquifers are strongly anisotropic. In many cases, the major
orientation of the anisotropy varies in space and cannot be aligned with the

overall model grid orientation.

Traditional FD methods that account for the variable orientation are often

numerically problematic and lead to inaccurate or physically unrealistic solutions.

The popular MODFLOW approach is able to circumvent the numerical problem
by simply ignoring the orientation of the anisotropy. The resulting solution,
though numerically well behaved, is inaccurate, especially in areas where the
hydraulic gradient is significantly different the natural orientation of the

anisotropy.

The solution obtained from the improved FD method is significantly more
accurate than those obtained from the traditional FD method and the MODFLOW
approach. The improved FD method is robust and can adapt to complex flow

conditions

43

Future research should focus on more thorough testing of the improved FD method in

three-dimensional (3-D) flow situations.

44

6 References

Anderman, E.R., and Hill, MC. (2000). MODFLOW-2000. the US. Geological Survey
modular ggound-water model — Documentation of the Hydrogeologic-Unit Flow
(HUF) Package: US. Geological Survey Open-File Report 00-342.

Anderman, E.R., Kipp Kenneth L., Hill Marcy C., Valstar Johan, and Neupauer
Roseanna M. (2002). MODFLOW-2000. the US. Geologicfl Survev modular
ggound-water model — Documentation of the model-layer variable- irection
horizontal anisotropy (LVDA) capability of the Hydrogeologic-Unit Flow (HUF)
Package: US. Geological Survey Open-File Report 02-409.

Bear, Jacob. (1972). llvrgmics of fluid in porous media New York, American Elsevier.

Dale F. Ritter. (1986). Process geomorphology. Wm.C. Brown Publishers, Dubuque,
Iowa.

Dash, Z.V., Zyvoloski, GA, and Robinson,B.A.. (2002). User’s manual for the FEHM
application a finite-element heat- and mass-transfer code. Los Alamos, New
Mexico, Los Alamos National Laboratory Report LA-13306-M.

Davis, S. N. (1969). Porosity and permeabilig of natural materials. Flow throqu porous
media. ed. R. J. M. De Wiest. Academic Press, New York.

 

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, M.G.. (2000). MODFLOW-
2000, the US. Geologi_c_afl survey modular ground-water model — User guide to
modularization concepts and thground-water flow process. US. Geological
Survey Open-File Report 00-92.

Lin, H.-C. J ., D.R. Richards, C.A. Talbot, G.-T. Yeh, H.-P. Cheng, N.L. Jones. (2001).
FEMWATER: A three-dimensional finite element computer model for simulating
d§_n_sity-de_pgrdent flow; and trapsport in va_r_'i_ablv sfiurated media, Version 3.0.
Technical Report CHL-01-***, US. Army Engineer Waterways Experiment
Station.

Mary P. Anderson, William W. Woessner. (1992). Applied groundwater modeling.
Academic Press, Inc.

Massimo Ferraresi, Alberto Marinelli. (1996). An extended formulation of the integrated
finite difference method for groundwater flow and tra_r_r§port. Journal of
Hydrology. Amsterdam : Elsevier Science B.V. Feb 1996. v.175 (1/4) p. (453-
471)

45

 

Michael J. Friedel. (2001). A model for simulating transient.vafirblv saturated,
coupled water-heat-solute tran_sport in heterogeneous, anisotropic,2-dimensional,
ground-water systems with vapiable fluid density. Water-Resources Investigations
Report 00—4105 Urbana, Illinois.

 

Molson, J .W., and Frind, ED. (2002). WATFLOW/W TC user gpide and
documengtion. Version 3: Waterloo, Ontario, Canada, Department of Earth
Sciences, University of Waterloo.

R.Allan Freeze, John A. Cherry. (1979). Groundwater. Prentice Hall Inc.

Shu-Guang Li. (2002). Interactive ggoundwater model (IGW). Michigan State
University: http://www.egr.msu.edu/~lishug.

Voss, CL, and Provost, A.M.. (2002). SUTRA - A model for saturated-unsaturated

variable-density ggoundwater flow with solute or energy tranaport. US.
Geological Survey Water-Resources Investigations Report 02-4231.

Wang, HF. and Anderson, MP. (1982). Introduction to gpoundwater modeling, finite

difference and finite element methods. San Francisco, W.H Freeman and
Company.

46

   

‘lllmillill'llllllllll‘l