Fundamentals of Groundwater

Flow
 Darcy’s Law
• Groundwater in its natural state is invariably moving. This movement is
governed by established hydraulic principles.
• Darcy's Law along with the equation of conservation of mass is equivalent to
the groundwater flow equation, one of the basic relationships of
groundwater hydrology
• Darcy established it based on the results of experiments on the flow of
water through beds of sand.
• Darcy's Law states that the flow rate through porous media is proportional
to the head loss and inversely proportional to the length of the flow path.
Q h
qv  K  K i
A l
where constant of proportionality K = hydraulic conductivity (m/s); and i =h/l
= hydraulic gradient = rate of head loss per unit length of medium.
• The negative sign indicates that the total head is decreasing in the direction
of flow because of friction or resistance
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 2
Darcy’s Law
• Darcy's law is a simple mathematical statement which neatly
summarizes several familiar properties that groundwater
flowing in aquifers exhibits
– if there is no hydraulic gradient (difference in hydraulic head over a
distance), no flow occurs (this is hydrostatic conditions),
– if there is a hydraulic gradient, flow will occur from high head
towards low head (opposite the direction of increasing gradient,
hence the negative sign in Darcy's law),
– the greater the hydraulic gradient (through the same aquifer
material), the greater the discharge, and
– the discharge may be different through different aquifer materials
(or even through the same material, in a different direction) even if
the same hydraulic gradient exists.
• It is analogous to Fourier's law in the field of heat conduction,
Ohm's law in the field of electrical networks, or Fick's law in
diffusion theory.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 3
Hydraulic Head

• The energy available for groundwater flow is given the

name hydraulic head.

• It consists of three components, related to elevation,

pressure and velocity.
p v2
h z 
g 2 g

• In groundwater flow, the velocity is so low that the energy

contained in velocity can be neglected .
p
h z
g

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 4

Velocity potential

• Velocity potential is a scalar function of space and time  =

f(x,y,z,t) such that partial derivative with respect to any
direction gives the flow velocity in that direction .
  
vx  ; vy  ; vz 
x y z
• In groundwater the velocity potential may be defined as
 p 
   K   z   Const   Kh  Const
 g 
• Equipotentials are defined as curves of constant potential.
• The discharge vector is directed normal to the
equipotentials

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 5

Direction of Groundwater Flow

• The direction of groundwater movement is from higher

head to a lower head, or the groundwater flow is in the
direction of increasing velocity potential.

• The maximum hydraulic gradient is perpendicular to the

equipotential lines.

• At a minimum, it takes three hydraulic head measurements

to determine the hydraulic gradient and the direction of
groundwater flow.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 6

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 7
 Hydraulic Conductivity
• A medium has a unit hydraulic conductivity if it will transmit in
unit time a unit volume of groundwater at the prevailing
kinematic viscosity through a cross section of unit area
measured at right angles to the direction of flow, under a unit
hydraulic gradient.
• It has units of velocity.
• The hydraulic conductivity of a soil or rock depends on a variety
of physical factors, including porosity, particle size and
distribution, shape of particles, arrangement of particles, and
other factors.
• In general for unconsolidated porous media, K varies with
square of particle size; clayey materials exhibit low values of K,
whereas sands and gravels display high values.
• The hydraulic conductivity in saturated zones can be
determined by a variety of techniques,
including calculation from formulas, laboratory methods, tracer tests,
auger hole tests, and pumping tests of wells.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 8
• The hydraulic conductivity depends on both properties of the
porous medium and the fluid.

• For many ground-water studies, water is the fluid of interest,

providing more or less constant values of density and viscosity
(neglecting temperature dependencies).

• Thus, measurements of hydraulic conductivity are useful in

comparing differences in hydraulic behavior of the actual
materials.

• In more generally at systems where the fluids other than water

are present (such as, air, oil, and gasoline), hydraulic
conductivity becomes an awkward parameter because the
density and viscosity of the fluid vary together with the
medium properties.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 9

 Intrinsic Permeability
• The intrinsic permeability of a rock or soil is a measure of its
ability to transmit fluid as the fluid moves through it.

• Darcy's equation in a form of intrinsic permeability where the

properties of the medium and the fluid are represented
explicitly
kg h
v
 l
where, k is the intrinsic permeability and  is the
dynamic viscosity of fluid (Pa.s or kg/m/sec).

• The intrinsic permeability is independent of the fluid moving

through the medium and depends only upon the medium
properties.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 10
 Relation between Hydraulic Conductivity and Intrinsic
Permeability
kg 
K or k K
 
• The intrinsic permeability is used primarily when the density
or the viscosity of the fluid varies with position.
• The dimension of k is m2, but it is so small that square
micrometers (µ m)2 = 10-12 m2 is used.
• In the petroleum industry it is expressed in darcy (1 darcy =
0.98710-12 m2).
• For water, at 20 C and 1 atm pressure  = 998.2 kg/m3 and 
= 1.00210-3 Pa.s or kg/m/sec, so

K  9.77  10 6 k or k  1.023  10 7 K
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 11
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 12
Relation between Darcy velocity and Interstitial velocity
• Darcy velocity v is the apparent velocity or fictitious velocity or
Darcy flux (discharge per unit area).
• It is not the velocity which the water traveling through the
pores is experiencing.
• It assumes that flow occurs through the entire cross section of
the material without regard to solids and pores.
• Actually water can flow though pores only and the pore spaces
vary continuously with location within the medium.
• Therefore the actual velocity is nonuniform, involving endless
accelerations, deceleration, and changes in direction.
• To define the actual flow velocity or interstitial velocity, one
must consider the microstructure of the rock material.
• For naturally occurring geologic materials, the microstructure
cannot be specified three-dimensionally; hence, actual
velocities can only be quantified statistically.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 13
• Actually, the flow is limited to the pores only so that the
average interstitial velocity or actual velocity or seepage
velocity (vs) through pore space can be determined by
applying continuity equation.

Q  As v s  Av

A v
vs  v 
As 

• The pore velocity is the velocity a conservative tracer/dye

experiences if carried by water through the aquifer. This
indicates that for a sand with a porosity of 33 % the vs = 3
v.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 14
 Validity of Darcy’s Law
• Darcy law is valid for any Newtonian fluid.
• Though it was established under saturated flow conditions, it may
be adjusted to account for unsaturated and multiphase flows.
• Darcy’s velocity is proportional to the hydraulic gradient so it
applies to laminar flow in porous media.
• As inertial forces increase, turbulence occurs gradually. The
turbulence occurs first in the larger pore spaces; with increasing
velocity it spreads to the smaller pores. For fully developed
turbulence the head loss varies approximately with the second
power of the velocity rather then linearly.
• The Reynold’s number serves as a criterion to distinguish
between laminar and turbulent flow. Reynold's number for
porous media flow may be expressed as
vd10
NR 

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 15
• where d10 is the effective grain size.
• Darcy’s law is valid for NR< 1 and does not depart seriously up
to NR = 10.
• Fortunately, most natural underground flow occurs with NR < 1
so Darcy’s law is applicable.
• Deviations from Darcy’s law can occur where steep hydraulic
gradients exist, such as near pumped wells; also, turbulent
flow can be found in rocks such as basalt and limestone that
contain large underground openings.
• Darcy's law has been found to be invalid for very slow water
flow through dense clay; there the effects of electrically
charged clay particles on water in the minute pores produce
non linearity between flow rate and hydraulic gradient.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 16

Derivation of Darcy’s law from N-S equations
The various forces that may influence the motion of a fluid are
due to gravity, pressure, viscosity, turbulence, surface tension and
compressibility.
For viscous/laminar flow the force due to viscosity is
predominate over the forces due to turbulence, surface tension
and compressibility.
If we consider only gravity, pressure and viscous forces, the
equation of Newton’s 2nd law for motion of a fluid becomes
Ma  Fg  Fp  Fv
where M = mass, a = acceleration, Fg = Body Force, Fp = Pressure
Force, and Fv = Viscous Force.
This equation is known as Navier-Stokes’ equation which is useful
in the analysis of viscous/laminar flow.
N-S equation per unit mass is a  F / M  F / M  F / M
a p v

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 17

in x – direction this eq. becomes

v x v x v x v x 1 p  2
 vx  vy  vz    vx
t x y z  x 

For porous media v = vx/η etc., this equation for 2 –D (x – y

plane) case becomes
1 v x v x v x v y v x 1 p  2
 2  2    vx
 t  x  y  x 

Similarly in y – direction it is
1 v y v x v y v y v y 1 p  2
 2  2    vy  g
 t  x  y  y 

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 18

Analytical solution for this form of N-S equations in x and y
directions is impossible. We need approximations to make
these equations amenable for solution.
From dimensional consideration
 1 vx  vx  vy
 vx  2
2
 2
;  vy 
2

C L Cd Cd 2
where L = characteristic length = d (diameter of solid minerals
particles). From the definition of velocity potential
 
vx  ; vy 
x y
Therefore
1     1      1      1 p  1 
          
 t  x   x x  x   y y  x 
2 2
 x  Cd 2 x
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 19
Since  is a continuous function of x, y, and t; the order of its
partial derivatives may be changed, thus
2
1     1     1     1 p  1 
2

  2    2     
 x  t  2 x  x  2 x  y   x  Cd 2 x

Likewise in y-direction
2
1     1     1     1 p  1 
2

  2    2      g
 y  t  2 y  x  2 y  y   y  Cd y
2

Integrating wrt to y
2
1          p  
2
1 1
  2    2       gy  F ( x, t )
  t  2  x  2  y    Cd 2

Considering steady state and neglecting inertial terms

p  
  gy  Const
  Cd 2

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 20

Defining Cd2 = k = Intrinsic permeability and k / = K =
Hydraulic conductivity, gives
k  p  p 
   y   Const   K   y   Const   Kh  Const
    
Taking derivative results
h
v y  K ; v y   Ki y
y
which is Darcy’s equation in y- direction. Similarly we can obtain
corresponding equation in x- direction
 h
  K ; v x   Ki x
x x

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 21

This derivation resulted into the following important aspects
• Darcy’s law is not merely experimental result but it has theoretical
basis.
• N-S equations are for viscous/laminar flow, so Darcy’s equation is
applicable in laminar flow range.
• N-S equations are equations of motion (momentum equation), so
Darcy’s equation is a momentum equation instead of merely
empirical equation.
• Basis of expression for velocity potential in groundwater as  = -Kh
+ Const
• Basis of expression for Intrinsic permeability of the porous medium
as k = Cd2
• Relationship between Intrinsic permeability and hydraulic
conductivity as K = k /
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 22
 Homogeneity and Isotropy
• If hydraulic conductivity is independent of position within a
geologic formation then the formation is termed
homogeneous.
• If hydraulic conductivity is dependent of the position within a
geologic formation then the formation is termed
heterogeneous.
• If a formation has a hydraulic conductivity that is independent
of the direction of measurement it is called isotropic.
• If a formation has a hydraulic conductivity that is dependent on
the direction of measurement it is called anisotropic.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 23

Generalized Darcy’s equations
• In isotropic medium one of the coordinate axes can be aligned
along the direction of hydraulic gradient and the flow becomes
one-dimensional in that direction and the Darcy’s law can be
applied.
• But the hydraulic conductivity of material at a point in an
anisotropic medium exhibits directional dependency.
• In anisotropic medium or 3-dimensional flow, the generalized
Darcy’s law is expressed as

v   Kh

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 24

In Cartesian coordinates

v  vxi  v y j  vz k
   
h   i  j  k h
 x y z 
hydraulic conductivity tensor
 K xx K xy K xz 
 
K   K yx K yy K yz 
 K zx K zy K zz 

where i, j, and k are unit vectors along x, y, and z directions
respectively.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 25

 Darcy velocity components in x, y, and z directions are given
by
h h h
v x   K xx  K xy  K xz
x y z

h h h
v y   K yx  K yy  K yz
x y z

h h h
v z   K zx  K zy  K zz
x y z

K xy  K yx ; K xz  K zx ; K yz  K zy

If the principal directions of anisotropy coincide with the in x, y,

and z directions, the nondiagonal components of hydraulic
conductivity tensor are zero and thus
h h h
v x   K xx v y   K yy v z   K zz
x y z
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 26
 Laplace Equation
Physically all flow systems extend in three dimensions. However,
in many problems the features of the groundwater motion are
essentially planar or two dimensional and thereby reducing the
complexity of the analysis.
For an incompressible fluid if the flow is 2-D and steady then
equation of continuity is
 vx  v y
 0
x y
h h
Inserting Darcy (momentum) equation vx   K ; vy  K
x y
in continuity equation results
  h    h 
  K     K   0
x x  y y

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 27

 2h 2h
 0  2h  0
x 2
y 2

This is Laplace equation, thus Laplace equation is the governing

equation for 2-D steady flow of incompressible fluid in porous
medium.
By using velocity potential function in continuity equation

       
   0
 x   x   y   y 

hence
 2  2
 0   2  0
x 2
y 2

which is again Laplace equation in terms of velocity potential.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 28
Stream function
Stream function is a scalar function of space and time  = f
(x,y,t)such that its partial derivative with respect to any direction
gives the velocity component at right angles (in the clockwise
direction) to this direction
 
vx  ; vy  
y x
Equating the respective potential and stream functions of vx and vy,
we get Cauchy-Reimann equations
   
 ; 
x y y x
Substitution of stream functions of vx and vy, in continuity equation
for 2-D steady flow shows that it satisfies the equation of continuity
 2  2
 0
 x y  y x
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 29
Taking derivative of the first equation wrt y and of the second
equation wrt x of Cauchy-Reimann equations and subtracting
results to Laplace’s equation
 2  2
 0   2  0
x 2
y 2

Thus both  and  satisfy Laplace’s equation and hence they are
the conjugate harmonic functions.
The curves= f (x,y)= constant are the orthogonal trajectories of the
curves  = f (x,y)= constant.
It should be noted that within a given region of flow the streamlines
and equipotential lines are unique.
Hence in solving groundwater problems we need concern ourselves
only with the determination of one of the functions, subject to the
imposed boundary conditions.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 30
The other function will follow directly from the relationship
between them. For example
 
d  dx  dy
x y

So    
    dy  dx 

 x y 
Similarly    
    dx  dy 
 y x 

Laplace equation indicates that for conditions of steady-state,

laminar flow, the form of the groundwater motion can be
completely determined by solving one equation, subject to the
boundary conditions of the flow domain.
The Laplace equation for the stream function or velocity potential
or head is linear, so that these functions may be superposed.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 31
Flownet

• Solving groundwater flow problems amounts to solving the

differential equation of Laplace with the appropriate
boundary conditions.
• Only few solutions to Laplace’s equation in three dimensions
exist.
• Fortunately, many practical flow problems either are two-
dimensional or can be approximated by a two-dimensional
analysis.
• If analytical solution is difficult to obtain then several
problems can be solved graphically.
• Flownet is a graphical solution of Laplace equation consisting
of a set of flow lines orthogonal to a set of equipotential lines.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 32
 Streamline
If the path of an individual particle in flow field is tracked then the
tangent at any point gives the direction of velocity at that point and
the locus of all points defines the path of flow of an individual
particle known as flow line or streamline.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 33

 vy dy
 tan  
vx dx

v y dx  v x dy  0
Substituting vx and vy,in terms of stream functions
 
dx  dy  0  d  0
x y

 d    Const
Thus along a streamlime the stream function is constant.
Therefore the curves  = f (x, y) equal to a sequence of constants,
are at all points tangent to the velocity vectors and hence a set of
streamlines.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 34
Another important physical property of the stream function in
terms of Discharge between two stream lines can be obtained by
considering the flow between the two streamlines. Elementary
discharge between two stream lines through dy-section

dq  Vds  v y dx  v x dy 
 
dx  dy  d
x y

2
q   dq   d   2   1
1

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 35

• Thus total discharge between two streamlines (flow channel) is
a constant, so it behaves like a solid tube.

• Once the streamlines of flow have been obtained, their plot not
only shows the direction of flow but the relative magnitudes of
the velocity along the flow channels i.e., the velocity at any
point in the flow channel varies inversely with the streamline
spacing in the vicinity of that point.

• The discharge vector is normal to the equipotentials, the flow

occurs in the direction of the steepest slope of the plot of
potentiometric contours.

• Discharge streamlines are curves whose tangent at any point

coincides with the discharge vector.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 36
For 2-D steady flow, the equipotential lines is  = f(x, y) =
constant. By total derivative
 
d  dx  dy  0
x y
dy vx
v x dx  v y dy  0  
dx vy
For 2-D steady flow, along the stream line
dy v y
v y dx  v x dy  0  
dx v x
Product of these two is -1 so they intersect each other
orthogonally.
The streamlines ( = constant) and equipotential lines ( =
constant) form a grid of curves called a flow net wherein all
intersections are at right angles.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 37
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 38
 Flownet – important points
• Streamlines ( = constant) and equipotential lines ( = constant) form
a grid of curves called a flow net wherein all intersections are at right
angles.
• Along a streamline the stream function is constant.  functions exist
only for irrotational flow.
• Total discharge (quantity of flow) between two streamlines (flow
channel) is a constant, so it behaves like a solid tube.
• A flownet for an isotropic and homogeneous system provides a
graphical solution to the Laplace equation
• Flownets are usually drawn only for a 2D flow field.
• In heterogeneous media the flownet is not perfect squares.
• In anisotropic media, streamlines do not intersect equipotential lines
at right angles except when flow is aligned with one of the principal
directions of hydraulic conductivity.
• Several guidelines and procedures are available in the literature for
drawing flownets for variety of problems.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 39
Transformation of Medium

An equivalent transformed section is a fictitious section

such that the seepage discharge is the same as that for
an actual anisotropic and/or inhomogeneous section
under same head loss.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 40

Transformation of Inhomogeneous into Homogeneous
Medium
An equivalent transformed section is a fictitious section such
that the seepage discharge is the same as that for an actual
inhomogeneous section of same dimensions (thickness and
cross-sectional area) under same head loss.

Let there are n homogeneous and isotropic layers of thicknesses

d1, d2, d3, …, dn and hydraulic conductivities K1, K2, K3, …, Kn
respectively in a stratified (inhomogeneous) medium.

The equivalent transformed section has the hydraulic

conductivities Kx along the stratification and Ky normal to the
stratification and dimensions Le = L and de = d1 + d2 + d3 + … + dn

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 41

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 42
 Flow Parallel to Stratification
Consider flow parallel to the stratification (x – direction).
Seepage discharges per unit width of the medium in each layer
are
q1  d1 K1i1 q 2  d 2 K 2 i2 q3  d 3 K 3i3 q n  d n K n in

The conditions for this case are

h h2  h1
i1    i2  i3  ...  in  i x
l L
q x  q1  q2  q3  ... qn
For the equivalent transformed section, let discharge, hydraulic
conductivity, and hydraulic gradient in the x-direction are qx, Kx,
and ix respectively, so that

q x  d e K x i x  (d1  d 2  d 3  ...  d n ) K x i x

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 43

Therefore
q x  (d1  d 2  d 3  ...  d n ) K x i x  d1 K1i1  d 2 K 2 i2  d 3 K 3i3  ...  d n K n in

Hence
d1K1  d 2 K 2  d 3 K 3  ...  d n K n
Kx 
d1  d 2  d 3  ...  d n

n n
K x   K jd j d j
j 1 j 1

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 44

Flow Normal to Stratification
Consider flow normal to the stratification (y – direction).
Seepage discharges per unit width of the medium in each layer
are
q1  LK1i1 q2  LK 2i2 q3  LK3i3 qn  LK nin

where hydraulic gradient in each layer is given by

h hL1 h hL 3 hLn
i1   ; i2  L 2 ; i3  ; in 
l d1 d2 d3 dn

hL1  d1i1 ; hL 2  d 2 i2 ; hL3  d 3i3 ; hLn  d n in

For the equivalent transformed section, let discharge, hydraulic

conductivity, and hydraulic gradient in the y-direction are qy, Ky,
and iy respectively, so that
q y  Le K y i y  LK y i y
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 45
 h hL hL
iy   
l d e d1  d 2  d 3  ...  d n
here the total head loss hL is the sum of head losses in each
layer i.e.
hL  hL1  hL 2  hL3  ... hLn
Therefore
(d1  d 2  d 3  ...  d n )i y  d1i1  d 2 i2  d 3i3  ...  d n in
Seepage equations give
qy q1 q2 q3 qn
iy  ; i1  ; i2  ; i3  ; in 
LK y LK 1 LK 2 LK 3 LK n

from the condition of continuity of flow in the y-direction

q y  q1  q2  q3  ... qn

and hence
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 46
 d1  d 2  d 3  ...  d n d1 d 2 d 3 dn
    ... 
Ky K1 K 2 K 3 Kn

n n
Ky  d j  (d j /Kj)
j 1 j 1

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 47

Comparison of hydraulic conductivities along and normal the
stratification
The hydraulic conductivity along the stratification is always
greater than the hydraulic conductivity normal to the
stratification i.e. Kx > Ky or
n n n n

K
j 1
j dj d
j 1
j  d j
j 1
 (d
j 1
j /Kj)

To prove it consider a two layer porous medium

K 1 d1  K 2 d 2 d  d2
 1
d1  d 2 d1 d 2

K1 K 2

 d1 d 2 
K1d1  K 2 d 2     d1  d 2 2
 K1 K 2 

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 48

 K 2 K1
 2  K12  K 22  2 K1 K 2
K1 K 2

and hence

K 1  K 2  2
0

which is always true so Kx > Ky.

This can also be proved for three layered or more layered soil
medium.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 49

Transformation of Anisotropic into Isotropic Medium
In Figure, let ‘s’ and ‘n’ represent the directions of the tangent
to the flow line and the normal to the equipotential line,
respectively.
In an isotropic flow medium the flow lines and equipotential
lines form an orthogonal system, hence the ‘s’ and ‘n’
directions are identical (β = 0).
In anisotropic flow, the direction of the stream lines in
general, will not coincide with the direction of the normal to
the equipotential lines and hence β ≠ 0.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 50

From Darcy’s law, the flow velocity along the steam line (‘s’
direction) is given by
dh
vs   K
ds
where K = coefficient of permeability in the s-direction; vs =
velocity in s-direction; dh/ds = hydraulic gradient in s-direction.
Hence, the velocity components in the x and y directions are
h h
vx  K x v y  K y
x y
The velocity components along x and y directions are

v x  v s cos  v y  vs sin 

therefore

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 51

 h h  v s
vx  K x  v s cos    cos 
x x K x
h h  vs
vy  K y  vs sin    sin 
y y K y
dx dy
Also from Figure  cos   sin 
ds ds
dh h dx h dy
Since  
ds x ds y ds
By substituting the corresponding values
dh  v s  v s  vs
  cos 2   sin 2 
ds K  Kx Ky

1 cos 2  sin 2 
 
K Kx Ky
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 52
From rectangular coordinates, we have

x  r cos  y  r sin 
By substituting the corresponding values

r2 x2 y2
 
K K x K y
This is equation of an ellipse, which has Major semi-axis = K x
and Minor semi-axis = K y

In case of a stratified medium with the x-direction parallel to the

bedding plane, the hydraulic conductivity in any direction,
making an angle ‘α’ with the x-axis, can be obtained from the
graphical construction as in Figure, is called the ellipse of
direction.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 53
 Direction of Ellipse

Therefore, the square root of the directional hydraulic

conductivity for a homogeneous and anisotropic layer when
plotted from a point will generate an ellipse.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 54
X-Transformation
Using velocity components from the velocity potential
functions in the continuity equation
  h    h 
  K x     K y   0
x x  y y
If the medium were homogeneous then
2 h 2 h Kx 2 h 2 h
Kx  Ky 0   0
x 2
y 2
Ky  x 2
y 2

2 h 2 h 2 h 2 h
 0   0

 x Ky Kx 
2
y 2
X 2
y 2

where
X  x Ky Kx

Substituting this in direction of ellipse equation, we get

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 55

 r2 X 2 y2
 
K K y K y

Ky
X y 
2 2
r2  R2
K

R  r K y K

If r  K  then R  K and expression is an equation of circle

y

which demonstrate that, by the transformation of scale X  x K K y x

in the x-direction, the ellipse of direction will be transformed in

to a circle where the hydraulic conductivity will be an invariant
with the direction.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 56
Y-Transformation
Using velocity components from the velocity potential
functions in the continuity equation
  h    h 
  K x     K y   0
x x  y y

If the medium were homogeneous then

2 h 2 h 2 h K y 2 h
Kx  Ky 0  0
x 2
y 2
x 2
Kx  y 2

2 h 2 h 2 h 2 h
 0   0
x 2
 y Kx Ky  
2
x 2
Y 2

where Y  y Kx Ky

Substituting this in direction of ellipse equation, we get

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 57
 r2 x2 Y 2
 
K K x K x

Kx 2
x Y 
2
r  R2
2

K

R  r K x K

If r  K  then R  K and expression is an equation of circle

x

which demonstrate that, by the transformation of scale Y  y K K x y

in the x-direction, the ellipse of direction can be transformed in

to a circle wherein the hydraulic conductivity is direction
independent
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 58
Equivalent Hydraulic Conductivity
If the medium is layered parallel to x- axis then Kx > Ky so X  x K K
y x

transformation results into a contraction along X axis to have the

circle while Y  y K x K y transformation results into an expansion
along Y axis to have the circle.
Thus by simple expansion or contraction of spatial coordinates a
given homogeneous and anisotropic flow region can be
transformed into a fictitious isotropic region where in Laplace’s
equation is valid and consequently the potential theory is
applicable. This fictitious flow region is called the transformed
section.
The equivalent hydraulic conductivity for the transformed section
can be computed by equating seepage discharges in anisotropic
(original) and isotropic (transformed/equivalent) mediums in x-
direction as
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 59
 dh dh
q  AB K x  A' B' K e
dx dX

in X-transformed where X  x Ky Kx but AB = A’B’ as there is no

change in y-dimensions thus
dh dh Kx
Kx  Ke  Ke  K x K y
dx dx Ky

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 60

Similarly for Y-transformed case the seepage discharges in x-
direction dh dh
q  AB K x  A' B ' K e
dx dx
where A' B'  AB K x K y but dh/dx is the same as there is no change
in x-dimensions thus
Kx
AB K x  K e AB  Ke  K x K y
Ky

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 61

Therefore from both transformations the equivalent hydraulic
conductivity is same

Ke  K x K y
Once the problem has been solved for the transformed section,
the solution for the natural medium can be obtained by applying
the inverse of the scaling ratio.

Thus, the effects of anisotropy can be taken into account by a

simple transformation of spatial coordinates.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 62

Different Boundaries in a Porous Medium

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 63

Boundary 1: Interface between soil and water or boundaries of
the reservoirs
Along the boundaries of the reservoir the pressure distribution
may be taken as hydrostatic. Therefore the pressure at any point
(x, y) on this boundary in Figure is
p( x, y )  patm   wh

but h  h1  y
p( x , y ) patm p( x , y ) patm
   (h1  y )  y  h1
w w w w

p( x , y ) patm
  K (  y)  C   K (  h1 )  C1 = constant
w w

So  is constant at Boundary 1 that means Boundary 1 is an

equipotential line
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 64
Boundary 2: Impervious boundary
At impervious boundaries the fluid can neither penetrate the
boundary nor leave gaps; thus the velocity component normal to
the boundary at any point must vanish thus

Vn  0 

n
0    constant
Thus  is constant along ‘n’ direction i.e. normal to the impervious
boundary

Vn  0   0   constant
t
Hence  is constant along ‘t’ direction (along the impervious
boundary) that means Boundary 2 is a flow/stream line. Likewise,
any streamline satisfies the condition for an impervious boundary
and may be taken as such.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 65

Boundary 3: Phreatic line or line of seepage
• The line of seepage is the upper most streamline in the flow
domain.
• It separates the saturated region of flow from that part of the
porous medium through which no flow occurs.
• The determination of its locus is one of the major objectives of
groundwater investigations.
• In addition to the requirement that the this boundary is a
streamline  = constant.
• The pressure at every point along its surface is constant and
equal to atmospheric pressure.
At any point on the streamline
p( x , y ) patm
  K (  y)  C     K (  y)  C
w w
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 66
 K atm
Therefore    Ky   C = constant    Ky = constant
w
which demonstrates that the velocity potential (and the total
head) along the line of seepage varies linearly with elevation
head. Hence
1  Ky1  2  Ky2 (2  1 )
so  y1  y2
K
similarly (3  2 )
 y2  y3
K

For successive equipotential lines

(2  1 )  (3  2 )  y1  y2  y2  y3  y1  y2  .......... =const
This requires constant vertical intercepts at the points of
intersection of the line of seepage with successive equipotential
lines of equal drops ( = constant).
Hence Boundary 3 is the top flow line ( is constant) and  varies
linearly.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 67
Boundary 4: Interface between soil and air (Seepage face)

It is a boundary where the seepage leaving the flow region enters

a zone free of both liquid and soil.

As the pressure on this surface is constant and atmospheric similar

to line of seepage thus

  Ky = constant

Therefore the velocity potential varies linearly but  = constant.

So Boundary 4 is not a flow line.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 68

Boundary 5: Flow Through Two Media (Flow crossing
anisotropy)

Consider a case where flow

passes from a region of
hydraulic conductivity K1 to
other region of K2. Due to
difference in conductivities
the direction of flow changes
from 1 in region K1 to 2 in
region K2
The velocity potentials in
regions K1 and K2 are given
by
 p1   p2 
1   K1 h1  C1   K1   y1   C1
  2   K 2 h2  C 2   K 2   y 2   C 2
w  w 
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 69
At the common point O at the interface, p1 = p2; y1 = y2 and C1 = C2
which may be assumed zero. With these conditions
1 K1

2 K 2
From continuity condition, the normal components of flow
approaching and leaving the boundary/interface must be equal;
hence the normal velocity components must be such that
v1n  v2 n or v1 cos1  v2 cos 2 hence dh1 dh2
K1 cos  1  K 2 cos  2
ds1 ds 2

But between the same two equipotential lines dh1 = dh2 and also
ds1  b sin  1 and ds 2  b sin  2 thus

cos  1 cos  2 K1 tan  1

K1  K2  
sin  1 sin  2 K 2 tan  2
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 70
Therefore the flowlines deviate towards the normal when they
cross from a coarse (higher K) to fine (lower K) medium and
deviate away from the normal when they cross from a fine (lower
K) to coarse (higher K) medium and vice-versa.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 71

 Theory of Groundwater Flow

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 1

Governing Equations in Cartesian Coordinates

• Groundwater in its natural state is invariably moving.

• This movement is governed by established hydraulic
principles
• The differential equations for groundwater flow can be
developed from principles of mass conservation and
momentum (Darcy’s) equation.
• Conservation of mass in a representative elementary
volume(REV) can be used to derive governing equations

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 2

A REV is defined as a volume exhibiting the average
properties of the porous media around a point P(x, y, z),
which is the center of the volume. Let, w is the density of
water and η is the porosity in the REV.

For reference purposes,

let's label the six faces of
the REV as x1 (ABCD), x2
(HGFE), y1 (EFBA), y2
(DCGH), z1 (HEAD), and
z2 (GFBC). Now, assume
that the water mass
fluxes through the six
faces are Mx1, Mx2, My1,
My2, Mz1, and Mz2,
respectively.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 3
 Continuity Equation

The mass conservation statement can be stated as

mass inflow rate - mass outflow rate = change of mass storage

with time

The initial mass within the REV is

M   wxyz
Let, also there exists a groundwater source of strength R
(volumetric source rate per unit volume), so the mass inflow rate
(Ms) due to source is

M s   w R( x, y, z, t )xyz
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 4
The mass conservation equation may be rewritten in mathematical
terms as
M
M x1  M y1  M z1  M s   M x 2  M y 2  M z 2  
t
The net inflow rate into the REV in the x-direction is the difference
between the inflow and the outflow rates would be
 M x1  M x1   w v x 
M x1  M x 2  M x1   M x1  x    x   xyz
 x  x x

Similarly, the net flow rates into the REV in the y- and z-directions
are
M y1  w v y 
M y1  M y 2   y   xyz
y y
M z1   w v z 
M z1  M z 2   z   xyz
z z
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 5
Substituting these in mass conservation equation and dividing both
sides by volume of REV (i.e. xyz)

   wv x    wv y    wvz   1   wxyz 
       w R 
 x y z   x  y z t

This is continuity equation for flow through porous media.

This can be manipulated for different cases.
Assuming, the density of the fluid does not vary spatially and
aquifer is incompressible, the density term vanishes from both
sides so
 v x   v y   v z  1  xyz  h
  R  Ss
x y z x y z t t
where Ss is the specific storage.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 6

Use of Darcy’s law (momentum equation) for velocities vx, vy and vz
in continuity Equation results

  h   h    h h
 K x    K y    K z   R  S s
x x y y z z t

This is the main equation of groundwater flow in saturated media.

It can be written in many forms that apply to a variety of different

conditions.

For example, (1) if there is no source or sink in the control volume

then R(x, y, z, t) = 0 and Eq becomes
  h   h   h h
 K x    K y    K z   S s
x x  y y z z t
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 7
Steady State flow
with source
  h   h   h
 K x    K y    K z   R( x, y, z )  0
x x  y y z z
without source
  h   h   h
 K x    K y    K z   0
x x  y y z z
If the porous medium is homogeneous and isotropic (Kx = Ky = Kz),
it simplifies to the Laplace equation

2h  2h 2h
  2 0
x 2
y 2
z

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 8

Unsteady State flow in homogeneous and isotropic porous
medium
with source

 2h  2h  2h R Ss  h
  2  
x 2
y 2
z K K t

without source

 2h  2h  2h Ss  h
  2 
x 2
y 2
z K t

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 9

Diffusion Form of Groundwater flow equation
Assuming constant specific storage and dividing both sides by Ss,
the equation is transformed into

  Kx  h    Ky  h    Kz  h  R  h
         
 x  Ss  x   y  Ss  y   z  Ss  z  Ss  t

  h   h    h R h
 Dx    D y    Dz   
 x   x   y   y   z   z  Ss  t

where Dx=Kx/Ss, Dy=Ky/Ss, and Dz=Kz/Ss are called hydraulic

diffusivities.

Diffusivities and hydraulic gradients in the x-, y-, and z-directions

are the determining factors.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 10

Multiplying both sides of by the aquifer thickness (b) gives

  h   h   h h

 Tx    T y    Tz   S
x x  y  y z z t

where Kxb = Tx; Kyb = Ty; Kzb = Tz are transmissivities and

storativity S = Ssb.

This form of the groundwater flow equation is solved in

MODFLOW.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 11

Compressibility of Aquifers
To treat the compressibility of aquifers the following
assumptions are made:
• The elemental volume is constrained in lateral directions and
undergoes change in length only in z-direction.
• The pore water is compressible.
• The solid grains of the aquifer are incompressible but the
pore structure is compressible.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 12

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 13
Noting that the mass in a REV for a compressible aquifer

M   wxyz
The change in M with time happens as, η and z may vary due to
vertical compression or expansion of the medium, and w may
change in time as well as in space.

But variations of the lateral dimensions of aquifers subject to

variable flow are negligible because of the constraints of these
aquifers by their surroundings.

So
M   wz    z    w 
 xy    w   w z  z  xy
t t  t t t 

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 14

Compressibility of Aquifer Material
Concept of vertical compressibility of the granular skeleton of
the medium, treated as continuum, is used.
The stress z on the intergranular skeleton in the vertical
direction is called intergranular pressure or stress.
Compressibility α = 1/Eb, where Eb is the bulk modulus of
elasticity of this skeleton defined by
change in pressure d z 1
Eb    
change in volume per unit volume d (z ) z 
So that

(z )  z
d (z )   z d z    z
t t

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 15

Compressibility of Voids
Volume of voids = ηxyz so the volume of solid grains Vs = (1-
η)xyz which may be considered as constant because the
compressibility of the individual grains is considerably smaller than
that of their skeleton and is also smaller than the compressibility of
water.
The total derivative of this quantity is zero, or
dVs  d (1   )xyz   0
Again x and y (lateral dimensions of the volume element) are
constant as they do not change in comparison to the change in the
vertical z and thus
 1    ( z )
z d (1   )  (1   )d ( z )  0  
t z  t
  z
 (1   )
t t
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 16
Compressibility of Fluid
Compressibility β of the fluid is the reciprocal of its bulk modulus
of elasticity dV f V f d V f 
      dp
dp Vf
Here p is the pressure experience by water in the pores and is
called pore pressure or neutral pressure. At any depth
intergranular pressure and pore pressure add to combined
pressure so  z  p  const  d z  dp
From conservation of mass  wV f   w0V f 0  const
Total differentiation of this equation gives
d V f  d  w 
 w d V f   V f d  w   0  
Vf w
Which can be manipulated as
 w p  z
 w   w 
t t t
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 17
Overall Compressibility

M  z p
  w   w (1   )   w xyz   w xyz    
t t t

h
p
 z  const p   w g h  z  const  
w g
Taking partial derivative and manipulating
p h p  w
 w g 
t t  w t
p 1 h  w 1 h
 wg or   w2 g
t 1   p t t 1  p t
Substituting
M h
 xyz    
1
 w2 g
t 1  p t
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 18
Comparing equations

 w g       w g    
1
Ss 
1  p
For practical values of p in groundwater flow and for bulk modulus β
for water, the ratio 1/(1- β p) may be replaced by one as a first
approximation.

Two parts of specific storage may be interpreted as

•  w g = water in storage released due to the compression of the
intergranular skeleton per unit volume and per unit decline of
head
•  w g = water in storage released due to the expansion of the
water per unit volume and per unit decline of head.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 19

It may be observed that β > 0 for a decrease in pore pressure and a
resulting expansion of the water.
In the computation of Ss, the compressibility of water cannot be
neglected because the term  w g may be of the same order of
magnitude as the term  w g
The compressibility of aquifer material is relatively important only
when the material is completely saturated with water and confined
between impervious strata.
When the flow is unconfined the compressibility of the material and
the water are relatively unimportant compared to unsteady
perturbations or vertical displacement of the free surface which
affect the flow pattern.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 20

Governing Equation Considering Compressibility
If flow is also compressible i.e. the density of the fluid varies
spatially, then the continuity equation can be expanded as
   w    w    w     v x   v y   v z   1 M
  v x  vy  vz    w       w R 
 x y z   x y z  xyz  t
We know that
h
p
 z  const or p   w g h  z   const
w g
1
Taking partial derivative, using dp  d w and manipulating
w 
p h p  w  w 1 h
 w g     w2 g
x x  w x x 1   p x
p h p  w  w 1 h
 w g     w2 g
y y  w y y 1   p y

p  h  p  w  w 1 2  h 
  w g   1      w 
g  1 
z  x   w z z 1   p  x 
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 21
Use of Darcy’s law for velocities vx, vy and vz and substituting other
terms in LHS and manipulating
  h
2
 h 
2
 h   h  
1 2 
 w g K x    K y    K z   1   w R 
1  p   x  y  z   z 
 
   h    h     h  1 M

 w   K x    K y    K z 
  
x x y  y   z   z   xyz  t

substituting M/t in RHS


 w g    h 
2
 h 
2
 h   h  
K x    K y    K z   1  R 
1  p    x  y  z   z 
 
   h    h     h   w g h
  K x
x    K    K 
  y  y  y   z  z  z      
   x      1   p t

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 22

which is the most general governing equation as it considers
compressibility of flow and aquifer, heterogeneity, anisotropy, and
unsteady condition.
Different particular cases can be deduced from it.
Since 1 – βp ≈ 1

   h 2   
2
    
 w g  K x    K y    K z  1   R 
h h h

 x  y  z   z 
 
   h    h     h  h
  x  K x  x    y  K y  y    z  K z  z     w g     t
 
      

For incompressible flow β = 0 and this equation reduces to earlier

derived equation.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 23

Solution and Boundary Conditions
A governing partial differential equation has an infinite number of
possible solutions, each of which corresponds to a particular case of
flow through a porous medium domain.
To obtain one particular solution (from infinite solutions)
corresponding to a specific problem of interest, it is necessary to
provide supplementary information e.g.

– Geometry of flow domain

– Values of physical coefficients
– Initial conditions
– Boundary conditions

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 24

The partial differential equation together with the supplementary
information defines an individual problem leading to a particular
solution.
Different boundary and initial conditions lead to different solutions
and hence these should be determined or assumed correctly on the
basis of available information and past experience.
Initial conditions include the specification of h at all points within the
domain D at some initial time usually denoted as t = 0.
This can be written schematically as specifying h  f x, y, z, 0

for all points x, y, z inside D; f is a known function.

The various types of boundary conditions in flow through porous

media are described in subsequent slides.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 25
1. Boundary of prescribed head/potential or Dirichlet boundary
or First type boundary

The head h is prescribed for all points on this boundary S

h  f1 x, y, z, t  on S
where f1 is a known function. A special case is
h  h0 or h  Const on S

wherein head is constant on the boundary or the boundary is an

equipotential surface.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 26

2. Boundary of prescribed flux or Neumann boundary or Second
type boundary

On this type of boundary, the flux (derivative of head/potential)

normal to the boundary surface is prescribed for all points i.e.
h
qn  C  f 2 x, y, z, t  on S
n
where f2 is a known function.

A special case is an impervious or no flow boundary where the flux

normal to the boundary vanishes everywhere.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 27

3. Semipervious Boundary or Cauchy boundary or Third type
boundary or Mixed boundary

On this type of boundary, a functional relationship between the flux

and head/potential is prescribed i.e.

qn  hx, y, z, t  on S

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 28

4. Unsteady Phreatic Surface with Accretion

The location and shape of the free surface are unknown, in fact,
their determination constitutes part of the required solution.
Since the pressure at all points of the phreatic surface S is
atmospheric thus

hx, y, z, t   Kz  Const on S
which in fact gives at any time t a relationship between the
coordinates of points of the phreatic surface.
Therefore it may be considered equivalent to F x, y, z, t   Const
describing the geometry of this surface
F x, y, z, t   hx, y, z, t   Kz  Const on S
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 29
The unsteady phreatic surface with accretion is a surface on which a
certain property is maintained constant. For such a surface
2 2
 h  h  h h   h h

Kx 
 
 Ky 
  Kz 
  
 1  R1  
 
x y zz   z t
The difficulty with this boundary condition stems from the fact that
the distribution h  x, y , z, t  and hence F  x, y , z, t  is unknown before
the problem is solved.
In fact, in order to determine h we must know the boundary’s
location F and in order to know where this boundary is, we have to
know h.
Iterative techniques are used to overcome this difficulty in numerical
techniques on the other hand analytical techniques are seldom
capable to solve a problem with this complicated boundary
condition except for some particular classes of problems.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 30
5. Seepage Face
Along a seepage face water emerges from the flow domain trickling
downward to the adjacent body of water.
The pressure along a seepage face is atmospheric pressure and
hence similar to phreatic surface the boundary condition is

hx, y, z, t   Kz  Const on S
However the geometry of boundary S along the seepage face is
known ie F x, y, z, t   Const except for its upper limit which is exit
point of unknown phreatic surface.
The location of this point is part of the required solution.
In unsteady flow the location of the upper limit of the seepage face
varies with time.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 31

6. Surface of Discontinuity in Conductivity

On any common point on the surface of discontinuity the elevation

and pressure are same and also the flux across the boundary should
be same from continuity requirement i.e.

h1 h2
h1 x, y, z, t   h2 x, y, z, t  and qn1  qn 2 ; K1  K2
n n

where n is distance measured along the normal to the boundary at

the common point.
Thus two boundary conditions to be satisfied which requires
simultaneous solution of two equations for h1 and h2.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 32

7. There may be additional types of boundary conditions, for
example if effect of capillary fringe is considered.
Finally any problem of groundwater flow can be stated
mathematically by

• Defining the flow domain; part of the boundary may be at

infinity; sometimes the domain boundary is unknown a priori
• Specifying the governing partial differential equation to be
satisfied by h (other variables p, , , or  can be used in place of
h) at all points within the flow domain and t > 0
• Specifying the boundary conditions that must be satisfied by h at
all points of the boundary and at all times
• Prescribing initial conditions at all the points of the flow domain,
when the flow is unsteady.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 33

• Once a problem is mathematically formulated, it can be solved by
an appropriate technique.
• A well posed mathematical problem satisfies existence,
uniqueness and stability (three basic requirements) of the
solution i.e. the solution must exist, the solution must be
uniquely determined and the solution should depend
continuously on the data (be stable).
• In majority of cases the governing partial differential equation,
boundary conditions, domain geometry or material properties or
their combination is not simple to be solved by analytical
techniques and hence numerical techniques are adopted.
• Analytical solutions are possible for simple cases and few of such
cases are dealt in the following section.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 34

Steady Unidirectional Flow in Confined Aquifer
Let groundwater flow with a velocity vx in the x-direction of a
confined aquifer of uniform thickness.
Then for steady flow Laplace equation reduces into the following
ordinary differential equation
d 2h
2
0
dx
h  C1 x  C 2
h0  h1
h  h0  x
K
vx
 h0  x
K
head varies linearly
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 35
Steady Flow in a Varying Thickness Confined Aquifer
Governing equation for a steady groundwater flow in a
homogeneous isotropic (wrt K) confined aquifer of varying
thickness in the x-direction becomes

  h  d  dh 
 Tx 0   bx 0
x  x  dx  dx 
where the thickness of aquifer at section x is
b0  b1
bx  b0  x
L

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 36

Integrating this ordinary differential equation
dh C1 C1 L  b b 
 h ln b0  0 1 x   C 2
dx bx b0  b1  L 
Boundary conditions h = h0 at x = 0 and h = h1 at x = L can be used
to find C1 and C2 so

h0  h1 b0  b1 C 2  h0  h0  h1 
ln b0
C1   ln(b0 / b1 )
L ln(b0 / b1 )

h0  h1
h  h0  ln b0  lnb0  (b0  b1 ) x / L 
ln(b0 / b1 )

h0  h ln(b0 / bx ) ln b0  ln bx
 
h0  h1 ln(b0 / b1 ) ln b0  ln b1

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 37

The discharge per unit width at any section x can be found by

dh h0  h1 (b0  b1 )
q  Ax .v x  bx .1.() K   KC1  K
dx L (ln b0  ln b1 )
This shows that q is independent of x as it should be from the
continuity of flow.
Time to travel from upstream to downstream reservoir

dx vx dx
vs    dt  
 dt vx
 L2 b0  b1 ln b0  ln b1
L
t
KC1  b dx  2 K h
0
x
0  h1 b0  b1

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 38

 Dupuit – Forchhemeir Assumptions and 2-D Equations for
Unconfined Flow
If there is no change of hydraulic head in the vertical direction, the
relevant equations can be reduced to a two-dimensional
groundwater flow equations by dropping z – directions terms.

Consider the prism ABCDEF in

which the bottom is parallel to the
direction of flow and the top
surface is the water table.
The datum plane for h coincides
with the upper surface of the
impervious stratum underlying
the aquifer.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 39

Applying mass balance for the REV similar to confined aquifer
except aquifer height h (equal to saturated thickness of aquifer =
height of water table from the bed rock)
M x M y M
M x1  M x 2  M y1  M y 2   w Rxy   x  y   w Rxy 
x y t
here recharge rate R is additional inflow per unit area into REV.
In an unconfined aquifer, direct analytical treatment is not
possible.
The difficulty arises from the fact that the water table in the 2-D
case represents a flow line.
The shape of the water table determines the flow distribution, but
at the same time the flow distribution governs the water table
shape.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 40

To obtain a solution, Dupuit assumed

(1) the velocity of the flow to be proportional to the tangent of the

hydraulic gradient instead of the sine, and

(2) the flow to be horizontal and uniform everywhere in a vertical

section.

These Dupuit – Forchhemeir assumptions, although permitting

solution to be obtained, limit the application of the results.

With these assumptions

h h
M x   w v x hy; M y   w v y hx; vx  K x ; v y  K y
x y

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 41

For incompressible flow
  h   h h
 K x h    K y h   R  S
x x  y y t

   h2     h2  h
 K x    K y   2 R  2 S
x x  y y  t

For water table aquifer storage coefficient= specific yield (S = Sy).

If medium is isotropic then Kx = Ky = K and we have

 2 h 2  2 h 2 2 R 2S y  h
  
x 2
y 2
K K t
This relation is known as Boussinesq Equation.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 42

Boussinesq equation can only be used when vertical hydraulic
gradient is small and negligible.

For general problems of flow unconfined aquifer, general equation

should be used.

If there is no recharge and flow is steady then

 h 2
 h2 2 2
  h 0
2 2

x 2
y 2

which is the Laplace equation in h2 for unconfined flow.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 43

Steady Unidirectional Flow in Unconfined Aquifer
For steady unidirectional flow in an unconfined aquifer without
recharge, the governing equation reduces to
d 2h2
2
0
dx
h 2  C1 x  C 2
dh
2h  C1
dx

dh
q   Kh 
dx
dh q
h 
dx K

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 44

 C1  2q / K C 2  h02
Thus the equation of water table is 2q
h   x  h02
2

K
which indicates that the water table is parabolic in form.
Neglecting the existence of a seepage face at the exit i.e. the
condition that for x = L, h = h1 may be used to determine q as
q
K 2
2L

h0  h12 
Slope at entry
dh q

dx Kh 0
The travel time between two points in the unconfined aquifer is
dx 2
x2
 x2
2q
KC1 x1 
dt  t  hdx  h02  x dx
vx q x1
K
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 45
In the direction of flow, the parabolic water table increases in slope and
two Dupuit assumptions become increasingly poor approximations to the
actual flow; therefore the actual water table deviates more and more from
the computed position in the direction of flow.
The fact that the actual water table lies above the computed one is that the
Dupuit flows are all assumed horizontal, whereas the actual velocities of
the same magnitude have a downward vertical component so that a
greater saturated thickness is for the same discharge.
At the downstream boundary a discontinuity in flow forms because no
consistent flow pattern can connect a water table directly to a downstream
free water surface.
The water table actually approaches the boundary tangentially above the
water body surface and forms a seepage face.
However for flat slopes, where the sine and tangent are nearly equal,
Dupuit’s solution closely predicts the water table position except near the
outflow.
The solution also, accurately determines q or K for given boundary heads.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 46
Steady Unidirectional Flow with Recharge in Unconfined Aquifer
For steady unidirectional flow with Recharge in an unconfined
aquifer without recharge, the governing equation reduces to

d 2 h 2 2R
2
 0
dx K

dh 2 R
2h  x  C1
dx K

R 2
h2   x  C1 x  C 2
K

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 47

Boundary conditions h = h0 at x = 0 and h = h1 at x = L can be used
to find C1 and C2 so C2  h02 and 2 2
h0  h1 RL
C1   
L K
Therefore
R h 2
 h 2
RL
h  h0  x 
2 2 2 0 1
x x
K L K
which shows that water table varies elliptically.
The water table first rises and after attaining hmax it falls.
At hmax the water table slope is zero and known as water divide
line.
At the water divide line seepage discharge is zero and on one side
of it the seepage is towards one stream and on the other side the
seepage is towards the other stream.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 48
Thus the location of the water divide line d can be computed using
dh/dx = 0 condition as
d
KC1 L
2R
 
K
2 2 RL

h02  h12 
h02  h12
 h  d L  d  
2 R
2
h max 0 d
K L
The discharge per unit width at any section x can be found by
dh h02  h12 RL
q  Ax .v x  h.1.() K K   Rx
dx 2L 2
The above equation can be used at x = 0 and x = L to find the base
flows in the streams.
It should be noted that total base flow into both streams is equal
to the total recharge RL between the streams.

Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 49

The travel time between two points for this case is
dx L2
h
dt   t  2  dx
vx L1
2 Rx  KC1
This is a complicated function difficult to integrate analytically, though
it may be reduced into elliptical integral form and then the solution can
be adopted from the standard formulae available in the literature.
For approximate computation of time the variation in the head from
water divide line to the stream may be assumed linear, which is
equivalent to assumption of constant hydraulic gradient.
For example, the travel time from water divide line to the left stream
comes out dx  d d  d2
K h
dt  t  dx 
dh  h0 K hmax  h0
K 0 max
dx
Similarly travel time on the other side of the water divide line can be
computed.
Prof B R Chahar, IIT Delhi Groundwater Hydrology - McGraw Hill 50