Introduction of MODFLOW
Yangxiao Zhou PhD Associate Professor in Hydrogeology UNESCO-IHE Institute for Water Education
General Information
• Developed by McDonald & Harbaugh of the USGS, 1983
• Public Domain available, 1988
• Block-Centered, 3D, modular structural, finite difference groundwater flow model • Most widely used groundwater model • Steady state or transient saturated flow • Currently MODFLOW 2000, 2005 available
Mathematical model
3D groundwater flow in the saturated heterogeneous and anisotropic porous media
h h h h K xx K yy K zz W Ss x x y y z z t where: Kxx, Kyy, Kzz = values of hydraulic conductivity along xyz axes (LT-1) h = total head (L) W = Sources and sinks (T-1) Ss = Specific storage (L-1) t = time (T)
Finite Difference Model
Model grid:
• Rows, columns, layers
• Count from the upper left corner
Finite Difference Model
Model grid:
• Cell-Centered
• Horizontal grid: rectangular cell varies with size
Finite Difference Model
Model grid:
• Cell-Centered
• Vertical layers: thickness varies
MODFLOW Packages
• Basic • Flow Package – BCF – LPF – HUF • River • Drain • Well • Recharge • Evapotranspiration • Changing Head Boundary • General Head Boundary • Horizontal Flow Barrier • Stream-Aquifer Interaction • Solvers • Output Control
MODFLOW Basic Package (BAS)
Functions: • Discretization of model domain; • Selection of major options and the designation of their input unit numbers; • Specifying initial and boundary conditions, and • Discretization of time.
MODFLOW Basic Package (BAS)
Definition of model grids and boundaries
• Number of rows, columns and layers • Cells outside model domain are marked as inactive • Controlled through IBOUND array
MODFLOW Basic Package (BAS)
Definition of stress periods and time steps
• Simulation time is divided into stress periods • The values of transient stresses are constant at each stress period • Each stress period is further subdivided into time steps.
MODFLOW Basic Package (BAS)
Initial heads
• • • • Provides initial set of heads for iterative solver Defines head values at specified head cells Provide initial conditions for transient modelling Typical approaches – Assign a constant value – Interpolate from scatter points of measured values – Use head solution from previous simulation
MODFLOW Flow Package (BCF)
• Each layer is assigned a layer type – Type 0: Confined layer: transmissivity and storage coefficient are constant. – Type 1: Unconfined layer: transmissivity varies, specific yield is used to calculated the change of the storage. – Type 2: Confined/Unconfined: transmissivity remains constant, storage parameter may alternates between storage coefficient and specific yield. – Type 3: Confined/Unconfined: transmissivity varies, storage parameter may alternates between storage coefficient and specific yield.
South North Layer 1 2 3 4 5
MODFLOW Flow Package (BCF)
• Layer data are entered depending on the type
– – – – – – Hydraulic conductivity K Top and bottom elevation Transmissivity Leakance Effective porosity Specific yield and storage
MODFLOW Recharge Package (RCH)
Areal recharge:
QRi,j = Ii,j*DELRj*DELCi
QRi,j is the recharge flow rate to the cell (i,j) expressed as a fluid volume per unit time;
Ii,j is the recharge flux to the cell (i,j) in units of length per unit time; can be constant or transient; The units of length and time are consistent with the units used in all other model parameters.
MODFLOW Recharge Package (RCH)
Options of application of recharge to model layers
MODFLOW ET Package (ET)
Evapotranspiratoion rate:
QET = QETM QET = QETM (h - hd)/d QET = 0 h > hs hd < h < hs h < hd
h: groundwater level in the cell (i,j) hs: ET surface QETM : maximum ET rate (length per unit time) when h > hs d: ET distinction depth hd = hs – d
MODFLOW ET Package (ET)
• Lump sum of E and T • Linear function • QETM •d
MODFLOW River Package (RIV)
Discretization of a river into reaches
MODFLOW River Package (RIV)
Cross section of a streamaquifer system (a) and Conceptual representation of stream-aquifer interconnection in simulation (b)
• If head is above river stage, flow is from aquifer to river • If head is below river stage, flow is from river to aquifer
MODFLOW River Package (RIV)
Exchange flux between the river and aquifer:
QRIV = CRIV*(HRIV-hi,j,k), hi,j,k > RBOT QRIV = CRIV*(HRIV-RBOT), hi,j,k < RBOT QRIV HRIV is the flow between the stream and the aquifer; is the stream level;
hi,j,k is the groundwater head in the cell underlying the stream reach; RBOT is the elevation of the bottom of the streambed layer; CRIV =K*L*W/M, is the hydraulic conductance of stream-aquifer; K is the conductivity of the streambed, L is the length of the stream reach, W is the width of the stream reach, and M is the thickness of the streambed.
MODFLOW River Package (RIV)
• Suitable for permanent rivers • Exchange of water will not cause large change of river stages
MODFLOW River Package (RIV)
River conductance:
L = Length of reach
Direction of flow K = Hydraulic conductivity of river bed material
M= Thickness of river bed W = Width of river
K * ( area of flow ) KLW C RIV ( length of flow ) M
MODFLOW River Package (RIV)
Case 1: Groundwater level is above the river stage:
Hijk
HRIV
RBOT
QRIV
CRIV
QRIV = CRIV * (HRIV - Hijk) Negative QRIV indicates that flow is from aquifer to river
MODFLOW River Package (RIV)
Case 2: Groundwater level is below the river stage:
HRIV
Hijk
CRIV
RBOT
QRIV
QRIV = CRIV * (HRIV - Hijk)
Positive QRIV indicates that flow is from river to aquifer
MODFLOW River Package (RIV)
Case 3: Groundwater level is below the river bottom:
HRIV
RBOT
CRIV
QRIV
H ijk
QRIV = CRIV * (HRIV - RBOT) Positive QRIV indicates that flow is from river to aquifer
MODFLOW Drain Package (DRN)
Underground drain (a) and open drain (b)
Used to simulate • Agricultural drains • Springs • Creek beds
MODFLOW Drain Package (DRN)
Simulation of drainage rate:
QDi,j,k = CDi,j,k*(hi,j,k - di,j,k) for hi,j,k > di,j,k
QDi,j,k = 0
for hi,j,k < di,j,k
CDi,j,k is a lumped (or equivalent) conductance describing all of head losses mentioned above; di,j,k is the drain elevation.
MODFLOW Drain Package (DRN)
• Discharge by springs and drains • Spring discharge measurements should be used to calibrate conductance
MODFLOW Well Package (WEL)
• • • • • Assigned to individual cells Q can be steady state or transient Extraction well (negative Q) Injection well (positive Q) Well package is used to simulate individual well rate, or well field rate • Well package is often used to simulate flow boundaries
MODFLOW Well Package (WEL)
The discharge of the multi-layer well:
Ql/Qw = Tl/ T where: Ql is the discharge from layer l to a multi-layer well in a given stress period; Qw is the total discharge of the multi-layer well in that stress period; Tl is the transmissivity of the layer l; and T is the sum of the transmissivities of all layers penetrated by the wells.
MODFLOW GHB Package (GHB)
• Assigned to individual cells • Used to simulate head dependent flow boundaries and lakes • Required parameters – Head hbi,j,k – Conductance Cbi,j,k
Qbi,j,k = Cbi,j,k * (hbi,j,k - hi,j,k)
MODFLOW GHB Package (GHB)
Lake Lake bottom sediments Cell
Flow Direction
Area = A
Ccell
Thickness = L
KA L
MODFLOW GHB Package (GHB)
• Use of simulating 3rd boundary conditions • Assume continuous hydraulic contact • Exchange of water between aquifer and sources/sinks
MODFLOW HFB Package (GHB)
Horizontal Flow Barrier
• Not part of original MODFLOW packages • Used to simulate low permeability barriers such as sheet pile walls, impermeable fault, or slurry trenches • Assigned to cell boundaries • Each instance is assigned a hydraulic characteristic = K/thickness
Plan View
Adjacent Cells
Cell Face
MODFLOW HFB Package (GHB)
Horizontal Flow Barrier
• Not part of original MODFLOW packages • Used to simulate low permeability barriers such as sheet pile walls, impermeable fault, or slurry trenches • Assigned to cell boundaries • Each instance is assigned a hydraulic characteristic = K/thickness
Plan View
Adjacent Cells
Cell Face
MODFLOW Stream Package (SFR)
MODFLOW Stream package:
The sources of inflow include a specified inflow at the beginning of the first reach of any segment (Qsri); the sum of tributary flow from upstream segments into the first reach of a segment (Qtrb); direct overland runoff to a reach (Qro), precipitation that falls directly on a reach (Qppt); and ground-water leakage to a reach calculated by the model (QLi).
MODFLOW Stream Package (SFR)
MODFLOW Stream package:
These losses include streamflow out of a reach (Qsro); specified diversions from the last reach in a segment (Qdiv); evapotranspiration from a reach (Qet); and leakage to the underlying aquifer (QLo).
For each reach, the sum of flows (in units of volume per time) into the reach is equal to the sum of flows out of the reach :
MODFLOW Stream Package (SFR)
MODFLOW Stream package:
Flow at the midpoint is computed as:
Use Manning’s equation to determine depth as a function of flow::
C is a constant, which is 1.0 for units of cubic meters per second or 1.486 for units of cubic feet per second; n is Manning’s roughness coefficient, dimensionless; A is cross-sectional area of the stream, units of length squared; R is hydraulic radius of the stream, in units of length; and S0 is slope of the stream channel, in units of length per length. w is width of channel, in units of length; and y is depth of water in stream, in units of length
MODFLOW Stream Package (SFR)
Irregular river reach Iterative method
MODFLOW Stream Package (SFR)
Test case: Semi-area Recharge from river leakage Steady simulation Transient simulation
MODFLOW Stream Package (SFR)
Test case: Semi-area Recharge from river leakage Steady simulation Transient simulation
MODFLOW Stream Package (SFR)
MODFLOW Stream Package (SFR)
MODFLOW Stream Package (SFR)
MODFLOW Stream Package (SFR)
MODFLOW Stream Package (SFR)
Input and output files for new STR package for MODFLOW 2000 Input: *.str *.gag Output: *.flw *.ag1 *.dvag9
MODFLOW Input Files
Name Basic Abbr. BAS Description Handles a number of administrative tasks, such as selection of packages, specification of boundaries, time-step length, initial conditions, and printing of results. Formulates finite-difference equations representing flow within porous medium Adds flow by wells to the finite difference equations. Adds flow by areally distributed recharge to the equations. Adds flow by rivers to the finite difference equations. Adds flow by drains to the finite difference equations Adds flow by ET to the finite difference equations. Adds flow by general-head boundaries to the finite-difference equations. Solves the system of finite-difference equations iteratively using the Strongly Implicit Procedure. Solves the system of finite-difference equations iteratively using Successive Overrelaxation. Input files BAS.DAT
BlockCentered Flow Well Recharge River Drain Evapotranspiration General- Head Boundary Strongly Implicit Procedure Slice Successive Overrelaxation
BCF WEL RCH RIV DRN EVT GHB SIP
BCF.DAT WEL.DAT RCH.DAT RIV.DAT DRN.DAT EVT.DAT GHB.DAT SIP.DAT
SSOR
SOR.DAT
MODFLOW Output Files
File name OUTPUT.DAT Format ASCII Description Information on inputs, running errors and water budget. Calculated heads at all grids for every time steps. Calculated drawdowns at all grids for every time steps.
HEADS.DAT DDOWN.DAT
Binary Binary
BUDGET.DAT
Binary
Calculated cell-by-cell flow components
MODFLOW Output Files
Visualisation of outputs with PMWIN
-
Contour maps of groundwater heads
Contour maps of drawdowns
-
Hydrograph
Animation of water table development Water budget of defined zones 