Adam Říčka
Department of Geological Sciences, Faculty of Science, Masaryk University,
Kotlarska 2, 611 37 Brno, Czech Republic
Groundwater flow and transport modeling
Introduction to ground-water modeling
Modeling case studies
t
h
SW
z
h
K
zy
h
K
yx
h
K
x
szzyyxx
Content of the workshop
• Introduction to basic principles of
gorundwater flow modeling
• Overview of important input data
• Set-up of model geometry nad boundary
conditions
• Model calibration and verification
• Practical using of Processing Modflow –
PMWIN
Aim of the lectures + workshop
• Model is any device representing a field situation
• Model types: physical (sand tanks), analogous (electric
field) and mathematical (governing equation)
Mathematical model: can be solved analytically and
numerically
Analytical model is too simplistic
Mathematical model is easy applicable and allows to solve
more complex flow problems
What is a model?
• Prediction:
most groundwater modeling efforts - predicting the
consequences of a proposed action
• Interpretation:
insight into the controlling parameters in a site-specific
setting, framework for formulating ideas about system
dynamics
• Generic geologic settings: helpful in formulating
regional regulatory quidelines and screening tool to
identify regions suitable for some proposed action
Why model?
Numerical model is a set of equations that describes the
physical and/or chemical processes occurring in a system.
_________________________________________________
Groundwater flow: determination of filtration velocity vector q (flow
direction and velocity) and hydraulic head h in modeled area
Solute transport: flow velocities
______________________________________________________
Components of numerical model:
Darcy’s low + Law of mass balance = Groundwater flow equation
Groundwater flow equation
3D Vector form
of Darcy law
Darcy’s law
hq
L
hh
kAQ 21
zzyzxz
yzyyxy
xzxyxx
kkk
kkk
kkk
zyx ),,(
.
,
,
z
h
k
y
h
k
x
h
kq
z
h
k
y
h
k
x
h
kq
z
h
k
y
h
k
x
h
kq
zzyzxzz
yzyyxyy
xzxyxxx
Hydraulic conductivity k is a Tensor
(9 components in anisotropic media)
Specific discharge q is a vector (3
components - anizotropic environment)
Hydraulic head h is a
scalar (1 component)
z
h
kq
y
h
kq
x
h
kq
zz
yy
xx
L
hh
kq 21
gradhkq .
Divided by flow crosssectional
area A
Groundwater flow equation
Itemized form of hydraulic
conductivity tensor
Hydraulic gradient
at anisostropic
media
Negative sign because fluids
flows from high pressure to
low pressure, If the change in
pressure is negative (where h1
> h2) then the flow will be in
the positive x direction
Law of Mass Balance = water balance equation
W
t
h
S
z
q
y
q
x
q zyx )()()(
Steady State Water Balance Equation
Inflow = Outflow
Transient Water Balance Equation
Outflow - Inflow = Change in Storage
Specific yield = Sy Storativity = S
W
t
h
Sq )(
Groundwater flow equation
Specific discharge q at anizotropic environment with
flow of compressible fluid, Transient regime
Enters or
leaves of
groundwater
Diveregence – scalar product of Hamilton operator and vector field (divergence is
a vector operator that measures the magnitude of a vector field's source or sink at
a given point, in terms of a signed scalar)
Darcy’s law + Law of mass balance = Groundwater flow equation
W
t
h
Sh )(hq W
t
h
Sq )(+ =
Equation of transient flow of
compressible fluid in anisotropic
and compressible environment
Simplification of Groundwater flow equation based on assumptions that:
• Tensor κ replaced by one-directional hydraulic conductivity k
• Water density ρ is constant and thus can be released
• Very often is possible to release transient term
• Sometimes may be z direction neglected
t
h
S
Groundwater flow equation
W
t
h
S
z
h
k
zy
h
k
yx
h
k
x
szzyyxx
Incompressible fluid in anisotropic incompressible medium
W
t
h
Sk
z
h
y
h
x
h
s2
2
2
2
2
2
Incompressible fluid in isotropic incompressible medium
W
z
h
k
zy
h
k
yx
h
k
x
zzyyxx
Incompressible fluid in anisotropic incompressible medium
Wk
z
h
y
h
x
h
2
2
2
2
2
2
Incompressible fluid in isotropic incompressible medium
W
t
h
S
y
h
bk
yx
h
bk
x
syyxx
W
t
h
Sk
y
h
b
yx
h
b
x
s
Wk
y
h
bk
yx
h
bk
x
yyxx
Wk
y
h
b
yx
h
b
x
If W=0 → Laplace equation
Groundwater flow equation
Transient flow
Steady state flow
3D 2D
3D
_____________________________________________________________________
2D
3D 2D
3D 2D
If is Analytical solution difficult – approximation by numerical methods
The finite differences method (FDM)
• Problem is solved by difference equations, which are discretized by nodes in grid
• Derivation is in nodes of grid replaced by differences – linear combination of functional
value (hydraulic heads, flow rates) in surrounding nodes
• System of linear algebraic equations with unknown values of differences in nodes
• System of equations is solved in district Ω bounded by boundary conditions
Finite difference method
Mesh centered cell-centered
Boundary Conditions
• Dirichlet boundary condition (boundary condition of 1. type): it
specifies the values of solution on the boundary of the domain,
specified head h=k
• Neumann boundary condition: it specifies the flow rate values
on the boundary of the domain, specified flow, mostly q=0
(boundary condition of 2. type)
• Cauchy boundary condition: flow rate is dependent on hydraulic
head on boundary, head dependent flow h=f(h)
• (boundary condition of 3. type)
• Initial conditions: in case of transient flow - hydraulic heads in
time=0
Groundwater flow equation + boundary conditions = numerical groundwater flow model
Finite difference method
Finite difference method
Cell i j k and indices for the six adjacent cells
4
1,1,,1,1
,
m
ji
m
ji
m
ji
m
jim
ji
hhhh
h
4
1
1,1,
1
,1,11
,
m
ji
m
ji
m
ji
m
jim
ji
hhhh
h
Numerical methods - iteration:
Problem is solved by repeating a process with the aim of approaching a desired goal. Each repettion is called an iteration and the
result of one iteration is used as the starting point for the next iteration.
• The iteration process starts from intial estimation (or measured).
• Calculation process proceeds for example In Gauss-Seidel interation from upper left corner to right and from up to down.
• In every node is calculated new value according to values in adjacent nodes
• New values are also calculated both new and old values in adjacent nodes
• Iterations proceed until reaching of Convergence criteria ε (i. e. required accuracy )
Gauss-Seidel Iteration
m
jih,jix ,
0
h
m
jih,
1
,
m
jih
2-D
3-D
initial guesses m
and next step m+1
1mm
hh
1mm
hh
initial guesses m
Iteration reached Convergence criteria
Finite difference method
Numerical methods – iteration:
Example of spreadsheet modeling
Governing groundwater flow equation – transient prediction of hydraulic head in 3D
domain for an anizothropy hydraulic conducticity field
Kxx, Kyy, Kzz = hydraulic conductivity
Ss = specific storage
W = flux term (sources or sinks)
h = hydraulic head
x, y, z = space coordinates
t = time
t
h
SW
z
h
K
zy
h
K
yx
h
K
x
szzyyxx
Governing flow equation is replaced by algebraic equations
written and solved for aech node of subdivided modeled region:
Finite differences:
Rectangular grid – cells
mesh-centered, block-centered
Finite elements:
Triangular elements – more accurate
diskretization
Basic principles of numerical modeling
Varies approach to model:
determinical - detailed knowledge of geology
stochastic – statistical evaluation
Model aplication:
conceptual model of groundwater flow – idealized
aquifer geometry
model as predictive tool –calibration, verification
Simulation:
groundwater flow – steady-state, transient
transport model – based on ground-water flowing model
Ground-water flow media:
pore media – the most sophisticated method, less problems
with hydraulic anisotropy
fractured media – difficulties coupled with fractured net
estimation
Numerical model conception
Conceptual model is a pictoral representation of
the gorundwater flow system:
Block diagram
Cross section
Simplification of the field system, complete
reconstrucion is not feasible
Conceptual model
Steps in formulating the conceptual model:
• Define the area of interest – identify the boundaries of
the model
• Definining hydostratigraphic units – comprise geologic
units of similar hydrogeologic properties (regional x local
scale), facies model, structural model
• Prepare a water budget – sources of water to the system
(recharge: precipitation, overland flow, from surface
water bodies), expected flow direction, outflows
(springflow, baseflow to stream, ET, pumuping)
• Defining the flow system – hydrologic and geochemical
information are used to conceptualize the movement of
groundwater through system - the location of recharge
and discharge areas
Conceptual model
Vztah mezi srážkami a vydatnostmi pramenů
0.01
0.10
1.00
10.00
1.1.2001 1.4.2001 30.6.2001 28.9.2001 27.12.2001
datum
vydatnost(l.s-1)
0.00
5.00
10.00
15.00
20.00
25.00
30.00
srážky(mm)
PB0290 PB0189 PB0487 Strážek Sejřek
Physical framework
• Geologic map and cross sections
• Topographic map – surface water bodies and divides
• Contour maps – elevation of the base of the aquifers and confining beds
• Isopach maps – thickness of aquifer and confining beds
• Map – extent and thickness of stream and lake sediments
Hydrogeologic framework
• Water table and potentimetric maps for all aquifers
• Hydrographs of groundwater head and surface water levels and discharge rates
• Maps and cross sections – hydraulic conductivity and/or transmissivity
distribution
• Hydraulic conductivity – distribution for stream and lake sediments
• Spatial and temporal distribution of: ET, Groundwater recharge, surface
water-groundwater interaction, groundwater pumping and natural groundwater
discharge
• Maps – morphology, geology, hydrogeology, climatic, landuse, settlement
• People activity - pumping or injection wells, build-up area
Data requirements for a Groundwater Flow Model
• Conceptual model determinates the dimension of the numerical model and
the design of the grid
• Grid – continous problem domain si discretized to finite difference blocks – cells
Finite difference grids: block centered x mesh-centered
Finite elements – more flexibility in designing a grid
• Defining model layers – hydrostratigraphic units (various conditions for
groundwater flowing)
• Orienting the grid
colinear with the principal direction of the hydraulic conductivity tensor
minimalize the number of cells outside the boundaries of the modeled area
grid colinear with regional flow direction
• Size of the nodal horizontal and vertical spacing – function of the
potentiometric surface and water table curvature, significant gradients requier finer
nodal spacing
• Grid refinement (regular and irregular grid) – the rule of thumb is to
expand the grid by increasing nodal spacing no more than 1,5 times the previous nodal
spacing
Coarse grid
Hydrogeological model: grid design
• Physical boundaries – formed by the physical presence
impermeable body of rock (outcrops), river, lake
• Hydraulic boundaries – streamline, groundwater divide
• Internal boundaries – sources and sinks within interior
of the grid (river, drain, lake, well etc.)
Hydrogeologic boundary condition:
I. type - specified head boundary, head = constant, often draining
stream
II. type - specified flow boundary, q = constant, often inflow outflow
boundary, q = 0, no flow boundary, often
groundwater divide, impermeable rock outcrop, impermeable
fault
III. type – head dependent flow boundary, q = f (h), mixed
function of stream
Correct selection of boundary condition si a critical step in model design
Hydrogeological model: boundaries
Water enters or leaves a model through the boundaries or
through sources and sinks
Sources and Sinks
• Injection and pumping wells
• Flux across the water table –
recharge x discharge (ET)
• Leakage – movement of water through a
layer that has a vertical hydraulic conductivity
lower than that of the aquifer. Leakage may
enter or leave tle aquifer depending on the
relative differences in heads between tle
aquifer and the source reservoir – river, lake –
III. type
Hydrogeological model: sources and sinks
Steady-state flow simulation
• Numerical iterations are calculated until the heads
and flow rates are unvarying
• Model calibration
Transient flow simulation
• Time dependent problems – transient variation of
water level, flow rates
• Begin with steady-state intial conditions
• More complicated for data manegment:
1. storage characteristics (storage coeffiecient for
confined aquifer and specific yield for unconfined
aquifer)
2. Intial conditions giving the head distributions in the
aquifer at the beginning of the simulation must be
specified in addition to boundary conditions
3. The time dimension must be discretized
Hydrogeological model: simulation flow type
Code selection:
• Finite elements, finite differences
• Saturated, unsaturated flow
• Transport – advective, multiphase
• Heat flow, fluid density
(complex packages – Fefflow, GMS, PMWIN)
Initiating Model Execution:
• Error criterion – nonconvergence (iteration residual error)
- inappropriate intial guesses of heads and parameters,
poorly discretizated time and space
• Error criterion for heads – convergence criteria
• Error criterion for tha water balance – checking of tle
residual error in tle solution allows total simulated inflows
and outflows as computed by tle Water Budget
Hydrogeological model: model execution
• Boundary conditions always affect a steady state
solution.
• Initial conditions should be selected to represent a
steady state configuration of heads.
• In general, the accuracy of the numerical solution
improves with smaller grid spacing and smaller time
step
• A water balance should always be included in the
simulation.
Modeling rules/guidelines
Calibrated model means that is capable to produce field-measured
heads and flows.
Finding of the right parameters (hydraulic conductivity, specific storage), boundary
conditions and stresses (recharge):
• First calibration (steady-state flow) and second calibration (transient flow)
• Mean annual water level or seasonal average of heads
Calibration values (Sample Information)
• Heads – errors caused by measure accuracy, scalling effect (average head in well with
the long screen, but model require point value), unmodeled heterogeneity, interpolations
errors (calibration values do not coincide with nodes)
• Fluxes – baseflow, springflow, infiltration from a losing stream or ET from the water table
– usually larger errors than errors asociated with head measurements
• It as advisable to use both heads and fluxes
Hydrogeological model: calibration process
Parameter estimates (Prior information) – used during the calibration
process
• Hydraulic conductivity and/or Transmissivity, storage parameters – usually
derived from aquifer tests
• Recharge – usually not available information – choose plausible rang of values
• Boundary conditions – uncertainty if do not correspond to the natural physical
boundaries of the aquifer, specified head boundary (III. type) will help achieve calibration
Calibration Techniques
• Manual trial and error adjustment of parameters – parameter values adjusted in
sequential model runs to match simulated heads and flows to the calibration targets
• Automated parameter estimation – specially developed codes guarantee the
statistically best solution and quantify the uncertainty in parameter estimates
Hydrogeological model: calibration process
Evaluating the Calibration
• Qualitative measure - comparsion between contour maps of
measured nad simulated heads
• Quantitative masure – scatter plot of measured against
simulated contours
• Deviation of points from
the straight line should be
randomly distributed
• If a trend is evident,
parameter values should be
adjusted
Hydrogeological model: calibration process
Expressing the average difference between
simulated and measured heads:
1. The mean error (ME) - the mean difference
between measured heads (hm) and simulated
heads (hs)
2. The mean absolute error (MAE) - the
mean for the absolute value of tle differences
in measured and simulated heads
3. The root mean squared (RMS) - error or
the standart deviation as the average of the
squared differences in measures and
simulated heads
n
i
ism hhnMAE
1
/1
i
n
i
sm hhnME
1
/1
5,0
1 2
/1
i
n
i
sm hhnRMS
Hydrogeological model: calibration process
Sensitivitiy analysis - quantify the uncertainity in the calibrated model
caused by uncertainity in the estimates of aquifer parameters
Changing two or more parameters also
might be examined to determine the
widest range of plausible solutions
Model verification - help establish
greater confidence in the calibration
Calibrated values are used to simulate a
transient response for which a set of data
exists (e.g. Record of groundwater levels
in response to drought or long-term
pumping)
Hydrogeological model: sensitivity analysis and model
verification
• Trace of out flow paths or pathlines by
tracking the movement of infinitely
small imaginary particles placed in the
flow field
• Help visualize the flow field and to track
contaminant paths
• Detect conceptual errors
Transport of particles by advection
gradhnKv e/
Hydrogeological model: particle tracking of
groundwater flow
• Unsuitable type of the numerical model – finite
elements x differences
• Deficit of important data – heads (lack of the
well), flow rates, precipitations (without gauge
stations) → poor calibration
• Heterogenity – hydraulic conductivity, character of
the heterogenities: linear heterogeneity (fractures)
or areal (rock bodies form)
• Drying and wetting of the cells
• Deep flow – recharge, depth of the groundwater
circulation – boundary condition
Hydrogeological model: problems
References
Anderson, M. P., Woessner, W. W. (1992): Applied groundwater modeling:
simulation of flow and advective transport.- Academic Press, San Diego, CA.
Cvachovcová, E. (2009): Matematický model proudění podzemní vody.- MS,
Diplomová práce, PřF MU, Brno.
Chiang, W-H., Kinzelbach, W. (2001): 3D-Geoundwater Modeling with PMWIN –
A Simulation System for Modeling Groundwater Flow and Pollution.- Springer.
Michael, G., McDonald, G., Harbaugh, A., W. (1988): A modular threedimensional
finite-difference ground-water flow model.- U.S. Geological Survey,
Open-File Report 83-875.
Říha, J. et al. (1997): Matematické modelování hydrodynamických a disperzních
jevů.- VUT Brno.