aáëÅêÉíÉ=ÑÉ~íìêÉ=ãçÇÉäáåÖ=çÑ=ÑäçïI=ã~ëë=~åÇ=ÜÉ~í

íê~åëéçêí=éêçÅÉëëÉë=Äó=ìëáåÖ=cbcilt

H.-J. G. Diersch
N
WASY Institute for Water Resources Planning and Systems Research, Berlin, Germany

fluid motion can be defined within such discrete fea-

N aáëÅêÉíÉ=ÑÉ~íìêÉ=ãçÇÉäáåÖ=çÑ=ÑäçïI=ã~ëë=~åÇ=ÜÉ~í=íê~åëéçêí=éêçÅÉëëÉë=Äó=ìëáåÖ=cbcilt

NKN qÜÉ= aáëÅêÉíÉ= cÉ~íìêÉ tures, e.g., Darcy, Hagen-Poiseuille or Manning-Strick-

ler laws. Both the geometric and physical
^ééêç~ÅÜ characteristics of the discrete feature elements provide
a large flexibility in modeling complex situations.
The discrete feature approach provides the crucial
Table 1.1 summarizes the most important characteris-
link between the complex geometries for subsurface
tics and typical applications for the used 1D and 2D (as
and surface continua in modeling flow, contaminant
well as 3D porous media) features.
mass and heat transport processes. In this holistic
approach a three-dimensional geometry of the subsur-
Apparently, the range of applications and the
face domain (aquifer system, rock masses) in describ-
dimension of the features require an unified approach,
ing a porous matrix structure can be combined by
where linear and nonlinear laws of fluid motion, porous
interconnected one-dimensional and/or two-dimen-
media and free fluid flows, phreatic and non-phreatic
sional features as shown in Fig. 1.1. In the finite-ele-
conditions as well as spatial (3D), plane (1D, 2D) and
ment context the three-dimensional mesh for the
axisymmetric (1D) geometries are embodied.
porous matrix can be enriched by both ’bar’ (channels,
mine stopes) and areal (overland, fault) elements.

NKP mêÉäáãáå~êáÉë
NKO qÜÉ= Na= ~åÇ= Oa= aáëÅêÉíÉ NKPKN cìåÇ~ãÉåí~ä= Ä~ä~åÅÉ= ëí~íÉJ
cÉ~íìêÉ=bäÉãÉåíë=rëÉÇ ãÉåí
FEFLOW2 provides 1D and 2D discrete feature ele-
ments which can be mixed with the porous matrix ele- The conservation of mass, momentum and energy is
ments in two and three dimensions. Different laws of described by the balance statement2 (symbols are listed

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=N
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
in the Appendix A ’Nomenclature’)

--------------
  -
+    v  +   j = f (1-1)
t

conserving the (extensive) quantity    . Individual

balance laws for    , j and f are summarized in
Tab. 1.2.

1D channel element
river

runoff surface

2D fracture element

unsaturated zone well

fault

3D porous matrix element

aquifer

saturated zone

Figure 1.1 Schematization of a subsurface modeling system by combining discrete feature elements with volume dis-
cretizations of the total study domain: 1D elements are used to approximate rivers, channels, wells and specific faults,
2D feature elements are appropriate for modeling runoff processes, fractured surfaces and faulty zones, and 3D ele-
ments represents the basic tessellation of the subsurface domain consisting of an aquifer-aquitard system and involving
unsaturated and saturated zones.

O=ö=cbcilt
 NKP=mêÉäáãáå~êáÉë

NKPKO cçêãë=çÑ=Ä~ä~åÅÉ=Éèì~íáçåë
Table 1.1 Used discrete feature elements
According to the applications for the discrete fea-
Fluid motion ture elements indicated above we are interested in four
Type Dimension Application forms of the governing balance equation (1-1):
law

1D, plane channels • form A: free fluid balance law

Darcy (phreatic, non- mine stopes • form B: vertically integrated free fluid balance
Hagen-Poi- phreatic) law
seuille • form C: porous medium balance law
Manning- 1D, axisym- pumping wells
metric abandoned • form D: vertically integrated porous medium bal-
Strickler
(phreatic, non- wells, bore- ance law
phreatic) holes
The form A is already represented by Eq. (1-1).
Darcy 2D, plane fractures
A vertical integration of (1-1) over a depth B can be
Hagen-Poi- (non-phreatic) faults
rigorously performed as described in2,3,8 leading to the
seuille
Manning- 2D, plane runoff form B:
Strickler (phreatic) overland flow
------------------
 B - top bottom
+    Bv  +    Bj  = Bf + j  – j  (1-2)
t
3D porous media
Darcy (phreatic, non- aquifer systems with the new exchange terms of the quantity  at the
phreatic)
top and bottom boundaries

Table 1.2 Balance laws

top 1 top

Quantity  j f
j  = ------
S  n   j +   w – v   dS

S
top

 (1-3)
mass 1
  j +   w – v   dS 
bottom bottom
j = ------  n
fluid mass  0 Q  S
S
bottom 
contami-
nant mass C jc rc

momentum v  g Notice, the balance quantities of Eq. (1-2) are now

energy   E + --- v 
1 2
  v + jT   g  v + QT  averaged over the depth B .
2 The transformation of the balance equation (1-1) to a
porous medium is performed by a spatial averaging

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=P
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
procedures referred to the representative elementary
 1 non-integrated form
volume (REV) composed by fluid and solid phases. It B =  for (1-8)
 arbitrary vertically integrated form
results finally to the form C of the basic balance
statement3
interface
the interface exchange term j  as
-----------------
  - interface
+    v  +    j  = f + j  (1-4)
t
interface  0 free fluid flow
j =  for (1-9)
 0 porous media flow
where an exchange term at the fluid-solid interface nat-
urally results
top bottom
and the top and bottom exchange terms j   j  as
interface 1 interface
j = ------
S  n   j +   w – v   dS (1-5)
 0 non-integrated form
interface top bottom
S  j  j  =  for (1-10)
 0 vertically integrated form
Notice, the balance quantities of the porous medium
conservation equation (1-4) are averaged over the REV
volume. NKPKP j~íÜÉã~íáÅ~ä=ÅçåîÉåíáçåë
Finally, the porous medium equation (1-4) can also Both Cartesian and cylindrical coordinate systems
be vertically integrated over the depth B , which yields will be employed. They are defined as
form D of the basic balance statement as

---------------------
 B -  x y z   3D
+    Bv  +    Bj  = (1-6)   
t  x y   2D
interface top bottom x =   for  (1-11)
Bf + j  + j – j  x   1D
 r  z   axisymmetry
  
It is obvious, the balance statement (1-6) of form D is
the most general form which encompasses all other The velocity vector v is accordingly
forms when we specify the porosity  as

 1 free fluid flow

 =  for (1-7)
 1 porous media flow

the depth B as

Q=ö=cbcilt
 NKP=mêÉäáãáå~êáÉë

2 2 2
  ---------
  - ---------  
 u  2- + --------- 2
+ 2- 3D (x y z) Cartesian
 Cartesian coordinates  x y z
 v  2 2
 w     2D (x y) Cartesian
  ---------
2
- + ----------
2
v =  for (1-12) 2  x y
 vr   =  (1-14)
  2 
 v cylindrical coordinates  ---------- 1D (x) Cartesian
  x 2
 vz 
  1    1  2   2 
 --- -----  r ------- + ----2 ---------2- + ---------
- cylindrical (r  z)
 r r r r  z
2

The scalar product   v is given by

 u v w NKPKQ dê~îáíó=~åÇ=î~êá~ÄäÉë
 ------ + ----- + ------- 3D (x y z) Cartesian
 x y z
 u v In the following we assume an exclusive action of
 ------ + ----- 2D (x y) Cartesian
 x y gravity in the form
  v =  (1-13)
 u
------
 x 1D (x) Cartesian
 gx ex
 1   rv r  1 v  v z g = – ge g = gy e = ey (1-15)
 --r- --------------
r
- + --- --------- + --------
r  z cylindrical (r  z)
 gz ez

2
and the derivative operation  for the different coor-
As a further useful variable the hydraulic head h
dinate systems is given, for instance, for the variable 
(piezometric head) related to the reference fluid density
as
 o is defined

p
h =  + z = --------- + z (1-16)
o g
and

p = o g  h – z  (1-17)
Thus,

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=R
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
 – o
p – g = p + ge =  o g   + z + --------------- e (1-18)
o
Table 1.3 Hydraulic radii for different applications
 – 
=  o g  h + --------------- e =  o g  h + e 
o
o Type r hydr

submerged rectangular cross-

A section
NKPKR eóÇê~ìäáÅ=ê~Çáìë
Bb -
B --------------------
The hydraulic radius is defined as the flow cross- 2 b + B
sectional area divided by the wetted perimeter b

submerged slit plane

flow area
r hydr = ---------------------------------------- (1-19) B
wetted perimeter
b
bB
B ------- = b---
2B 2
Table 1.3 lists the hydraulic radii for interesting cases.
open rectangular cross-section
C
NKPKS cêÉÉ= EéÜêÉ~íáÅF= ëìêÑ~ÅÉ= ÅçåÇáJ B Bb
----------------
íáçå b
b + 2B

A phreatic surface represents a macroscopic moving open wide channel (b > 20B)
material interface between two fluids, e.g. air and D plane
water. A material surface F = F  x t  = 0 is governed B
B -
by the kinematic equation ---------------------- B
b 1 + 2B  b
F
------ + w  F = 0 (1-20) submerged circular cross-section
t E
R
2
R - R
--------- = ---
2R 2
The outward unit vector normal to F is defined as

F and accordingly
n = ------------ (1-21)
F
F  t
w  n = – --------------- (1-22)
F

S=ö=cbcilt
 NKP=mêÉäáãáå~êáÉë

---------------------
 B -     h
= B ------------------ +  ------ (1-26)
where F denotes the magnitude of the vector F . t t t
For the vertical integration along the thickness B we
can express the geometries of the top and bottom sur-
faces in the forms (Fig. 1.2) NKPKT sáëÅçìë=ëíêÉëëÉë=çå=ëìêÑ~ÅÉë
top top top
F  F  x t  = z – b  x y t  = 0  The viscous stresses on a surface  (note  can
bottom bottom bottom
 (1-23)
F F  x t  = z – b  x y  t  = 0  indicate a top and bottom surface as well as a fluid-
solid interface) result from exchange relationships (1-
3) and (1-5) if replacing the general flux vector j by
and
the viscous stress tensor of fluid  (cf. Table 1.2), viz.,
top bottom
B = B  x t  = b  x y t  – b  x y t  (1-24)  1 
 = ------
S n    + v  w – v   dS (1-27)

S
top
F = 0 n

Here  stands for the stress on the surface  with nor-

B (x,y,t) mal n . It represents a surface force per unit area
bottom h depending on the orientation of the surface10. For
F = 0 instance, let us consider the stress components on a pla-
top nar top surface as illustrated in Fig. 1.3. Assuming
z b
bottom
 x y t  b  x y  t 
y additionally a rigid and impermeable surface
x datum level
( w = v  0 ) with a constant stress property on the unit
area S the surface stresses are explicitly given by
Figure 1.2 Surface conditions. top top
 n  (1-28)

top
For a free surface the top elevation z = b  x y t  is With n
top
=  0 1 0  the stress components become
identical to the hydraulic head h = h  x y t  . Accord-
ingly, the thickness is given by
= 0 xx + 1 yx + 0 zx =  yx 
top
x

bottom top 
B = h–b (1-25) y = 0 xy + 1 yy + 0 zy =  yy  (1-29)

= 0 xz + 1 yz + 0 zz =  yz 
top
z
If the bottom geometry is stationary the storage term in
(1-6) becomes

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=T
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
  B  h
=   BS o + S s  ------ (1-32)
top t t
y =  yy

0
n
top
= 1
where the compressibility S o and storativity S s are
0
introduced as
top
x =  yx
y S o =  +  1 –   
 (1-33)
top
z =  yz Ss =  
x
z

Figure 1.3 Surface forces related to components of Notice, for a free fluid we have to set  = 1 and
the viscous stress tensor  . So =  .
Neglecting the density effects in the divergence term of
(1-30) by applying the Boussinesq approximation2 the
fluid mass balance equation (1-30) yields
NKQ _~ëáÅ=_~ä~åÅÉ=bèì~íáçåë
h 
NKQKN cäìáÇ=ã~ëë=ÅçåëÉêî~íáçå S ------ +    Bv  = BQ  
t  (1-34)
S =  BS o + S s  
The fluid mass conservation is described by speci- 
fying Eq. (1-6) with Table 1.2 as

  B  +    Bv 
t
= BQ  (1-30) NKQKO cäìáÇ=ãçãÉåíìã=ÅçåëÉêî~íáçå

The fluid momentum conservation is specified from

which can be employed for all flow problems under
Eq. (1-6) with Table 1.2 as
discussion when setting  and B appropriately. Notice,
the sink/source term Q  includes both interfacial and 
surfacial flux conditions (cf., Eq. (1-5)). ----  Bv  +    Bvv  = –   Bp  (1-35)
t
The storage term in (1-30) +    B'  + Bg + B  
interface
+
top
–
bottom

-
  B  = B -----  B
+ B ----- +  ------ (1-31)
t t t t where the stress tensor is splitted into the equilibrium
(pressure) and non-equilibrium (deviatory) parts as
can be expanded with regard to the hydraulic head h
and one gets with (1-26)  = – pI + ' (1-36)

U=ö=cbcilt
 NKQ=_~ëáÅ=_~ä~åÅÉ=bèì~íáçåë

v
interface ----- 0  v   v  0 (1-40)
Notice, in Eq. (1-35) the exchange term  t
vanishes for free fluid motion and the terms
top bottom
  are dropped if the equation is not vertically As the result, one yields a general momentum equation
integrated. for porous media (we consider the non-integrated form
In the following we assume the Newton’s viscosity law top
with B  1  = 
bottom
= 0 ) as
(including the Stokes’s assumption10) which is written
in the form   p – g  = 
interface
+  v
2
(1-41)

1
' = 2 d – ---    v I (1-37) Furthermore, the drag forces due to fluid viscosity can
3 2
usually be dropped  v  0 with respect to the drag
with the strain-rate tensor term of momentum exchange 
interface
at interfaces of
phases. The interfacial drag term of momentum
1 T interface
d = ---  v +  v   (1-38) exchange  can be derived as a linear friction
2
relationship of the form3
For an incompressible fluid with a divergenceless (so- interface –1
 = – k   v  (1-42)
called solenoidal) velocity   v = 0 the momentum
equation (1-35) leads to the well-known Navier-Stokes
equation where the permeability k represents an inverse friction
tensor due to the viscous drag at the interfaces of fluid-
v 2 solid phases.
B ----- +  Bv   v = – B  p – g  + B v (1-39)
t Finally, the momentum equation (1-41) reduces to the
interface top bottom well-known Darcy equation of the form
+ B   + – 

k
v = – ------  p – g  (1-43a)
from where specific forms can be derived as follows. 
or with (1-18)

NKQKOKN a~êÅó=Ñäçï=áå=éçêçìë=ãÉÇá~
v = – Kf   h + e  

Commonly, in a porous medium the velocity v is k o g 
K = ------------ 
sufficiently small, that means the Reynolds number o  (1-43b)
based on a typical pore diameter is of order unity or 
o 
smaller. As the result, the inertial effects in the momen- f = -----
- 
 
tum equation (1-39) can be neglected

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=V
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
valid for flow in a porous medium. dp
2
------ – g x =  --------u
d
(1-45)
dx dy
2

NKQKOKO mä~åÉ=~åÇ=~ñáëóããÉíêáÅ=é~ê~ääÉä=
EmçáëÉìáääÉF=Ñäçï Integrating (1-45) with the boundary conditions
u  0  = u  b  = 0 it yields
A flow is called parallel when inertial terms of the
Navier-Stokes equation (1-39) vanishes. That means, a 1 dp
u = – ------  ------ – g x y  b – y  (1-46)
fluid particle is subjected to zero acceleration, accord- 2 dx
ingly, it moves in pure translation with constant veloc-
ity. It follows that pathlines must be straight lines and and we obtain the average velocity in the aperture b as
that the velocity of each particle may depend only on
coordinates perpendicular to the direction of flow. Such b
2
1 b dp
flow fields occur between two parallel plates or in a u = ---  u dy = – ---------  ------ – g x (1-47)
b 12  dx 
circular tube as depicted in Fig. 1.4. y=0

a) b)
y r and the discharge Q

3
b dp
b
R Q = ub = – ---------  ------ – g x (1-48)
12 dx
u vz z

x which is called the cubic law of the Hagen-Poiseuille

Figure 1.4 a) 2D plane and b) axisymmetric Poiseuille flow. The relationships (1-47) can be expressed by the
flow. hydraulic radius r hydr if replacing the dimension b  2
for the slit flow according to Table 1.3 (type B)
For 2D parallel laminar flow (Fig. 1.4a) we have 2
r hydr dp
u = – -----------  ------ – g x (1-49)
3 dx
u
v = v u = u y  v = w = 0 (1-44)
w Similarly, for the axisymmetric flow in a circular tube
(Fig. 1.4b) with

and the momentum equation (1-39) in the x-direction

becomes (we consider the free fluid case with no verti-
cal integration)

NM=ö=cbcilt
 NKQ=_~ëáÅ=_~ä~åÅÉ=bèì~íáçåë

vr
As seen the Hagen-Poiseuille’s laws of laminar fluid
v = v vz = vz  r  vr = v = 0 (1-50) motion for 1D and 2D plane flow (1-47) and for axi-
vz symmetric flow (1-53) represent linear relationships
with respect to the pressure gradient and gravity
 p – g  . In a generalized form one yields finally
one solves in the z-direction the momentum equation with (1-18)

dp   v
------ – g z = --- -----  r --------z (1-51) v = – Kf   h + e 
dz r r  r  (1-56)
2
r hydr  o g  r hydr = b  2 a = 3 for 1D/2D plane
With dv z  dr = 0 at r = 0 and v z  R  = 0 the integra- K = -------------------- I with 
a o  r hydr = R  2 a = 2 for axisymmetry
tion of (1-51) gives

1 dp
v z = – ------  ------ – g z  R – r 
2 2
(1-52)
4 dz NKQKOKP i~ïë=çÑ=ÑäìáÇ=ãçíáçå=Ñçê=çîÉêä~åÇ=
~åÇ=ÅÜ~ååÉä=Ñäçï
Then, the average velocity for the Hagen-Poiseuille
flow in a circular tube becomes Basically, the fluid motion for overland and channel
flow is described by the vertically integrated Navier-
2 R
2 Stokes equation (1-39) according to
1 R dp
v z = ---------2   v z r dr d = – ------  ------ – g z (1-53)
R 8  dz  v
=0 r=0 B ----- +  Bv   v = – B  p – g  (1-57)
t
2 top bottom
+ B v + B   – 
and we get for the discharge through the tube

4
R dp which is a formulation of the well-known De Saint-
Q = R v z = – ---------  ------ – g z
2
(1-54)
8 dz Venant equations1. Over a wide range of practical over-
land and channel flow (Fig. 1.5) at low-to-moderate
The relationships (1-53) can be expressed by the velocity/flow regimes the inertial terms in the govern-
hydraulic radius r hydr if replacing the dimension R  2 ing momentum balance equation (1-35) can be ignored
for the tube flow according to Table 1.3 (type E) compared with the gravitational terms, friction and
pressure effects. Furthermore, the interior viscous
2 effects can be neglected over the shear stress effects at
r hydr dp
v z = – -----------  ------ – g z (1-55) the surfaces1,9. Assuming this,
2 dz

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=NN
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
v 2 following momentum equation can be derived
-----  0  v   v  0  v  0 (1-58)
t
 p – g  +  o gS f = 0 

the momentum equation (1-57) reduces to v v-  (1-62)
S f = ---------------
2  
 r hydr 
top bottom
 p – g  –  + = 0 (1-59)

where different specific laws for the friction slopes S f

z can be specified as summarized in Table 1.4 for isotro-
g pic roughness coefficients.


h  top Table 1.4 Various friction laws

w in
v d str
es s
Law   Sf
 bottom
datum level slop
e frict
x,y ion
 v v-
Newton-Taylor --g- 1 -------------
 gr hydr
Figure 1.5 Open channel flow.

v v-
Chezy C 1 -----------------
top 2
The shear effect  at the top (free) surface can be C r hydr
caused by wind stress. For the present application we
neglect influences caused by wind stress: v v-
Manning-Strickler M 43 ------------------
2 43
M r hydr
top
 0 (1-60)

On the other hand, the shear effect at the bottom can be Instead of using the pressure p as primary variable
expressed by a friction slope relationship of the form the hydraulic head h or the local water depth  (cf.
Eqs. (1-16) and (1-18)) are alternative formulations of
bottom o g v v (1-62), viz.,
 = -------------------
2 
(1-61)
 r hydr
 o g  h + S f + e  = 0 (1-63a)
and
representing friction laws, where v = v  v ,  is a
friction factor and   1 is a constant. As the result, the

NO=ö=cbcilt
 NKQ=_~ëáÅ=_~ä~åÅÉ=bèì~íáçåë

 o g   + S f – S o + e  = 0  It can be easily shown that the velocity v in the rela-

 (1-63b) tionship (1-67) tends to zero if the gradient h van-
S o = – z  ishes (provided that   0 )

2
where fluid density effects are included in the e - r hydr
term. lim v = – lim ------------------- I h = 0 (1-68)
h  0 h  0 4 2
h
Equation (1-63a) can be used to derive a diffusion-
type flow equation7,9. Since, exemplified for 2D
NKQKP `çåí~ãáå~åí= ã~ëë= ÅçåëÉêî~J
v
2
= u +v =
2 2 2 
 r hydr
2
S fx +
2
S fy (1-64) íáçå
and with (1-62)
The balance equation for a contaminant mass results
from Eq. (1-6) and Table 1.2 in form of
2 2 2 2
u + v- u + v-
S fx = -------------------- u S fy = -------------------- v (1-65)
2 
 r hydr
2 
 r hydr ------------------
 BC 
+    BCv  +    Bj c  = Br c (1-69)
t

we find with (1-63a): S fx = –  h  x + e x  ,

which can be employed for all the interesting mass
S fy = –  h  y + e y 
transport problems when specifying  and B appropri-
2 ately. Notice, the reaction term r c includes both inter-
r hydr h facial and surfacial mass transfer conditions (cf., Eq.
u = – ----------------------------------------  ------ + e x (1-66)
2 2  x  (1-5)).
 h
4 ------ + ------
 h
 x  y
2 The reaction term can be splitted into a first-order
r hydr h-
 -----
v = – ---------------------------------------- + e y reaction rate and a zero-order production term2, respec-
2 2  y 
 h
4 ------ + ------
 h tively,
 x  y
r c = – C + Q c (1-70)
or more general
The mass flux j c is expressed by the Fickian law in
v = – K  h + e  

form of
2
r hydr  (1-67)
K = ------------------- I 
4 2
 j c = – D  C 
h  (1-71)
D = Dd I + Dm 

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=NP
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
 =  +  1 –    C  
The hydrodynamic dispersion tensor D consists of the 
(1-76)
molecular diffusion part D d I and the mechanical dis-  d =  +  1 –   d    C   C  
persion part D m . In a porous medium D m is commonly dC 
described by the Scheidegger-Bear dispersion relation-
ship as in which the sorption function   C  can be specified
for Henry, Freundlich or Langmuir isotherms2.
vv
D m =   T v I +   L –  T  ------------ (1-72)
v

In a free fluid flow there is a large variety for Dm in

NKQKQ båÉêÖó=ÅçåëÉêî~íáçå
dependence on laminar and turbulent flow conditions.
The energy balance equation is derived basically from
For instance, in a fluid-filled tube under laminar flow
Eq. (1-6) and Table 1.2 under the assumption of a ther-
conditions D m can by estimated by Taylor’s analysis11
mal equilibrium between fluid (f) and solid (s) phases.
We obtain finally2
 R2 v  v  v
D m =  ------------- ------------ (1-73)
 48D d  v 
---- f f s s f f
 B   E +  1 –    E   +    B E v  (1-77)
t
f f s s
+    Bj T  = B   Q T +  1 –   Q T 
Using the Fickian law (1-71) and incorporating the
continuity equation (1-30) Eq. (1-69) yields2
which can be applied to all the interesting heat trans-
C
B ------- + Bv   C –    BD  C  (1-74) port problems when specifying  and B appropriately.
t f s
Notice, the thermal sink/source terms Q T Q T include
+ B  Q  +  C = BQ c both interfacial and surfacial heat transfer conditions
(cf., Eq. (1-5)).
Considering additionally sorption effects in the porous
medium the following contaminant mass transport Using the state relation for the internal energy2
equation can be derived on the basis of Eq. (1-74)2  
dE = c dT for  = s f (1-78)
C
B d ------- + Bv   C –    BD  C  (1-75)
t and the Fourierian heat flux as
+ B  Q  +  C = BQ c

with the retardation relationships

NQ=ö=cbcilt
 NKR=dÉåÉê~äáòÉÇ=jçÇÉä=bèì~íáçåë

j T = –   T 

 (1-79)
 = 
cond
+
disp f s f f
=   +  1 –   I +  c D m  NKRKO `çåí~ãáå~åí=ã~ëë

The governing contaminant mass transport equation
one yields the following balance equation for the ther- (1-75) can now be specified for the different flow con-
mal energy2 ditions and discrete feature elements. Table 1.6 summa-
rizes the different terms and expressions for both
f f s s T f f porous media and free fluid conditions.
 B   c +  1 –    c   ------ +  c Bv  T
t
f f
–    B  T  + B c Q   T – T o  (1-80)
NKRKP eÉ~í
f f s s
= B   Q T +  1 –   Q T 
The specified terms for the governing heat transport
to be solved for the system temperature T . equation (1-80) are summarized in Table 1.7 for both
porous media and free fluid conditions.

NKR dÉåÉê~äáòÉÇ= jçÇÉä= bèì~J

íáçåë
NKRKN cäçï

The fundamental flow equation represents a combi-

nation of the fluid mass conservation equation (1-34)
and the fluid momentum conservations for porous
media (1-43b), Poiseuille flow (1-56) and overland/
channel flow (1-67). As the result, Table 1.5 summa-
rizes the governing equation for the used discrete fea-
ture elements for 1D, 2D and 3D in dependence on the
problem cases under consideration. For the Poiseuille
flow and overland/channel flow standard geometric
forms of the fractures are implemented in FEFLOW.
Different geometries can be input by means of correc-
tions in the corresponding hydraulic parameters as
thoroughly described in Appendix D.

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=NR
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
Table 1.5 Flow model equations
h
L  h  = S ------ –    Kf  B   h + e   – Q = 0
t
S Kf  B Q

Case Darcy Poiseuille Overland Darcy Poiseuille Overland Darcy Poiseuille Overland

k o g 2 2
bB -----------
- f hydr  o gI f 
bB r----------------------- r hydr I
1DPP b  BS o + S s  b  B + 1  b  B + 1  o bB ------------------
- f bBQ  bBQ bBQ 
3 o 4 2
(1D plane phre- h










K






atic) K K

2 2
hydr  o gI f 
k o g bB r----------------------- r hydr I
1DPN bB -----------
- f bB ------------------
- f
bBS o bB bB o 3 o bBQ  bBQ bBQ 
4 2
(1D plane non- h










K






phreatic) K K

2 2
Ss 2 k o g r hydr  o gI f 2 r hydr I
1 1 R ----------- 2
1DAP R  S o + ---- R   + --- R   + --- - f R ----------------------- 
2 2 2 2 2 2
o  R ------------------
- f R Q  R Q  R Q 
B B B 2 o 4 2
(1D axisym- h










K






metric phreatic) K K

2 2
2 k o g f r hydr  o gI f 2 r hydr I
1DAN 2
R S o R 
2
R 
2
R -----------
-  R
2
-----------------------  R ------------------
- f 2
R Q 
2
R Q 
2
R Q 
o 2 o 4 2
(1D axisym-  h














metric non- K K K
phreatic)

2 2
k o g f
B ----------- hydr  o gI f 
B r----------------------- r hydr I
2DPP BS o + S s B + 1 B + 1 -  B ------------------
- f BQ  BQ  BQ 
o 3 o 4 2
(2D plane phre- h















atic) K K K

2 2
k o g f
B ----------- hydr  o gI f 
B r----------------------- r hydr I
2DPN BS o B B -  B ------------------
- f BQ  BQ  BQ 
o 3 o 4 2
(2D plane non- h















phreatic) K K K

3DP 2 2
(3D phreatic k o g f r hydr  o gI f r hydr I
So   ------------  -----------------------  ------------------
- f Q  Q Q
o 3 o
and non-phre- 4
h
2















atic) K K K

NS=ö=cbcilt
 NKR=dÉåÉê~äáòÉÇ=jçÇÉä=bèì~íáçåë

Table 1.6 Contaminant mass transport model equations

C
L  C  = S ------- + q  C –    BD  C  + C – Q = 0
t
S q BD  Q

free free free

Case porous porous porous free fluid porous free fluid porous
fluid fluid fluid

1DPP bB d bB bBv bBv bB  D d I + D m  bB  D d I + D m  bB  Q  +   bB  Q  +   bBQ c bBQ c

(1D plane phreatic
and non-phreatic)

1DAP 2
R  d R
2 2
R v
2
R v
2
R   D d I + D m 
2
R  D d I + D m 
2
R  Q  +  
2
R  Q  +  
2
R Q c
2
R Q c
(1D axisymmetric
phreatic and
non-phreatic)

2DPP B d B Bv Bv B  D d I + D m  B  Dd I + Dm  B  Q  +   B  Q +   BQ c BQ c

(2D plane phreatic
and non-phreatic)

3DP d 1 v v   Dd I + Dm  Dd I + Dm Q  +  Q +  Q c Qc
(3D phreatic and
non-phreatic)

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=NT
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
Table 1.7 Heat transport model equations
T
L  T  = S ------ + q  T –    B  T  +   T – T o  – Q = 0
t
S q B  Q

free
Case porous free fluid porous free fluid porous free fluid porous free fluid porous
fluid

f f f f f
1DPP bB   c + f f f f f f bB    +
f bB   I + f f f f bB   Q T + f f
bB c bB c v bB c v f f
bB c Q  bB c Q  bB Q T
(1D plane phre-  1 –   c 
s s s
 1 –   I +  c Dm   1 –   Q T 
s s

atic and non- f f

 c D m 
phreatic)

2 f f 2 f 2 f f
1DAP R   c + 2 f R   I + R   Q T +
2 f f
R  c
2
R  c v
f f 2 f f
R  c v R    + 2 f f
R  c Q 
2 f f
R  c Q 
2 f
R  Q T
f
f f
(1D axisym-  1 –   c 
s s
s  c Dm   1 –  
s s
QT 
 1 –   I +
metric phreatic f f
and  c D m 
non-phreatic)

f f f f f
2DPP B   c + f f f f f f B    +
f B I + f f f f B   Q T + f f
B c B c v B c v f f
B c Q  B c Q  B Q T
(2D plane phre-  1 –   c 
s s s
 1 –   I +  c Dm   1 –   Q T 
s s

atic and non- f f

 c D m 
phreatic)

f f f f f
3DP  c + f f f f f f   + f f f f f f f  Q T + f f
c  c v cv  I +  c Dm  c Q   c Q  QT
s
(3D phreatic  1 –   c
s s
 1 –   I +  1 –   Q T
s s

and non-phre- f f
 c D m
atic)

NU=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

NKS cáåáíÉ= bäÉãÉåí= cçêãìä~J given. In (1-82) n corresponds to the normal unit vec-
íáçåë tor (positive outward),  1 and  2 are prescribed
boundary values of  on  1 and  2 , respectively.
NKSKN j~ëíÉê= Éèì~íáçåI= ÄçìåÇ~êó
The finite element formulation is based on the weak
ÅçåÇáíáçåë=~åÇ=ïÉ~â=ëí~íÉãÉåí form of the basic equation (1-81). Introducing a spatial
weighting function w we get
The governing balance equations as listed in Tables
1.5, 1.6 and 1.7 can be generalized by the following 
master equation  w  S ------
t
- + q   d
 (1-83)


L    = S ------- + q   –    D    – Q  = 0
t
=  w     D    + Q  d
(1-81) 
Q  = –  + Q
Applying partial integration and the divergence theo-
which has to be solved for flow (  = h ), contaminant rem (Green’s theorem) to the weak statement (1-83)
mass (  = C ) and heat (  = T ) for 1D, 2D and 3D and inserting the Robin-type BC (1-82) the following
discrete feature elements. weak form for the finite element method finally results

D 
Let   R and  0 T t  be the spatial and temporal
 w  S ------- + q   + w   D    d (1-84)
domain, respectively, where D is the number of space t

dimension (1, 2 or 3) and T t is the final time, and let
 =  1   2 denote the boundary of  , where  1
+  wa d =  wQ d +  w  a2 – b  d
2  2
and  2 are two disjoint portions of the total boundary,
 , the following boundary conditions (B.C.’s) have
to be appended to (1-81):
NKSKO pé~íá~ä=ÇáëÅêÉíáò~íáçå
 = 1 on 1 
 (1-82) In the finite element context a spatial semi-discreti-
– n   D    + a   2 –   = b on 2  h
zation  of the continuum domain  is achieved by
e
the union of a set of non-overlapping subdomains  ,
where on  1 we have Dirichlet BC and on  2 it repre- the finite elements, as
sents a more general form of a Robin type BC in which
more specific Neumann and Cauchy type BC’s are h
    
e
(1-85)
involved. If a = 0 a Neumann BC of 2nd kind results,
e
while for b = 0 a common Cauchy BC of 3rd kind is

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=NV
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
e
On any finite-element domain  , the unknown vari- polynomials that are piecewise-continuously differen-
able  (and dependent coefficients) are replaced by a tiable and square integrable (but whose second and
continuous approximation that assumes the separability higher derivatives need not to exist).
of space and time, thus
Using the Galerkin-based finite element method
h where the test function w becomes identical to the trial
  x t     x t  = N i  x  i  t  (1-86)
space N , Eq. (1-84) leads to the following global
matrix system of M equations
where i = 1  M designates nodal indices, M is the
total number of nodes, N i is the nodal basis function, ·
O+K–F = 0 (1-87)
called the trial space, and x are the spatial coordinates
(1-11). Note that the summation convention is used for
repeated indices. For the present analysis the basis with its components written in indicial notation
functions Ni are based on C 0 (continuous) piece-wise
e 
O ij =  Oij =   SNi Nj d 

e e e


e 
K ij =  K ij =    N i q  N j + N i   D  N j  + N i N j  d +  aN i N j d 
 (1-88)
e e e e
2 


e
F i =  F i =   N i Q d +  N i  a 2 – b  d 

e e e 2
e 

where the subscripts i j = 1  M denote nodal indi- cally only solved by numerical schemes. For stability
ces. The superposed dot in (1-87) means differentiation reasons implicit (A-stable) two-step techniques are pre-
with respect to time t , viz., ferred.

· d  Considering   t  within the finite interval

  t  =  -----   t   (1-89)  t n t n + t n  , where the subscript n denotes the time
 dt 
plane and t n is a variable time step length, the func-
tion   t  is defined as
NKSKP qÉãéçê~ä=ÇáëÅêÉíáò~íáçå
n
 =   tn  (1-90)
The spatially discretized equation (1-87) as a com-
mon first-order differential equation in time can practi-

OM=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

at the previous (old) time plane and as O-

 ------- O-
 n + 1 =  -------  n
 t n + K   t n – K  1 –    (1-95)
n+1
 =   t n + t n  (1-91) + F
n+1 n
 + F 1 – 

at the new time plane.

NKSKPKO mêÉÇáÅíçêJÅçêêÉÅíçê=ãÉíÜçÇ
NKSKPKN  JjÉíÜçÇ The predictor-corrector method is thoroughly
described elsewhere3,4,5,6. For the present analysis the
Introducing a weighting coefficient  0    1  , we fully implicit backward Euler (BE) scheme with a first-
can write order accuracy and the semi-implicit nondissipative
trapezoid rule (TR) with a second-order accuracy are
  t n + t n  =   t n + t n  +  1 –    t n  
enforced. The time derivatives are approximated, for
F  t n + t n  = F  t n + t n  +  1 –  F  t n   (1-92) the BE scheme, by
· · · 
  t n + t n  =   t n + t n  +  1 –    t n  
n+1 n
·n+1  –
 = -------------------------- (1-96)
t n
Using a backward difference approximation for
· ·
  t n + t n  and a forward difference for   t n  one
obtains and for the TR scheme, by

n+1 n ·n+1 2 n+1 n ·n

·  –  = --------   – – . (1-97)
  t n + t n  = -------------------------- (1-93) t n
t n

Inserting (1-96) and (1-97) into (1-87) results in

Common time stepping schemes result if choosing  in
an appropriate manner, viz., n
O-
 ----------  1 ·n
+ K  = O ----------- +  --- – 1  + F
n+1 n+1
(1-98)
 t n  t n   
 = 0 explicit scheme 

 = 12 trapezoid rule (Crank-Nicolson scheme)  (1-94)
 with    1--2- 1  for the TR and BE scheme, respec-
 = 1 implicit scheme 
tively.

Inserting (1-92) into (1-87) the following matrix equa-

tion finally results

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=ON
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
NKSKQ cáåáíÉJÉäÉãÉåí=Ä~ëáë=çéÉê~íáçåë alized (local) coordinates (see Fig. 1.6). The coordinate
transformation (or mapping) that bridges a computa-
D
A fundamental aspect of the finite-element method tional (transform)  -space and the Euclidean space R
is the use of master elements where all element-data is
inner products and integrations are performed in gener-

1 x
R
(+1)
(-1)


e e
e :   x = x   

(-1,1) (1,1)
x
2 y
R


(-1,-1) (1,-1)


(-1,-1,1) (-1,1,1)

3
R (1,1,1)
(1,-1,1) x

(-1,-1,-1)  z
(-1,1,-1)

(1,-1,-1) (1,1,-1) y

D
Figure 1.6 Finite elements with one-to-one mapping onto R ( D = 1,2,3).

e :   x = x    Based on (1-99) it is convenient to express the basis


 functions N i in local  -coordinates for each element e,
 –1    1  viz.,
 (1-99)
 =  –1    1 

 –1    1 

OO=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

e
Ni = Ni  x  =  Ni    (1-100)
with x =  x 1 x 2 x 3  according to (1-11). To evaluate
e
the flux vector divergence terms in (1-88) the inverse
x = Ni xi Jacobian is required

The mapping  e is one-to-one and onto its range pro-  N

--------i 
vided the transformation Jacobian J is nonsingular,   
3  
where in the R space – 1  N i 
N i = J  --------  (1-102)
  
 N 
   --------i 
     
  J 11 J 12 J 13
x  
J = ------ =   x 1 x 2 x 3  = J 21 J 22 J 23 (1-101)
   where
  J 31 J 32 J 33
 
  

x x x 3
-------1- -------2- --------
  
x
= -------1- x 2 x 3
-------- --------
  
x 1 x 2 x 3
-------- -------- --------
  


  J 22 J 33 – J 32 J 23   J 13 J 32 – J 12 J 33   J 12 J 23 – J 13 J 22 
 1- 3
 ------  J 31 J 23 – J 21 J 33   J 11 J 33 – J 13 J 31   J 21 J 13 – J 23 J 11  in R
 J
  J 21 J 32 – J 31 J 22   J 12 J 31 – J 32 J 11   J 11 J 22 – J 12 J 21 
–1  
J = ------ =  (1-103)
x  1 - J 22 – J 12 2
 ------ in R
 J –J J
21 11

 1- 1
 ------ in R
 J

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=OP
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
with the determinant of J

 3
 J 11  J 22 J 33 – J 32 J 23  – J 21  J 12 J 33 – J 13 J 32  + J 31  J 12 J 23 – J 13 J 22  in R
 2
J =  J J –J J in R (1-104)
 11 22 21 12
 J in R
1
 11

The master element matrices appearing in (1-87) and

e
(1-88) are to be integrated over element volumes 
e
and surfaces  . The integration in local coordinates
becomes for a ’volume’ element

 dx d y d z = J d d d 
  D
 dxdy = J d d  Cartesian R  D = 1 2 3 
 
d =  dx = J d  (1-105)



 rdrddz = 2 J r d d axisymmetric

and for an ’areal’ element in Cartesian coordinates of

D
R  D = 1 2 3  space:

OQ=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

 x- x
 ----- -------
   i j k J 12 J 23 – J 13 J 22
 y- d d = det J J J d d =
y-  ------
 ----- 11 12 13 J 13 J 21 – J 11 J 23 d d at  = 1
  
 J 21 J 22 J 23 J 11 J 22 – J 12 J 21
z
------
z
-------

  

 x x
 ------- ------
   i j k J 22 J 33 – J 23 J 32
 y y
 -------  ------ d d = det J 21 J 22 J 23 d d = J 23 J 31 – J 21 J 33 d d at  = 1
  
 J 31 J 32 J 33 J 21 J 32 – J 22 J 31
 z z-
------- -----
  

 x x
 ------ ------
   i j k J 12 J 33 – J 13 J 32

 y y at  = 1
------  ------ d d = det J 11 J 12 J 13 d d = J 13 J 31 – J 11 J 33 d d
  
 J 31 J 32 J 33 J 11 J 32 – J 12 J 31
d =  z
------
z
------
(1-106)
  



 x-
 -----
  d = J 11
 d =
2 2
J 11 + J 12 d at  = 1
 y J 12
------
 


 x
-------
  d = J 21
 d =
2 2
J 21 + J 22 d at  = 1
 y- J 22
------
 



  = –1

 =1



«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=OR
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
2
and in cylindrical coordinates of R (meridional) space

 r
 ------
  r d d = 2 J 11 rd = 2 J 2 + J 2 rd at  = 1
 11 12
z- J 12
 -----
 
d =  (1-107)
 r
 ------
  r d d = 2 J 21 rd = 2 J 2 + J 2 rd at  = 1
 z J 22
21 22
 ------
 

where r = Ni   r i . Commonly, for 2D and 3D element a Gaussian quadrature rule, e.g.,

entities the volume and area integrals are evaluated via

1 1 1 n Gauss n Gauss n Gauss

 f  x  d =   f    d =     f      d d d =     wi wj wk f  i j k 
 e e e –1 –1 –1 e i=1 j=1k=1
(1-108)
1 1 n Gauss n Gauss

 g  x  d =   g    d =    g     d d =    wi wj g  i j 
 e e e –1 –1 e i=1 j=1

where the n Gauss is the number of Gauss points, w i are NKSKR ^ëëÉãÄäó= çÑ= íÜÉ= ÇáÑÑÉêÉåí= ÑÉ~J
weighting coefficients and the indices  i j k  indicate íìêÉ= ÉäÉãÉåíë= íç= íÜÉ= ÖäçÄ~ä= ëóëíÉã
the positions of the evaluation points in their local
ã~íêáñ
coordinates  . The functions f  .  and g  .  in the inte-
grands are marked by an asterisk if the volume and sur-
face integrals are expressed in the  -coordinates NKSKRKN kÉÉÇë=Ñçê=ÅççêÇáå~íÉ=íê~åëÑçêã~J
according to (1-105), and .
íáçå
For 1D elements the integrals of (1-88) can easily
The global matrix equations (1-87) written in the
evaluated in a direct analytical manner so as shown in
form
Appendix B for a channel element with a linear basis
function N .

OS=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

· conflicts. We take into consideration that all flow and

O+K –F = 0
 transport processes are invariant with respect to a rota-
O = O
e  tion (orthogonal transformation) of the global coordi-

e  nates x . Accordingly, we can arbitrarily rotate x to the
 x -coordinates by using a suitable matrix of directional
e  (1-109)
K = K  cosines a as
e 

F = F
e  x = a  x

e  a 11 a 12 a 13 x
x (1-110)
y = a 21 a 22 a 23 y
represent the standard discrete system resulting from z a 31 a 32 a 33 z
the summation (assembly) of the elemental (e) matrix
contributions. The integrals of the matrices and vectors
e e e
O , K and F for each element e are performed in the Taking an appropriate rotation of the global x – y – z -
local coordinates  for the corresponding Euclidean coordinate system in such a way that the resulting local
D
space R as stated above. x – y – z -system becomes aligned to the orientation
3
of the 2D or 1D elements in the R space, there will be
Under normal conditions 1D finite elements are no more an elemental contribution to the z -direction
1 2
mapped to the R space, 2D elements to the R space for 2D elements and elemental contributions to the y -
3
and 3D elements to the R space. In such case the map- and z -directions for 1D elements (see Fig. 1.7).
ping is strictly one-to-one, that means 3 global coordi-
nates x y z are transformed to 3 local coordinates The advantages of this coordinate transformation
   in 3D, 2 global coordinates x y to 2 local coor- are that the corresponding Jacobian J
dinates   in 2D and 1 global coordinate x to 1 local
coordinate  in 1D. However, when 1D and 2D dis- x
J = -------- (1-111)
crete feature elements are generally mapped onto a 3D 
global space, the number of local coordinates  will be
less than the number of global coordinates x and the becomes again an invertible square matrix and the stan-
Jacobian J (1-101) will not be any more an invertible dard metric procedure can be maintained in the assem-
square matrix (e.g., for the  –  -system of a 2D fea- bly process for the global matrix system (1-109). To
ture element mapped into the global x – y – z -system ease the computations the x – y – z -coordinate system
the third row of J contains zeros, may, in fact, be different for every element e.
J 31 = J 32 = J 33 = 0 , because the  -coordinate does
not exist in 2D elements).

There is a simple way to overcome this mapping

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=OT
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt

e e
 e :   x = x   
3
z 2D element in the R space
y 
(-1,1) (1,1)

y
z x


(-1,-1) (1,-1)
y
x

3
1D element in the R space
y
z

x

z (+1)
(-1)

x

Figure 1.7 Global x – y – z -coordinate system, rotated elemental x – y – z -coordinate system and local  –    -coor-
3
dinate system for 2D and 1D elements in the R space.

NKSKRKO dÉåÉê~äáòÉÇ=ÇáëÅêÉíÉ=ëóëíÉã 
e Te
x = a  x x=a  x
Since the rotation matrix a forms an orthonormal 
e e
 = T  
e e Te
 = T  
e  (1-112)
basis we can transform between the x - and the x -sys- 
F = T  F 
e e e e Te e
tem according to F = T  F

OU=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

where e
K = T
Te e
  K  T 
e
(1-115b)

e e
a
e
0 0 
where q and q represent elemental ’flux’ compo-
e e
nents for the local and the global coordinate system,
T = 0 a 0 (1-113) respectively. On inserting (1-112), (1-115a) and (1-
e
0 0 a 115b), using the same transformation rules for the stor-
 age matrix O and assuming that in general the right-
hand side vector F can also be directionally dependent,
e the global matrix system (1-109) yields finally the form
is a diagonal directional cosines matrix built up of a
matrices in a number equal to that of the nodes in the ·
O+K–F = 0 
element e 
O =   T   O  T   
Te e e
and
e 

Te e e  (1-116)
Te e e e K =   T   K  T   
a 0 0  a 11 a 21 a 31
Te e 
Te Te 
T = 0 a 0 a = ae ae ae (1-114)
12 22 32 Te e
F =   T  F  
0 0 a
Te
e e e 
a 13 a 23 a 33 e 

e
e e
where the rotation matrix T is evaluated at element
are the transposes of matrices T and a , respectively. e
level e. Practically, T is only required for mapping 2D
3
and 1D feature elements in the general R space, while
We can usually assume that any directional proper- 3D elements need not to rotate to the x -system (the
ties of the discrete system (1-109) are available in the e
rotation matrix becomes unity T  I ) and can be
local x -coordinates. Then, the local (elemental) sys- directly mapped onto the local  -space via the Jaco-
e e
tem matrix K = K x is transformed into the global bian J (1-101).
e
matrix K according to
However, there are important special cases to be
e e e e e e
q = K   = K   T    considered here. Firstly, if the material properties of the
(1-115a) e e
e
q = T
Te e
 q =  T
Te e e
  K  T     = K  
e e e square matrices O and K are independent of the
coordinate directions (isotropic conditions) then we
have
with

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=OV
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
T
Te e e
  O  T   = O
e tors u 1 and u 2 (Fig. 1.8), which are parallel to the
Te e e e
(1-117) local  - and  -axes, respectively. They can be found
T   K  T   = K by the following shape-derived relationships
Te e
because T  T = I due to the condition of orthogo-  x-
e e e e  -----
nality and accordingly O = O , K = K .  J 11

 y = at  = 1
 ------ J 12
Secondly, if the sink/source and boundary condition  
terms incorporated in the right-hand side vector F rep-  z- J 13
 -----
resent direction-independent quantities so as occurred  
in the balance equations summarized in Tables 1.5, 1.6 
 x
and 1.7, the vector F consists of nodal scalars and can  ------
be directly evaluated (no rotation), viz.,   J 21
 y = J
u1 =  ------ 22
at  = 1 (1-119)
e  
F =  F (1-118) 
 z
------
J 23
e
 

 x-
NKSKRKP aÉíÉêãáå~íáçå=çÑ=íÜÉ=ÇáêÉÅíáçå~ä=  -----
  J 31
e 
ÅçëáåÉë= a =çÑ=ÉäÉãÉåí=É  y
------ = J 32 at  = 1
 
 J 33
e
The directional cosines a are only required for  z-
-----
3  
mapping 2D and 1D discrete feature elements in the R
space. If we commonly assume that the 3D continuum
domain  with its boundary  is completely filled
by 3D finite elements (e.g., hexahedral or pentahedral
isoparametric elements), the 2D fracture and 1D chan-
nel elements share the nodal points of the 3D mesh and
their geometric extents are aligned to surfaces, edges or
diagonals of the 3D matrix elements (Fig. 1.8).

For 2D fracture elements forming sur-

faces of the 3D matrix element it is convenient to
derive the directional cosines directly from the shape of
the 3D element. We can construct the 2 directional vec-

PM=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

y  y 
z 
u2

x 
u1 x 
u2
u1

z 

2D fracture element 2D fracture element

at top surface at a side surface
z z
y y
x 
n x
x x

z 

u1 u1

y 

1D channel element 1D channel element

at an edge m through the diagonal

Figure 1.8 Exemplified mapping of 2D fracture elements and 1D channel elements aligned to surfaces, edges and diago-
nals, respectively, for a 3D finite matrix element. Local and global coordinates.

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=PN
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
Notice, for 2D feature elements (fracs) we need
 x
 ------ only two directional vectors (i = 1,2), the remaining
  J 21 e
directional cosines a 3j are meaningless.
 y = J at  = 1
 ------ 22
 
 J Often we can assume that the 2D feature elements
z 23
 ------ are perfectly flat, i.e., they represent noncurved 2D
 
 geometries which occur for arbitrarily oriented linear
 x- triangles or for vertical linear quadrilaterals in the 3D
 -----
  J 31 space. Instead of using the above shape-derived expres-
 y sions (1-119) and (1-120), in such cases it is convenient
u2 =  ------ = J 32 at  = 1 (1-120)
  to derive the directional vectors u i in a direct manner
 z- J 33 (see Fig. 1.9).
 -----
 
 We specify the x -axis along the edge nm of the 2D
 x-
 ----- feature element. The vector u 1 is accordingly given by
  J 11

 y
------ = J 12 at  = 1
  xn – xm
 J 13 u1 = yn – ym (1-123)
 z-
-----
  zn – zm

These directional vectors can be easily used to compute

The second directional vector u 2 derived by simple
the directional cosines according to
vector algebra (as summarized in the Appendix C)
T yields
e ui  ej i = 1 2
a ij = cos  u i e j  = -------------------------- for (1-121)
ui ej j = 1 2 3 q  u1
u 2 = q –  ---------------- u 1 (1-124a)
 u 1  u 1




= 1

with the auxiliary vector q formed along the adjacent

with the basis (unit) vectors side lm of the 2D element as

1 0 0 xl – xm
e1 = 0 e2 = 1 e3 = 0 (1-122)
q = yl – ym (1-124b)
0 0 1
zl – zm

PO=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

k
z l
n
l

y
z z
u2
x n
y
m u1 y
x

u1
u2
m

x
Figure 1.9 Directional vectors u i (i = 1,2) for a linear triangular element and a vertically-oriented linear quad-
rilateral element.

e 
xn – xm
and the directional cosines a ij (i = 1,2; j = 1,2,3) can be  x
easily computed by using (1-121). u1 = yn – ym  = y

 z (1-125a)
zn – zm 
For 1D channel elements the same procedure 
e
can be applied to determine a ij (for i = 1; j = 1,2,3). u 1 = x + y + z 
2 2 2
e
Here, only one row a 1j of the rotation matrix is of
interest. Taking into consideration that 1D feature ele-
x
a 11 = -------------------------------------------- 
e
ments can be rather arbitrarily placed at mesh nodes
2 2 2
(which are not necessarily connected in one element x + y + z 
and oriented along edges) the following direct evalua- 
y
a 12 = -------------------------------------------- 
e
e
tion procedure can be used to compute a 1j for a 1D lin- 2 2 2
(1-125b)
ear channel element spanning between the two nodes n x + y + z 

z
a 13 = -------------------------------------------- 
and m (cf. Fig. 1.8): e
2 2 2
x + y + z 

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=PP
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
NKSKRKQ aÉãçåëíê~íáîÉ=Éñ~ãéäÉ

To demonstrate the assembly process for the differ-

ent feature elements let us consider the simple model as
shown in Fig. 1.10 consisting of only one 3D element
(e = 1) connected with both a 2D fracture element (e =
2) and a 1D channel element (e = 3).

) x y z
=2 Node
(e [unit] [unit] [unit]
ment
e ele
tur 5 1 1. 1. 0.
frac
2 D
2 11. 1. 0.

)
=3
3 7. 11. 0.
4

(e
4 1. 1. 7.

ent
em
3D matrix element (e = 1) 5 11. 1. 9.
l el
ne
6 7. 11. 11.
han
3
c
1D

z
y
1 2
x
Figure 1.10 3D prismatic matrix element connected with a 2D fracture element on the top surface and a 1D
channel element through a diagonal.

The assembly of the resulting matrix system has to be

performed according to (1-116), viz.,

PQ=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

·
O+K–F = 0 
T1 1 1 T2 2 2  T3 3 3
O =  T   O  T   +  T   O  T   +  T   O  T   
 (1-126)
T1 1 1 T2 2 2 T3 3 3
K =  T   K  T   +  T   K  T   +  T   K  T   

T1 1 T2 2 T3
F =  T  F  +  T  F  +  T  F 
3 

For the 3D triangular prism (e = 1) it is assumed that for coordinate transformation into the x – y – z -coor-
the elemental material properties can be specified in dinates. Accordingly,
the global coordinates ( x – y – z ). Here, there is no need

1 00 
1 T1 0 10 1 1 1 1 1 1
T = T = O = O K = K F = F (1-127)
0 01


For the 2D triangular fracture element (e = 2) the direc-

tional vectors are given by

 10 
 
 6 10 4  0 
10 6 6  10 – 0.5385
2 
 ------------------------------------
u1 = 0 q = 10 u 2 = 10 – 0 = 10 (1-128)
 
2 4 4  10  2 2.6923
 10 0 2  0 
 
 2 










t = 0.6538

where the formula (1-124a) is used. We note that

u 1 = 10.198 and u 2 = 10.370 . It leads to the follow-
ing directional cosines (1-121)

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=PR
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt

0.9806 0 0.1961 
0.9806 0 0.1961
2 2 – 0.0519 0.9843 0.2596
a = – 0.0519 0.9843 0.2596 T = (1-129)
. . .
. . .


and we get the matrices to be assembled for the ele-

ment (e = 2) as


  
  2 2 
 0.9806 – 0.0519 .   O xx O xy 0  0.9806 0 0.1961  
  
0 0.9843 . 2 2
 O yx O yy 0 – 0.0519 0.9843 0.2596  if anisotropic
T2 2 2 
 T   O  T   =  0.1961 0.2596 .  . . . 
  0 0 0 
    
   

 O2 if isotropic

(1-130)
T2 2 2
similar for  T   K  T  

 2
F x
 0.9806 – 0.0519 . 
 0 0.9843 . 2
T2 2  F y if directional
 T  F  =  0.1961 0.2596 .
 0
 
 
 2
F if scalar

4
Taking the u 1 -vector for the 1D channel element u 1 = – 10 u 1 = 14.036 (1-131)
(e = 3)
9

PS=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

we obtain

0.2850 – 0.7125 0.6412 

0.2850 – 0.7125 0.6412
3 3 . . .
a = . . . T = (1-132)
. . .
. . .


and finally the matrices for the element (e = 3) as

  
  O3 0 0 
 0.2850 . .   xx
0.2850 – 0.7125 0.6412  
 – 0.7125 . .  . . . 
T3 3 3   0 0 0  if anisotropic
 T   O  T   =  0.6412 . .  . . . 
0 0 0
  
   
 

 3
O if isotropic
T3 3 3 (1-133)
similar for  T   K  T  

 3
 0.2850 . .  F x
 – 0.7125 . .
T3 3  0 if directional
 T  F  =  0.6412 . .
0

  

 3
F if scalar

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=PT
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
oÉÑÉêÉåÅÉë ^ééÉåÇáñ=^
1. Chandhry, M. H. & M. E. Barber, Open channel flow. In: The kçãÉåÅä~íìêÉ
handbook of fluid dynamics (ed. by R. W. Johnson), CRC Press/
Springer, Chapter 45, 1998.
2. Diersch, H.-J. G., FEFLOW - Reference Manual. WASY Ltd., In the above the symbols have the following meaning:
Berlin, 2002.
3. Diersch, H.-J. G., Modeling and numerical simulation of geohy- Latin symbols
drodynamic transport processes (in German). Reprint, WASY
Ltd. Berlin, 1991. 2
A = flow area,  L  ;
4. Diersch, H.-J., Finite element modelling of recirculating density-
driven saltwater intrusion processes in groundwater, Advances in a = coefficient to specify boundary
Water Resources 11 (1988) 1, 25-43. conditions;
5. Diersch, H.-J. G. & O. Kolditz, Coupled groundwater flow and a = rotation matrix;
transport: 2. Thermohaline and 3D convection systems, Advances B = thickness or depth,  L  ;
in Water Resources 21 (1998) 5, 401-425.
b = aperture or surface elevation,  L  ;
6. Diersch, H.-J. G. & P. Perrochet, On the primary variable switch-
ing technique for simulating unsaturated-saturated flows, b = coefficient to specify boundary
Advances in Warer Resources 23 (1999) 3, 271-301. conditions;
–3
7. Gottardi, G. & M. Venutelli, LANDFLOW: Computer program C = concentration,  ML  ;
for the numerical simulation of two-dimensional overland flow. C = Chezy roughness coefficient,
Computers & Geosciences 23 (1997)1, 77-89. 1  2 –1
L T  ;
8. Gray, W.G., Derivation of vertically averaged equations describ- f s
ing multiphase flow in porous media. Water Resources Research cc = specific heat capacity of fluid and
2 –2 –1
18 (1982)6, 1705-1712. solid, respectively,  L T   ;
9. Di Giammarco, P., E. Todini & P. Lamberti, A conservative finite D = space dimension, (1, 2 or 3);
elements approach to overland flow: the control volume finite D = tensor of hydrodynamic dispersion,
element formulation. Institute for Hydraulic Construction, Univ. 2 –1
L T  ;
Bologna, Italy. 2 –1
10. Panton, R. L., Incompressible flow. 2nd Edition, J. Wiley & Sons,
Dd = molecular diffusion,  L T  ;
New York, 1996. Dm = tensor of mechanical dispersion,
2 –1
11. Taylor, G., Dispersion of solute matter in solvent flowing slowly L T  ;
–1
through a tube. Proc. R. Soc. London A, 219 (1953), 186-203. d = strain-rate tensor of fluid,  T  ;
E = internal (thermal) energy density,
2 –2
L T  ;
e = gravitational unit vector,  1  ;
ei = basis vectors,  1  ;
F = material surface;
f = specific rate of temporary
production;

PU=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

f. = general function;  = retardation,  1  ;

f = viscosity relation function,  1  ; d = derivative retardation,  1  ;
–2
g = gravitational acceleration,  LT  ; r = radius,  L  ;
g. = general function; rc = homogeneous chemical reaction
–2 –3 –1
g = gravity vector,  LT  ; rate,  ML T  ;
h = hydraulic head,  L  ; r hydr = hydraulic radius,  L  ;
I = unit (identity) tensor,  1  ; S =  BS o + S s  , storage term,  1  ;
J = Jacobian matrix; So = bed slope, inclination of the bottom
j = flux vector; plane to the horizontal x and y
directions,  1  ;
jc = Fickian mass flux vector, –1
–2 –1
 ML T  ; So = compressibility,  L  ;
jT = Fourierian heat flux vector, Sf = vector of friction slopes at channel
–3
 MT  ; bottom,  1  ;
j = surface or interface exchange, Ss = storativity,  1  ;
–1 –2 T T o = temperature and reference
 ML T  ;
K =  k o g    o , tensor of hydraulic temperature, respectively,    ;
–1 Tt = final time,  T  ;
conductivity,  LT  ;
K = tensor function specifying different t = time,  T  ;
–1
laws of flow motion,  LT  ; x = coordinate vector,  L  ;
2
k = tensor of permeability,  L  ; x y = Cartesian coordinates,  L  ;
M = number of nodes; z = axial or vertical coordinate,  L  ;
–1
M = Manning roughness coefficient, v = velocity vector of fluid,  LT  ;
1  3 –1
L T  ; w = velocity vector of surface or
–1
N = basis function; interface,  LT  ;
n = outward-directed normal unit vector w = spatial weighting function;
at surface;
n Gauss = number of Gauss points in each local
coordinate direction; Greek symbols
–1 –2
p = fluid pressure,  ML T  ;
–3 –1
Qc = mass source term,  ML T   = constant of friction slope
–1 –3
QT = source/sink of heat,  ML T  ; relationship,  1  ;
–1
Q = fluid mass sink/source,  T  ;  L  T = longitudinal and transverse
q = flux vector; dispersivity, respectively,  L  ;
R = radius of circular tube,  L  ; i = portions i of boundary  ;
D –1
R = space of dimension D ;  = fluid compressibility,  L  ;

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=PV
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
t n = time step length at time plane n,  = generalized friction factor;
T ;  = Newton-Taylor roughness
 = total boundary; coefficient,  1  ;
2
S = projected area of surface,  L  ;  = pressure head or local water depth,
3 L ;
V = volume of REV,  L  ;
 = porosity (= volume fraction of fluid C  = sorption function,  1  ;
phase),  1  ;  = balance quantity;
 =  – 1    1  , local coordinate,  1  ;  = domain;
 =  – 1    1  , local coordinate,  = azimuthal angle,    ;
1 ; –1
 = Nabla (vector) operator,  L  ;
 = local coordinate vector,  1  ;
 =   –  o    o , density ratio or Subscripts
buoyancy coefficient,  1  ;
 = weighting coefficient,  0    1  ; corr corrected;
–1
 = chemical decay rate,  T  ; e elemental;
 = coefficient of skeleton compressibi- n time plane;
–1
lity,  L  ;
i j k nodal or spatial indices;
 = tensor of thermal hydrodynamic
–3 –1 o reference value;
dispersion,  MLT   ;
f s
 = thermal conductivity for fluid and
–3 –1
solid, respectively,  MLT   ; Superscripts
  o = dynamic viscosity and reference
dynamic viscosity of fluid,  phase index;
–1 –1
respectively,  ML T  ; D space dimension;
 =  – 1    1  , local coordinate,  1  ; e elemental;
  o = fluid density and reference fluid f fluid (water) phase;
–3
density, respectively,  ML  ; n time plane;
 = viscous stress tensor of fluid, s solid phase;
–1 –2
 ML T  ;
T transpose of a matrix;
' = deviatory stress tensor of fluid,
–1 –2  surface index;
 ML T  ;
bottom
 shear stress at bottom surface,
–2 –2
 ML T  ;
interface –2 –2
 interfacial shear stress,  ML T  ;
top
 = shear stress at top surface,
–2 –2
 ML T  ;

QM=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

^ééÉåÇáñ=_
Furthermore, we have for the element, cf. (1-100),
^å~äóíáÅ=Éî~äì~íáçå=çÑ=ã~íêáñ= ÉäÉãÉåíë
x = N1 x 1  + N2 x 2  (B3)
ENJUUF=Ñçê=~=Na=ÅÜ~ååÉä=ÉäÉãÉåí

We consider the following linear 2-node element e and with and (B4)

x N 1 N 2 x e
J = J 11 = ------ = --------- x  1  + --------- x  2  = -------- (B4)
Ni 1    2

and
N1 N2
–1 1 2
J = ------- = -------- (B5)
J x e

0 Then, the divergence terms (1-102) becomes with (B2)

1 e 2
N 1 1
 = -1   = +1 x --------- – --------
x = x(1) x = x(2) N 1 –1  = x e
xe = J (B6)
N 2 N 2 1
--------- --------
 x e

with the basis functions at the nodes 1 and 2

According to (1-105)
1
N 1 = ---  1 –   x e
2 d = dx = J d = -------- d (B7)
(B1) 2
1
N 2 = ---  1 +  
2
and its derivatives the matrices (1-88) become for element e

N
---------1 = – 1---
 2
(B2)
N 2 1
--------- = ---
 2

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=QN
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt

1 e
N1 N1 N1 N2 e 2 2 x S x e
d = -----   1 –    1 –   -------- d = -------------- 2 1
e e S e
O =  S
N2 N1 N2 N2 4 2 2 2 6 12
(B8)
e –1  1 –    1 +  

e e e e e
K = K1 + K2 + K3 + K4
1
N 1 N 1 N 1 N 2 e x e e
d = ------------  –  1 –    1 –   -------- d = ----- – 1 1
e e q q
K1 =  q
N 2 N 1 N 2 N 2 2x e –  1 +    1 +   2 2 –1 1
e –1
1
N 1 N 1 N 1 N 2
D
e x e D
e
d = --------2-  1 – 1 -------- d = -------- 1 – 1
e e
K2 =  D
N 2 N 1 N 2 N 2 x e – 1 1 2 x e – 1 1
e –1 (B9)
1 e
N1 N1 N1 N2 e 2 2 x  x e

d = ------   1 –    1 –   -------- d = ---------------- 2 1
e e e
K3 =  N2 N1 N2 N2 4 2 2 2 6 1 2
e –1  1 –    1 +  
 = 1 = – 1
e 2 2
N1 N1 N1 N2 a  1 – 1   1 – 1 
= a 1 0
e e e
K4 = a = -----
N2 N1 N2 N2 4 2 2 0 1
 1 – 2   1 + 2 
 = 2 = 1

e e e
F = F1 + F2
1 e
N1 e x e Q x e
d = ------   1 –   -------- d = --------------- 1
e e Q
F1 =  Q
N2 2 1 +  2 2 1
e –1 (B10)
 = 1 = –1
e e
N1  a 2 – b  1 – 1
=  a 2 – b  1
e e e e e
F2 =  a 2 – b  = ----------------------------
N2 2 1 + 2 1
 = 2 = 1

QO=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

Finally, the discretized matrix equation (1-88) can be

summarized as

e  · e e  e
 S x e 2 1   1   q e – 1 1 D
e  x e 2 1 
e 1 0   1
1 – 1
  -------------
6
-    +  -----
1 2  e 2
·   –1 1
+ --------
x e –1 1
+ ----------------
6 1 2
+a   –
0 1   e 
e
 2  2 (B11)
e
Q x e 
– --------------- 1 –  a  2 – b  1  =  0 
e e
2 1 1

^ééÉåÇáñ=`
y
aÉêáî~íáçå= çÑ= íÜÉ= çêíÜçÖçå~ä= ÇáêÉÅJ l
íáçå~ä=îÉÅíçêë= u i =Ñçê=~=Oa=ÛÑä~íÛ=ÉäÉãÉåí n
x
u2
q q–v p
u1
For a linear 2D ’flat’ element we can find the fol-
z
lowing vectors for a typical element with the nodes y v
mnl  k  : m

x
xn – xm xl – xm
p  u1 = yn – ym q = yl – ym (C1) the parameter t can be easily found if inserting (C3)
zn – zm zl – zm into (C2)

 q – tp   p = 0 
Since the vector u 2  q – v is perpendicular to the vec- 
qp (C4)
tor p  u 1 the dot product yields t = ---------- 
pp
q – v  p = 0 (C2)
As a result, we get
Using the parametric description of the vector v as
qp
v =  ---------- p (C5)
v = tp (C3) pp

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=QP
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
and finally (with p  u 1 ) On the other hand, the friction laws for the Hagen-
Poiseuille flow in form of Eq. (1-56) and for the Man-
qp  ning-Strickler flow in form of Eq. (1-67) (with Table
u 2  q – v = q –  ---------- p  1.4) are specified by the hydraulic aperture b and the
pp 
 (C6) friction parameter M , respectively (cf. Table 1.9).
q  u
u 2 = q –  ---------------- u 1 
1
 u 1  u 1 Table 1.9 Frictional input parameter

Law Parameter

Hagen-Poiseuille hydraulic aperture, b

^ééÉåÇáñ=a
Manning-Strickler roughness, M
fåéìí= çÑ= ÅêçëëJëÉÅíáçå~ä= Ç~í~I= ëí~åÇ~êÇ
áãéäÉãÉåí~íáçåë=çÑ=ÜóÇê~ìäáÅ=ê~Çáá=~åÇ Based on the input parameters b , A or B in
íÜÉáê= êÉä~íáçåë= íç= ÇáÑÑÉêÉåí= íóéÉë= çÑ FEFLOW the following relationships of the hydraulic
radius r hydr are implemented according to the dimen-
Ñê~ÅíìêÉë= Ñçê= íÜÉ= e~ÖÉåJmçáëÉìáääÉ= ~åÇ
sion of the fracture elements. They represent so-called
íÜÉ=j~ååáåÖJpíêáÅâäÉê=ÑêáÅíáçå=ä~ïë standard implementations as summarized in Table 1.10.

In FEFLOW a reduced data set is used to input the Table 1.10 Implemented standard hydraulic radii
geometric relationships of the fractures. There is no
direct input of the hydraulic radii r hydr for the different r hydr
types of fracture elements applied to both the Hagen-
Poiseuille and the Manning-Strickler law of fluid Dimension expressed by appropriate input
motion. Instead, the geometry of the cross-sections is of fracture parameters
input by the flow area A for 1D elements and by the elements
thickness/flow depth B for 2D fracture elements (see Hagen- Manning-
Table 1.8). Poiseuille Strickler

Table 1.8 Geometric input parameter 1D A4

(plane) (submerged qua-
Dimension Fracture type/case Parameter b2 dratic cross-sec-
(Table 1.5) (submerged slit tion)
plane, type B of
2D Table 1.3) B2
1D 1DPP, 1DPN, 1DAP, flow area, A
(plane) (submerged slit
1DAN
plane, type B of
2D 2DPP, 2DPN thickness, B Table 1.3)

QQ=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë
3 –3 –3
 o = 10 kg m ,  o = 1.3  10 Pa s , and
–2
As seen from Table 1.10 the input parameter in form g = 9.81 ms . It results a factor of
6 –1 –1
of the hydraulic aperture b for the Hagen-Poiseuille f o =  o g   o = 7.55  10 m s . A hydraulic radius
law and the geometric input parameter in form of the which is different from the standard geometry, and
flow area A or thickness B for the Manning-Strickler parameters, which are different from the standard
law are used to express r hydr internally in FEFLOW. parameter factor f o , can be derived from the identity
Accordingly, FEFLOW assumes that in the case of
using a Hagen-Poiseuille law a submerged slit geome- 2 2
r hydr b- g
try is standard (both in 1D and 2D). In case of the Man- ----------- f = ----- f with f = ------ (D2)
a 12 o 
ning-.Strickler law a submerged quadratic cross section
2
( A = B ) is considered the standard geometry in 1D
A corrected hydraulic aperture b corr can be obtained
and a submerged slit plane geometry is standard in 2D.
from (D2) as
The question arises how is it possible to differ from the
standard geometry types? For example, instead of a
2 3 3 for plane geometry
quadratic cross-section for a 1D fracture element, a b corr = ---f- ---------- r hydr a =  (D3)
rectangular plane or an axisymmetric geometry is to be fo a 2 for axisymmetric geometry

used.
where r hydr is the actual (true) hydraulic radius, which
To specify hydraulic radii r hydr which are different can be taken from Table 1.3, and f is the true parameter
to the embodied standard geometries (Table 1.10) one factor, where dynamic viscosity, gravity and density
can input corrections in the hydraulic aperture b for can be specified different from the standard settings in
the Hagen-Poiseuille law and in the Manning rough- f o . Table 1.11 summarizes the corrected hydraulic
ness coefficient M for the Manning-Strickler law. apertures b corr for interesting applications.
These corrected parameters can be derived in the fol-
lowing manner.

The standard hydraulic conductivity K for the

Hagen-Poiseuille law is according to Eq. (1-56) and
Table 1.10:

2
r hydr  o g b o g
2
K = ----------- --------- I = ------ --------- I (D1)
a o 12  o

where a = 3 for plane geometry. Furthermore, the fol-

lowing standard parameters are used here for water:

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=QR
aáëÅêÉíÉ= ÑÉ~íìêÉ= ãçÇÉäáåÖ= çÑ= ÑäçïI= ã~ëë= ~åÇ= ÜÉ~í= íê~åëéçêí= éêçÅÉëëÉë= Äó= ìëáåÖ
cbcilt
Table 1.11 Corrected apertures b corr for different On the other hand, the standard hydraulic conduc-
applications in the case of Hagen-Poiseuille law for tivity K for the Manning-Strickler law is according to
1D and 2D fracture elements a) Eq. (1-67) with Tables 1.4 and 1.10:

Type b corr A 2  3 ------------------

 M  ------- I -
for 1D
  4 4 2
submerged rectangular cross- 23 I  h
section K = Mr hydr ------------------- =  (D4)
A 4 2  B 23 I -
h 
 M  --2- ------------------ for 2D
Bb f  4
h
2
B ----------------- ----
 b + B  fo
b
Accordingly, from (D4) we can find a corrected Man-
submerged slit plane ning coefficient M corr in the following form to input a
B standard-different geometry of the fracture elements in
b
b ---f- 1D and 2D
fo
B
no correction is
needed if f = fo  4  2  3
  ------- for 1D
open rectangular cross-section 23  A
M corr = Mr hydr   (D5)
C   --2- 2  3
  B for 2D
B 2Bb - f 
--------------- ----
b + 2B f o
b
where M is the true (physical) Manning roughness
open wide channel (b > 20B) coefficient. Tables 1.12 and 1.13 summarize the correc-
D plane
f tions M corr  M for the Manning coefficients for 1D and
B 2B ---o- 2D fracture elements, respectively, which are required
f
b when fracture cross-geometries different to the stan-
dard geometry have to be input.
submerged circular cross-section
E ------3- R ---f- =
R 2 fo
1.224745R ---f-
fo
o g g
~F o
- = 7.55  10 6 m – 1 s – 1
fo = -------- f = ------


QS=ö=cbcilt
 NKS=cáåáíÉ=bäÉãÉåí=cçêãìä~íáçåë

Table 1.12 Corrected Manning roughness Table 1.13 Corrected Manning roughness
coefficient Mcorr for different applications in the coefficient M corr for different applications in the
case of Manning-Strickler law for 1D fracture case of Manning-Strickler law for 2D fracture
elements elements

Type M corr  M Type M corr  M

submerged rectangular cross- submerged rectangular cross-

A section 23 A section
2Bb -
 -------------------------
  b + B  A b - 2  3
 ------------
B B
 b + B
b
Note if b = B and b
2
A = B then
M corr = M , i.e.,
submerged slit plane
no correction is
B
needed
b 1
submerged slit plane B
B no correction is
2b  2  3
 ------- needed
b
 A
B open rectangular cross-section
C
open rectangular cross-section
C
B 2b  2  3
 ---------------
4Bb - 2  3 -
 ----------------------------  b + 2B
B   b + 2B  A b

b open wide channel (b > 20B)

D plane
open wide channel (b > 20B) 23
B 2 = 1.5874
D plane
4B  23
 ------- b
B  A
b

submerged circular cross-section

E
R 2R  2  3
 -------
 A

«aef=ö=ïïïKãáâÉéçïÉêÉÇÄóÇÜáKÅçã cbcilt=ö=QT