Modelling groundwater flow with

MODFLOW: an introduction
v.1.5 released on 20/09/2019

Laura Foglia1, Iacopo Borsi2, Giovanna De Filippis3, Rudy Rossetto3

1. Technical University Darmstadt (Germany) & UC Davis, California

2. TEA Sistemi S.p.A. – Pisa (Italy)
3. Scuola Superiore Sant’Anna – Pisa (Italy)
 DOCUMENTATION LICENSES

Please attribute FREEWAT

with a link to http://www.freewat.eu/

Except where otherwise noted,

this slides are licensed under a
Creative Commons Attribution 4.0 International License:
http://creativecommons.org/licenses/by/4.0/

FREEWAT - Free and Open Source Software Tools for Water Resource Management
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant
agreement No 642224
 DATA LICENSES

Please attribute FREEWAT with a link to

http://www.freewat.eu/

Except where otherwise noted,

data used in this tutorial are licensed under a
Open Database License:
http://opendatacommons.org/licenses/odbl/1.0/

Any rights in individual contents of the databases are

licensed under the Database Contents License:
http://opendatacommons.org/licenses/dbcl/1.0/
FREEWAT - Free and Open Source Software Tools for Water Resource Management
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant
agreement No 642224
 This lecture presents an introduction to
MODFLOW, highlights the important
features and mainly provides
references to the most useful literature
and webpages! It does not pretend to
be a complete MODFLOW manual!

FREEWAT - Free and Open Source Software Tools for Water Resource Management
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant
agreement No 642224
 Overview of MODFLOW

MODFLOW: public-domain three-dimensional finite-difference

groundwater simulation code by the U.S. Geological Survey.

MODFLOW can be run as a MS-DOS command line executable and

reads several ASCII text files as input describing grid, properties,
boundary conditions, initial conditions, time steps, stresses, solution
and output options. Output is to ASCII text and binary files.

A graphical user interface such as FREEWAT can provide images of

model features and can support visualization and graphical
manipulation of model input and output. A GUI can create the
MODFLOW input files and execute MODFLOW. You never need to see
the text files or know the commands that are necessary to run
MODFLOW ... until something goes wrong!

1
 MODFLOW
Documentation and reading

Read the MODFLOW documents in order

MODFLOW 1988 1996 2000 (GWProcess OBSProcess) 2005
All these and much more can be downloaded at
http://water.usgs.gov/software/lists/groundwater/

Read for the essence of how the codes perform.

You do not need to read detailed descriptions of input or source code
at this time, but can go back for those details when needed in the
future

Note the conceptual model features of MODFLOW discussed in the

1988 are very important and are NOT repeated in later documents

2
 MODFLOW
Documentation and reading

Further Reading for MODFLOW modelers

There are many additional simulation package and many data-

processing tools.

For example, Anderson and Woessner, 1992, Applied Groundwater

Modeling is an excellent textbook for groundwater-flow modeling

And the new 2015 edition!!!

3
 Keep up-to-date

New Developments and MODFLOW Software updates

USGS Software Web Site

http://water.usgs.gov/software/lists/groundwater/

BOOKMARK THAT PAGE!

Let’s open a new browser session there now

4
Scroll Down to
http://water.usgs.gov/ogw/modf
low/MODFLOW-2005-Guide/ 5
6
And the new MODFLOW-OWHM
One Water Hydrologic flow Model

7
 Online MODFLOW Manual

USGS Online Guide to MODFLOW-2000/2005

http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/guide.h
tml

ESSENTIAL BOOKMARK!

The online guide provides references to MODFLOW documents,

definitions and formats for variables, and many helpful tips.

**Note you can download this to your computer and use it off line,
but it will not be updated unless you download and replace it

8
9
 What is MODFLOW?
Widely used ground-water flow simulation program that runs on any
platform (Windows, Sun, Unix, Linux,…).

Mostly written in standard FORTRAN (GMG is C++)

Solves the ground-water flow equation with different possible

properties, boundary conditions, and initial conditions

First version, 1983, McDonald and Harbaugh. Written to serve USGS

needs. Education emphasized.

Public domain (Free)

Open source (Anyone can check and change the source code)
• Changed versions are sometimes commercial –
it’s up to the developer
Well documented - Modularly constructed (More later) 10
 What is MODFLOW?

Cited in statutes, legally tested

Many other programs use results from or are based on MODFLOW:

– Public domain/open source
• MT3DMS (multi-species solute or heat transport, some
reactions, dual porosity) (Chunmiao Zheng, U Alabama)
• MODPATH (particle tracking) (Dave Pollock, USGS)
• SEAWAT (density-dependent transport using MODFLOW and
MT3DMS) (Chris Langevin, USGS)
• PHT3D (= MT3DMS + PHREEQC) (Henning Prommer)
– Commercial
• MODHMS, MODFLOW-Surfact (Integrated sw/gw/unsat)
• GUI’s: PMwin, VisualMF, Groundwater Vistas, GMS, Triwaco, …

11
 What is MODFLOW?

Program organized into MODules

– Activate the capabilities you need; no overhead from other
capabilities (execution time, RAM)
– The structure is clear and documented for adding additional
capabilities such as new equations
– Modularity in ‘Processes’ and ‘Packages’

What is MODFLOW-2005?

 Internal computer storage redesigned to support storage of

multiple models necessary for local grid refinement and facilitate
linkages to other models (GSFLOW = MODFLOW + PRMS).
 Parameter-estimation, sensitivity analysis, uncertainty now done in
UCODE_2005
12
 Summary of MODFLOW capabilities
• Single Phase, Saturated Flow
• Constant Density
• Porous Media (Darcy's Law applies)
• 1, 2, Or 3 Dimensional (2D areal or cross-section OR Quasi-3D)
• Steady State or Transient Flow
• Heterogeneous, Anisotropic (aligned with grid), Layered Units
• Block-Centered Finite-Difference with variable grid spacing
• Layers can be Confined or Convertible
• Boundary conditions include: Dirichlet, Neuman, Cauchy, and
Phreatic Surface
• Boundary condition packages include Wells, Recharge,
Evapotranspiration,
• Rivers, Drains, Thin Barriers to Horizontal Flow, Re-wetting, and
others.
• Many Solver options
• Monitoring of simulated values comparable to field observations
• Sensitivity analysis and parameter estimation to be accomplished
by UCODE 13
 Selected MODFLOW capabilities
Transport
– Link-MT3DMS (LMT). Produces files for use by MT3DMS (Zheng)
– SEAWAT v4. Variable density using MODFLOW and MT3DMS.
Improvements in representation of viscosity, chemical diffusion, constant
head boundaries.
– Ground-Water Transport (GWT) Package (currently for MF2000 only)
Solvers
– Link-Multi-Grid (LMG). Typically up to 10 times faster than PCG2, but uses
a lot of RAM. Free to USGS users.
– Geometric Multi-Grid (GMG). Can be very fast!
– PCG with Improved Nonlinear Control (PCGN).
Packages for improved simulation of saturated-zone
processes
– Multi-Node well (MNW)
– Hydrogeologic Unit Flow (HUF) Package updates: LVDA: Model-layer
variable-direction horizontal anisotropy. KDEP: K varies with depth. SYTP:
partial accounting for free surface; numerically stable. 14
 Selected MODFLOW capabilities

Local Grid Refinement

– MODFLOW-LGR: Accurate local grid refinement. Can have multiple refined
areas. Coming soon: LGR with particles tracked between regional and
local grids.
Interactions with processes above the water table and at
land surface
– FARM Process: Sophisticated accounting for ET, diversions, etc related to
agricultural demands.
– Streamflow-Routing with unsaturated flow beneath streams (SFR2)
– Unsaturated zone using (1) kinematic wave (UZF) and (2) Richard’s
equation (VSF)
– GSFLOW: Integrated surface processes represented using PRMS; separate
program
Optimization for ground-water management
– Ground-Water Management (GWM) Process
15
 Who is MODFLOW?
Collaborative open-source development with roots at the USGS

Some USGS developers

– Arlen Harbaugh (MODFLOW, Reston, Virginia, USA)
– Ned Banta (MODFLOW-2000, Lakewood, Colorado, USA)
– Mary Hill (SA/PE/UA, MODFLOW-2000, UCODE, MMA, Boulder, Colorado, USA)
– Steffen Mehl (local grid refinement (LGR), SA/PE/UA, now at CalSU-Chico, USA)
– Stan Leake (compaction and subsidence, TMR, Tucson, Arizona, USA)
– John Hoffman (compaction and subsidence, TMR, Tucson, Arizona, USA)
– Dave Prudic (gw/sw interaction, STR, SFR, GSFLOW, Carson City, Nevada, USA)
– Rick Niswonger (gw/sw interactions, SFR, GSFLOW, Carson City, Nevada. USA)
– Paul Barlow (ground-water management, MODMAN, Reston, VA, USA)
– Randy Hanson (FARM Process, MNW, San Diego, USA)
– Alden Provost (HUF, Reston, VA)
– Dave Pollock (particle tracking, MODPATH, Reston, Virginia, USA)
– Chris Langevin (transport, saltwater intrusion, SEAWAT, Miami, Florida, USA)
– Lennie Konikow (transport extended from MOC3D, GWT, Reston, Virginia, USA)
– George Hornberger (transport extended from MOC3D, GWT, Reston, Virginia, USA)
Some non-USGS developers
– Chunmiao Zheng (transport, MT3DMS, University of Alabama, USA)
– Eileen Poeter (UCODE, MMA, Colorado School of Mines, IGWMC, Golden, CO, USA)
– Evan Anderman (ADV, HUF, now at EvanAnderman.com, photography)
– Henning Prommer (MODFLOW+PHREEQC, CSIRO, Perth, Australia)
– Wolfgang Schmid (FARM Process, U. of Arizona, USA)
– David Ahlfeld (ground-water management, GWM, U. of Massachusetts, USA)
– You…??? 16
 MODFLOW development

•Originally organized in modules

•Modules grouped into packages that perform calculations either

specific to the behavior of a geohydrologic feature or a numerical
modeling task

•Packages allow
•examination of specific hydrologic features independently
•facilitates development of additional capabilities

•Originally solely a ground-water flow model

17
 MODFLOW development

•Scope broadened to allow capabilities such as transport and

parameter estimation

•Expansion of modular design required (MODFLOW-2000)

•addition of “Process”

•MODFLOW-2005 is similar in design to MODFLOW-2000

•Incorporates different approach for managing internal data

•Fortran modules are used to declare data that can be shared
among subroutines
•MODFLOW subroutines were originally called modules
•Generic term module has been eliminated and replaced by the
term subroutine

18
 Original modular structure (1988)
•BAS - basic package
•general tasks - gridding, constant head and no-flow boundaries, initial conditions,
time stepping
•OC - output control package
•controls the information and format of results
•BCF - block centered flow package
•layer types, grid dimensions, material properties
•WEL - well package
•locations and flow rates of wells
•RCH - recharge package
•recharge rates and locations
•RIV - river package
•locations, river bed material properties, and river stages
•DRN - drain package
•location, material properties surrounding drains, and elevation of drains
•EVT - evapotranspiration package
•parameters describing evapotranspiration rate with depth to water table
•GHB - general head boundary package
•locations, local material properties, and elevation of specified heads
•SOLVERS
•SIP - strongly implicit procedure package
•SOR - slice-successive over-relaxation package
19
 Packages written after original MODFLOW
•PCG2 - preconditioned conjugate-gradient 2 package
•alternative matrix solver
•STR1 - stream routing package
•differs from the river package in that the surface water stage varies based on the
surface water flow and the Manning equation
•BCF2 - block-centered flow 2 package
•allows for re-wetting of cells that have gone dry
•BCF3 - block-centered flow 3 package
•a supplement to the BCF2 package, allowing alternative interblock transmissivity
formulations
•HFB1 - horizontal flow barrier package
•simulation of thin, vertical, low permeability features that impede horizontal flow
•TLK1 - transient leakage package
•simulates transient leakage and storage changes in confining units of quasi-3D
models
•GFD1 - general finite difference flow package
•substitutes for the BCF package, allows user to enter conductance rather than
calculating with MODFLOW
•IBS1 - interbed storage package
•simulates compaction of compressible, fine-grained units within or adjacent to
aquifers in response to pumping
•CHD1 - time-variant specified head boundary package
•allows time varying specified head 20
 Packages written for MODFLOW-2000
and since
•rapidly growing long list
•earlier packages are listed here,
•refer to the
•USGS MODFLOW and related programs web page, OR
•use the USGS OnLine Guide for MODFLOW

GWF1 - ground water flow process (GWF in name file) finite difference simulation
of saturated porous media flow
OBS1 - observation process (OBS in name file) monitors value of head or flow at
specified locations
SEN1 - sensitivity process (SEN in name file) calculates the change in simulated
head and flows at observation locations PES1
- parameter estimation process (PES in name file) estimates values of
parameters by nonlinear regression to minimize the weighted sum of
squared residuals for observations
DIS - discretization package (DIS in name file) gridding, defining division of space
and time for the numerical solution
MULT - multiplier file (MULT in name file) defines the spatial distribution of
multipliers in the grid that act on parameter values specified in those
zone 21
 Packages written for MODFLOW-2000
and since
ZON - zone file (ZONE in name file) defines the spatial distribution of zones in
the grid where specified parameters apply
BAS6 - basic package (BAS6 in name file) constant head and no-flow boundary
conditions; and initial conditions
OC - output control package (OC in name file) controls the information and
format of results
BCF6 - block centered flow package (BCF6 in name file) defines material
properties with some parameters being dependent on grid
dimensions (e.g. transmissivity), thus this package ignores the
discretization information in DIS for some purposes -- the parameter
method of inputting data cannot be used -- method of interblock
conductance calculations can be selected
LPF1 - layer property flow package (LPFin name file) an alternative to BCF6
defines material properties with all parameters independent of grid
dimensions (e.g. hydraulic conductivity) -- the parameter method of
inputting data can be used -- method of interblock conductance
calculations can be selected
HFB6 - horizontal flow barrier package (HBF6 in name file) represents thin
barriers that occur between model cells by defining their hydraulic
conductivity divided by their thickness and specifying where they
occur 22
 Packages written for MODFLOW-2000
and since

WEL6 - well package (WEL in name file) locations and flow rates of wells
RCH6 - recharge package (RCH in name file) recharge rates and locations
RIV6 - river package (RIV in name file) locations, river bed material properties,
and river stages
STR6 - stream routing package (STR in name file) differs from the river package in
that the surface water stage varies based on the surface water flow
(calculated as specified flow and ground water flux to/from stream) and
the Manning equation
DRN6 - drain package (DRN in name file) location, material properties
surrounding drains, and elevation of drains (this update allows a fraction
[0-1] of the drain outflow to be returned to a specified cell)
EVT6 - evapotranspiration package (EVT in name file) parameters describing
evapotranspiration rate with depth to water table
GHB6 - general head boundary package (GHB in name file) locations, local
material properties, and elevation of specified heads
CHD6 - time-variant specified head boundary package (CHD in name file) allows
time varying specified head
23
 Packages written for MODFLOW-2000
and since
SOLVERS (SIP SOR PCG DE4 LMG in name file)
SIP5 - strongly implicit procedure package
SOR5 - slice-successive over-relaxation package
PCG2 - preconditioned conjugate gradient package
DE45 - direct solution by alternating diagonal ordering package
LMG1 - multigrid solver speeds execution for large grids and high
degree of heterogeneity
ADV2 - advective transport observation package (ADV2 in name file) allows use
of travel time observations for parameter observations
RES1 - reservoir package (RES in name file) simulates leakage between
reservoir and aquifer as reservoir area changes in response to stage
changes
FHB1 - flow and head boundary package (FHB in name file) allows flow and
head boundary conditions that vary at times other than starting and
ending times of stress periods
IBS6 - interbed storage (subsidence) (IBS in name file) simulates compaction
related to hydraulic head decline
HUF1 - hydrologic-unit flow package (HUF in name file) calculates effective
hydraulic properties for cells based on geometric description of
hydrologic units 24
 Packages written for MODFLOW-2000
and since

LAK3 - lake package (LAK in name file) allows variation of lake stage based on
water budgets
ETS1 - evapotranspiration package with segment ET function (ETS in name file)
allows function describing evapotranspiration rate with depth to water
table to be piece-wise linear
DRT1 - drain package with return flows (DRT in name file) allows user to allocate
proportions of drain flow to be recharge to specified cells
LMT6 - link to MT3D (LMT in name file) allows printing of file to be read by
MT3D for contaminant transport
SFR - Streamflow-Routing package (SFR in name file) is used to simulate streams
in a model (provides greater flexibility in how streams are specified than
STR)
UZF - Unsaturated Zone Flow Package (UZF in name file) simulates vertical flow
of water through the unsaturated zone to the saturated zone

25
 Basic MODFLOW input and output
BASIC INPUT ITEMS INCLUDE:
Grid
Time stepping
Hydraulic parameters
Boundary Conditions
Stresses
Solution parameters

BASIC OUTPUT ITEMS:

Hydraulic Heads
Drawdown
Flow rates
Mass Balance
Optional info at specified times
Iteration information

Underscore files describing various model inputs/outputs

Binary files containing output in compressed form 26
 Main MODFLOW input file

NAME file (.nam)

The NAME file specifies the names of the input and output
files used for a MODFLOW simulation, associates each file
name with a FORTAN unit number, and identifies the
packages that will be used in the model.

An ASCII text file is needed for every process and package

defined in the MODFLOW NAME file. This class will focus
on the commonly used packages.

27
 MODFLOW basics:
the Finite-Difference method
For complex systems and boundary conditions, there are no analytical solutions
to the groundwater-flow equation, and numerical methods are needed.
MODFLOW uses finite-difference methods.

User
- DISCRETIZE modeled area into N grid cells
- Assign aquifer properties to each grid cell
- Assign boundary conditions

MODFLOW
- Set up finite-difference solution to the GROUNDWATER-FLOW
EQUATION for each grid cell (N EQUATIONS)
- Either H or Q will be KNOWN at each block N UNKNOWNS
- SOLVE for H or Q in all blocks using matrix algebra and
iterative solutions

  h    h    h  h
 bK x    bK y    bK z   Ss  W ( x, y, z, t )
x  x  y  y  z  z  t
28
 The Finite-Difference grid
Complex geologic material distributions are simplified.
Physical aquifer properties likely vary within and between layers.
Square or rectangular grid cells arranged in rows and columns.
Active (closed circle) and inactive cells (open circles).

Values assigned to grid cells

represent the geometry, aquifer
properties, boundary conditions,
initial condition, and stresses on a
groundwater system to build a
representation of field conditions.

---- Aquifer boundary Layers may correspond to

● Active cell horizontal geohydrologic
intervals
○ Inactive cell
Δrj Width of cell in row direction (j indicates column
number)
Δci Width of cell in column direction (i indicates row
number)
Δvk Thickness of the cell 29
Δr Δc Δv Volume of cell with coordinates (i,j,k)
 Finite-Difference cells

• Mesh-centered

• Block-centered
• Easier math for boundaries
• MODFLOW

30
 MODFLOW: defining layers
• Layer Representation Options
• Constant layer thickness /variable properties
• Expedites modeling
• Rough approximation
• Compatibility with another function

• Variable layer thickness /constant properties

• More representative of field conditions

31
 MODFLOW: defining layers
• One layer
• layer represents a single
hydrostratigraphic unit or aquifer

• Quasi-3D
• Hydrogeologic units horizontal
• Leakance

• Fully 3D
• Dipping units
• Aquifers and confining units explicit 32
 MODFLOW: grid orientation
• Grid drawn on an overlay of a map of the area to be
modeled

• If possible orient the grid so that the x and y axes are

colinear with Kx and Ky and vertical axis is aligned with
Kz

• For finite difference, try to minimize the number of

nodes that fall outside the boundaries of the modeled
area

• Set boundaries far from the area of interest so imposed

stresses to the interior of system don’t reach the
boundaries 33
 MODFLOW: discretize time

Discretize Time

• TIME STEPS: temporal equivalent of grid cells

• Small when stresses change and increase
in length to a constant, convenient size
until the stresses change

• STRESS PERIODS: groups of time steps during

which stresses do not change
• Temporal data compiled at these
increments
34
MODFLOW: discretize time

35
 MODFLOW: discretize time

• Difficult to decide on initial time step size

• MODFLOW requires the time period, number of steps
and a multiplier to gradually increase steps

Multiplier is typically 1.1 to 1.5

36
 Flow model creation: time

Steady state
– Inputs = outputs. No change in storage
– No time dimension: easier to visualize
– Errors in model setup more clear in results
Transient
– Requires (often steady-state) initial conditions
– Requires a value for storage
– Stresses are defined using stress periods (time interval of input)
– Each stress period is divided into time steps (time interval of head
calculation).

– Lengthy calculation times can produce large output files

– For some tips on when to “go transient”, see H.M.Haitjema (2006)

Role of Hand Calculations in Ground Water Flow Modeling,
Ground Water.
37
 How small is small enough?

• YOU KNOW YOUR DISCRETIZATION IS APPROPRIATE WHEN:

THE ANSWER REMAINS THE SAME FOR: SMALLER TIME STEPS,
STRESS PERIODS, AND SMALLER CELL SIZES

• TIME
• Easy to test smaller time steps
• Stress periods require recompiling stress data (may be time
consuming) and updating any packages with stresses specified

•SPATIAL
• Unless you have an automated grid generator / input file creator,
then the time requirements and logistics of rebuilding the model
with smaller cell sizes renders the task unreasonable
• Important to use smaller grid sizes from the beginning of
numerical model development because you will never be able to
test this issue.
• In reality, few if any modelers check this.
38
 Solving the Finite-Difference equations

The finite-difference equations for head are solved for each grid cell in
a model resulting in a matrix of equations to be solved.

A variety of matrix solution techniques are utilized to solve the finite-

difference equations. These techniques iteratively approximate the
solution. There are a few essential concepts in order to use a solution
wisely

DIRECT SOLVER – Exact solution of matrixes

Direct solutions used rarely because the problems we solve are large
and require direct solutions of the entire matrix with available
computer memory.

39
 Finite-Difference solution
ITERATIVE SOLVER – Iterative (approximate) solution of equations starting from
initial conditions

Iterative numerical solutions of the groundwater-flow equation matrix are

usually employed, and several options are available with MODFLOW.

Iterative solvers begin from an input initial condition (initial head at n=0) and
repeatedly approximate (iterate) the solution to the matrix, until there is very
little change from the previous approximation and the solution converges.

Convergence tolerance is often specified by the user. Sometimes results from

subsequent iterations change very slowly, and there is a need to speed up the
process. Other times the approximation causes large changes that repeatedly
overshoot the answer in opposite directions, and there is a need to slow down
the process.

Often the code provides the user the opportunity to specify an acceleration or
relaxation (two terms for the same variable) parameter to speed up or stabilize
the convergence process. Occasionally a particular solver requires specific input
items and you will need to carefully read the manual to understand the options,
but nearly all solvers require a tolerance and relaxation. 40
 Finite-Difference solution

CONVERGENCE IS "ARBITRARY" i.e. it depends on the tolerance value:

GENERALLY we assume that the user specified a reasonable tolerance

and that the solution will be fairly accurate (have a small mass balance)
at that precision level

For heads A REASONABLE TOLERANCE is small relative to the total

head change across the system and small relative to the head
difference between cells in low gradient areas.

41
 Finite-Difference solution

IF TOLERANCE IS VERY SMALL, LACK OF CONVERGENCE MAY NOT BE

SIGNIFICANT because the solution may be very precise, that is the mass
balance may be very small.

For example, a user may specify a tolerance for head change between
iterations that is smaller than the precision of the computer ... this
tolerance could never be met but the solution would be very precise. In
a case like this, it is usually desirable to increase the tolerance so that
the problem converges to avoid the unfounded concerns of those who
feel uncomfortable with the words "did not converge!" In reality the
solution will not be any better. One could simply look at the model
output to see how much head change occurred for the last iteration
and use a value just slightly larger. Then the solution would be identical
but the output would say "solution converged!"
42
 Finite-Difference solution

IF TOLERANCE IS VERY LARGE, PRESENCE OF CONVERGENCE MAY

NOT BE SIGNIFICANT because the solution may be very imprecise,
that is the mass balance may be large.

For example, a user may specify a tolerance for head change

between iterations that is large compared to the total head change
across the model or to the difference in head between cells in areas
of low gradient. The model would converge after a small number of
iterations but the flow field may appear erratic, even with flow
apparently generated out of "nowhere" because one cell ends with a
high head and flow moves out to all surrounding cells and there is no
source to "feed" that radial outflow. Such conditions result in a poor
mass balance.

43
 Acceleration or Relaxation
If progress toward convergence is proceeding too slowly, use a factor larger than
one. If progress toward convergence is over shooting the answer (maximum head
change alternates sign with each iteration), use a factor less than one.

Typically the maximum head change varies in location. A good code will report the
node at which the maximum occurred. If your model is not converging efficiently, it
can be useful to note where the maximum head change is occurring. Sometimes
evaluation of your input will reveal that you have a "typo", for example you may
have typed the exponent on a hydraulic conductivity incorrectly resulting in a large
contrast of conductivity between adjacent cells, or you may have added an extra
"0" to a boundary head for one cell causing a huge influx at one cell. 44
 Example excerpt from
BEGINNING OF SIMULATION MODFLOW output
SOLVING FOR HEAD
ITER: 0 RES: 2.5980E+07 CFAC: 1.000
ITER: 1 RES: 1.8517E+06 CFAC: 0.071
ITER: 2 RES: 4.9647E+05 CFAC: 0.268
ITER: 3 RES: 9.6388E+04 CFAC: 0.194
ITER: 4 RES: 2.0855E+04 CFAC: 0.216
ITER: 5 RES: 5.9308E+03 CFAC: 0.284
---------------------------------------------------------
PCG ITERATIONS : 5
DAMPING : 0.5662E-07
L2-NORM OF RESIDUAL : 5.9308E+03
MAX HEAD CHANGE : 1.7661E+07
MAX HEAD CHANGE AT (COL,ROW,LAY) : ( 24, 110, 12)
---------------------------------------------------------
..
.
.LATER TIME IN SIMULATION
ITER: 0 RES: 8.6464E+02 CFAC: 1.000
ITER: 1 RES: 7.0492E+01 CFAC: 0.082
---------------------------------------------------------
PCG ITERATIONS : 1
DAMPING : 0.1000E+01
L2-NORM OF RESIDUAL : 7.0492E+01
MAX HEAD CHANGE : 8.5938E-01
MAX HEAD CHANGE AT (COL,ROW,LAY) : ( 11, 49, 6)
---------------------------------------------------------

-------------------------------
TIME STEP : 5
STRESS PERIOD : 2
GMG CALLS : 4
PCG ITERATIONS : 4
------------------------------- 45
TOTAL PCG ITERATIONS : 17880
 Mass balance – another consideration
for convergence
ULTIMATELY THE SOLUTION MUST CONSERVE MASS
Although your conceptual model may be based on a balanced hydrologic
budget, the numerical model will have mass imbalances if the heads are
not calculated accurately. The fluxes will be based on gradients between
cells. If the heads are not accurate, the gradients will produce
inconsistent fluxes, thus the inflows and outflows (note that flow into
storage is an outflow and flow in from storage is an inflow) will not
balance.

EXAMPLES of MODELING PROBLEMS THAT CAN LEAD TO MASS

IMBALANCE INCLUDE:
• large convergence tolerance
• poor conceptual model formulation (i.e. set up is hydraulic nonsense)
• large contrasts in properties between adjacent cells
• large difference in size of adjacent cells
• large time steps or time steps so small that the defined flows are
below the accuracy of the computer solution
46
 Mass balance – another consideration
for convergence
MOST CODES PROVIDE A TABLE TO ASSESS MASS BALANCE

Generally the table includes a summary of inflow and outflow to and

from various types of features as well as the overall flow. This is
another good place to pause and assess the model results to be sure
everything makes sense. It is also valuable to evaluate the inflow or
outflow to each cell where a head or head dependent flux has been
specified. Such evaluation can reveal errors that could otherwise go
unnoticed yet have negative repercussions on results. The mantra for
this is repeated below, always monitor the flow at constant-head nodes
and the flux at specified-flux nodes.

A SMALL MASS BALANCE IS REQUIRED

substantially <1% (0.5% is the high end of acceptable)

47
 Example excerpt from
MODFLOW output
VOLUMETRIC BUDGET FOR ENTIRE MODEL AT END OF TIME STEP 5 IN STRESS PERIOD 16
------------------------------------------------------------------------------
CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T
------------------ ------------------------
IN: IN:
--- ---
STORAGE = 61384087785.7633 STORAGE = 5644891.4979
CONSTANT HEAD = 0.0000 CONSTANT HEAD= 0.0000
WELLS = 0.0000 WELLS= 0.0000
DRAINS = 0.0000 DRAINS= 0.0000
RIVER LEAKAGE = 105133899486.2039 RIVER LEAKAGE = 3559330.3534
HEAD DEP BOUNDS = 0.0000 HEAD DEP BOUNDS= 0.0000
RECHARGE = 2.7380E+12 RECHARGE = 7194487.6472
SPECIFIED FLOWS = 0.0000 SPECIFIED FLOWS = 0.0000
STREAM LEAKAGE = 69744176600.6383 STREAM LEAKAGE= 1955526.7305
ET SEGMENTS = 0.0000 ET SEGMENTS = 0.0000
MNW = 321347679.5662 MNW = 47991.3395
TOTAL IN = 2.9746E+12 TOTAL IN = 88402227.5684
OUT: OUT:
---- ----
STORAGE = 42534787126.9181 STORAGE = 729950.6762
CONSTANT HEAD= 0.0000 CONSTANT HEAD = 0.0000
WELLS = 40637515208.1547 WELLS = 17999510.6501
DRAINS = 16552965696.6484 DRAINS = 373047.8186
RIVER LEAKAGE= 13184784416.5005 RIVER LEAKAGE = 975905.1566
HEAD DEP BOUNDS=18640052988.0677 HEAD DEP BOUNDS = 399313.7562
RECHARGE= 0.0000 RECHARGE = 0.0000
SPECIFIED FLOWS= 55214439264.0000 SPECIFIED FLOWS = 1219104.0000
STREAM LEAKAGE=725689791735.0275 STREAM LEAKAGE = 7224961.1519
ET SEGMENTS = 1.8580E+12 ET SEGMENTS = 48937653.1254
MNW = 4123872601.1030 MNW = 541647.0605
TOTAL OUT = 2.9746E+12 TOTAL OUT = 88401093.3958
IN - OUT = 10027659.6802 IN - OUT = 1134.1726
PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY = 0.00 48
 Mass balance – another consideration
for convergence
SOMETIMES A LARGE MASS BALANCE IS NOT A PROBLEM:

For example, if there is a large imbalance for an early time step, perhaps even
200%, it may not be important, because models generally start with very small
time steps and gradually increase the length of each time step, so a 200% error
in the first time step may only reflect lost accounting for a gallon of water while
the problem as a whole involves 10s of thousands of gallons and subsequent
time steps have reasonable mass balances.

MOST IMPORTANTLY ONE MUST CAREFULLY EVALUATE THE MODEL OUTPUT TO

DETERMINE WHERE WATER ENTERS AND LEAVES THE MODEL AND ASSESS IF
THE FLOW PATTERN AND VOLUMES ARE REASONABLE

ALWAYS MONITOR:
Q @ CONSTANT H CELLS
AND
H @ SPECIFIED Q CELLS
49
 MODFLOW: defining
Boundary Conditions
• Correct selection of boundary conditions is a critical step in model design.

• For steady-state simulations, the boundary conditions largely determine the

resulting simulated flow patterns.

• For transient simulations, boundary conditions can effect simulation results

when transient stresses reach the boundary, and realistic boundaries should
be selected.

• Physical boundary examples and type:

• Impermeable rock (no-flow or specified-flux)
• Faults (no-flow or specified-flux)
• Lakes (constant-head), constant-flux, or head-dependent)
• Streams (constant-head, constant-flux, or head-dependent)
• Springs (constant-head, constant-flux, or head-dependent)
• Saltwater interface (no-flow)

• Hydraulic boundary examples and type:

• Groundwater divides (no-flow)
• Flow lines (no-flow) 50
 It is most desirable to terminate a model extent at physical hydrologic
boundaries.
However, model extents may need truncation to maintain the desired level of
detail and still have the model execute in a reasonable amount of time.
Consequently models sometimes have artificial boundaries. For example, heads
may be fixed at known water-table elevations at a county line, or a flow line or
groundwater divide may be set as a no-flow boundary.
BOUNDARY TYPE PHYSICAL EXAMPLES ARTIFICIAL USES
CONSTANT Fully Penetrating Surface- Distant Boundary (Line of
or Water Features unchanging hydraulic- head
SPECIFIED HEAD contour)
SPECIFIED FLUX Precipitation/Recharge Flow line
Pumping/Injection Wells Groundwater Divide
Impermeable material Subsurface Influx
HEAD DEPENDENT FLUX Rivers Distant Boundary (Line of
Springs (drains) unchanging hydraulic- head
Evapotranspiration contour)
Leakage From a Reservoir or
Adjacent Aquifer
51
 Types of Boundary Conditions

BOUNDARY BOUNDARY TYPE

MATHEMATICAL MATHEMATICAL
CONDITION &
NAME EXPRESSION
NAME OTHER NAME

Constant/
Type 1
Specified Dirichlet H=C
specified head
Head

No-Flow or Type 2
Neumann
Specified Flux specified flux

Head-
Type 3
dependent Cauchy Q = f(h, K, A)
mixed condition
Flux
After:
Definition of Boundary and Initial Conditions in the Analysis of
Saturated Ground-Water Flow Systems – An Introduction,
O. Lehn Franke, Thomas E. Reilly, and Gordon D. Bennett, USGS - TWRI Chapter B5, Book 3, 1987
And Anderson and Woessner, 1992, Applied Groundwater Modeling 52
 Back to the world of Users
(instead of programmers)

What is required for a simulation?

How do we start building the model?

Tell the program what capabilities to use
– Name file (NAM)

Package input files for each process (only the GWF Process is always required)
– Basic (BAS6) (can define constant head BC’s here)
– Discretization (DIS)
– Hydrogeologic info (for example, LPF or BCF)

Solver. Here we use

– For example, Preconditioned Conjugate Gradient (PCG) 53
 Activating capabilities

Turn Packages on and define input files using the NAME file

Example:
# GW Flow process input files
bas6 41 test0.bas
lpf 42 test0.lpf
wel 43 test0.wel
pcg 44 ../data/test0.pcg
.
.
 Using FREEWAT, the NAME file gets created automatically, but it is
important to know that the NAME is the place where the names of all
the files used in the model are stored.

54
 Basics of data input
Two types of input data: List data (for example as in the RIVER package
that is presented later) and Array data (for example as in the
RECHARGE package)

 List data  Array data

– Data input using lists  Data input in arrays with
of cells one row for each row of the
model grid and one column
layer row column …… for each column of the
model grid. Sometimes
Example: repeat one array for each
model layer.
1 3 43 …….
# # # # # ………
2 62 53 …….
# # # # # ………
# # # # # ………
# # # # # ………
55
 Before starting with the model
development…
…aspects of flow model creation

• Conceptual model
• Base map
• Grid design
– Areal
– Model layers (thickness can be variable)
• Boundary conditions
• Aquifer properties
• Pumping wells
• Recharge
• Time
56
 If you need any assistance, please contact
Laura Foglia, UC Davis (Davis, California) - lfoglia@ucdavis.edu
Iacopo Borsi , TEA Sistemi SpA – iacopo.borsi@tea-group.com
Giovanna De Filippis – Scuola Superiore Sant’Anna (Pisa - Italy) -
g.defilippis@santannapisa.it
Rudy Rossetto, Scuola Superiore Sant’Anna (Pisa - Italy) – r.rossetto@santannapisa.it
FREEWAT development received funding from the following projects:

1. Hydrological part was developed starting from the project SID&GRID, funded by Regione Toscana through EU POR-FSE 2007-2013
(sidgrid.isti.cnr.it) (2010-2013)
2. Porting of SID&GRID under QGIS has been performed through funds provided by Regione Toscana to Scuola Superiore S.Anna -
Project Evoluzione del sistema open source SID&GRID di elaborazione dei dati geografici vettoriali e raster per il porting negli
ambienti QGis e Spatialite in uso presso la Regione Toscana (CIG: ZA50E4058A) (2015)
3. Saturated zone solute transport simulation capability has been developed within the EU FP7-ENV-2013-WATER-INNO-DEMO
MARSOL. MARSOL project received funding from the European Union's Seventh Framework Programme for Research, Technological
Development and Demonstration under grant agreement n. 619120 (www.marsol.eu) (2014-2017)
4. FREEWAT was developed within the EU H2020 project FREEWAT - Free and Open Source Software Tools for Water Resource
Management. FREEWAT project received funding from the European Union’s Horizon 2020 research and innovation programme
under grant agreement n. 642224 (www.freewat.eu) (2015-2017)
5. Integration of SFT (StreamFlow Transport) and LKT (Lake Transport) packages of MT3D-USGS is being performed at Scuola
Superiore Sant'Anna within the project SMAQua (SMart ICT tools per l'utilizzo efficiente dell'AcQua) - co-financed by Regione Toscana,
ASA S.p.A. and ERM Italia S.p.A. (2018-2020)