An Introduction to MODFLOW

1. An Introduction to MODFLOW

2. 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. • MODFLOW escaped! • 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) • Latest version: MODFLOW-2005 (Harbaugh, 2005)

3. What is MODFLOW? • Cited in statutes, legally tested • 250,000 hits on Google for MODFLOW • 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) • Phreeqc connections (PHT3D) (Henning Prommer) • Commercial • MODHMS-Surfact (Integrated sw/gw/unsat) • GUI’s: Visual MODFLOW, Groundwater Vistas, GMS, PMWin, …

4. 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? • Latest release of USGS MODFLOW • 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 from UCODE_2005

5. 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…???

6. MODFLOW-2005 Processes • Processeseach solve a fundamental equation. Of importance in this class are • Ground-water Flow (GWF) Kh = S(h/t) … • Observation (OBS) y = y′ + e

7. GWF Packages • Packages each represent a type of system feature. Of importance in this class are • Package that defines model layers and properties: • Layer-Property Flow (LPF) Package • Packages used to add/remove water at a specified rate: • Well (WEL) • Recharge (RCH) • Packages that add/remove water based on head in the aquifer: • General-Head Boundary (GHB) • River (RIV)

8. How Processes and Packages Interact In MODFLOW, subroutines are named using the three-letter identifiers for processes and packages. For example, GWF1LPF6RP

9. Back to the world of users instead of programmers --What is required for a simulation? • 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 (here, LPF) • Solver. Here we use • Preconditioned Conjugate Gradient (PCG)

10. Activating capabilities • Turn Packages on and define input files using the NAME file • Example: # GW Flow process input files bas6 41 tc1.bas lpf 42 tc1.lpf wel 43 tc1.wel pcg 44 ../data/tc1.pcg . . .

11. Basics of Data Input • Array data • Data input in arrays with one row for each row of the model grid and one column for each column of the model grid. Sometimes repeat one array for each model layer. # # # # # ……… # # # # # ……… # # # # # ……… # # # # # ……… . . . • List data • Data input using lists of cells layer row column …… Example: 1 3 43 ……. 2 62 53 …….

12. What is MODFLOW? MODFLOW is a calculation program Input files (plain text or binary) Output files (plain text or binary)

13. Oftenuse MODFLOW through a (Graphical) User Interface Results Maps Model

14. MODFLOW capabilities used in class • Class exercise • MODFLOW Packages used • Layer-Property Flow (LPF) • Recharge (RCH) • River (RIV) • General-Head Boundary (GHB) • Advective Transport (ADV) • Preconditioned-Conjugate Gradient (PCG) • MODFLOW Processes used • Ground-Water Flow (GWF) • Observations (OBS) • UCODE_2005 capabilities used • Sensitivity • Parameter-Estimation

15. 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 Here, describe selected aspects of capabilities used in class

16. Head-dependent boundaries Generally use many cells to define a feature. Here, shaded cells are used to simulate flow to compare to measured flow Q2-Q1. Other cells would be used to define the rest of the river. From Hill+, 2000

17. A Head-dependent boundaries For each finite-difference cell n: Qn = (KA/M)n (Hn – hn) • (KA/M)n = Cn = conductance of assumed distinct streambed • Often define Cn with parameters • Cn=FnP1 • Additive: Cn= Fn1P1 + Fn2P2 • Hn= water-body stage • hn= simulated head Areal view of typical cell n Cross-section of typical cell n Cell center

18. GHB EXPLANATION 3 q the simulated flow rate at one cell (L /T) n q (A) n (nega tive for flow out of the ground - water Positive q system) n indicates K the hydraulic conductivity (L/T) of, for n Slope = - Cn = - (K A )/D flow into n n n example, the riverbed or lakebed the subsurface D the thickness (L) of, for example, the riverbed n or lakebed A the area of the water body within the finite - n h n q = 0 2 n difference cell (L ) C the conductance calculated using K , D , and n n n A . n h is the simulated hydraulic head in the ground - n Negative q n water system adjacent to the head - dependent indicates boundary (L); and flow out of the H is the water level in the water body or the n (C) subsurface H elevation of the drain (L) n E is the botto m of the streambed n q n (C) Positive q (B) n indicates q n Slope = - Cn = - (K A )/D flow into n n n h the n q = 0 n subsurface Slope = - Cn = - (K A )/D n n n h n q = 0 n Negative q n Negative q indicates n indicates flow out of H flow out of n the the subsurface subsurface E H n n Packages that represent head-dependent boundaries GHB • q=C(H-h) • Important here: • GHB: General-Head Boundary • RIV: River RIV RIV DRN

19. RIV Package with hn below RBOTn (hn<RBOTn) Cell center

20. Pumping wells • Well (WEL) Package • List input: layer, row , column, rate (negative means flow out of the ground-water system) • Rate can be defined using parameters • Problem: If a well intersects many model layers, how much water comes from each layer?

21. 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.

22. Parameters • In the MODFLOW model for the class problem, parameters are used to define the following model inputs • Layer-Property Flow (LPF) Package • Horizontal hydraulic conductivity of model layers (HK) • Vertical hydraulic conductivity of an implicit confining unit (VKCB) • Recharge (RCH) Package • Recharge rate (RCH) • River (RIV) Package • Riverbed conductance (RIV) • Values of defined parameters can be controlled using the PVAL file. This makes it easy for users of UCODE_2005, etc.

23. C of river bed equals the value in the package input file times the factor in the package input file. Here, the factor is 1000. C of river bed = 1000 × 0.0012 The factor can be different for different cells. Parameters  Model Input • River Package input file PARAMETER 1 18 18 54 MXACTR IRIVCB K_RB RIV 1.200000E-03 18 1 1 1 100. 1000. 90. . . . 1 18 1 100. 1000. 90. 0 1 ITMP NP -- Stress Period 1 K_RB

24. C of river bed equals the value in PVAL times the factor in the package input file. Here, the factor is 1000. C of river bed = 1000 × 0.00116999 The factor can be different for different cells. Parameters and PVAL file • River Package input file • PVAL file PARAMETER 1 18 18 54 MXACTR IRIVCB K_RB RIV 1.200000E-03 18 1 1 1 100. 1000. 90. . . . 1 18 1 100. 1000. 90. 0 1 ITMP NP -- Stress Period 1 K_RB 6 HK_2 1.523554700000E-5 HK_1 4.619000000000E-4 VK_CB 9.903220000000E-8 K_RB.0011699900000000 RCH_2 38.39840000000000 RCH_1 47.55430000000000 Easy to use UCODE_2005 to change parameter values in PVAL file.

25. Flow model creation: Observations MODFLOW’s Observation Process allows simulated values to be compared to observations. • Here, use it for the following observations • Head observations • At a cell • Changes in head over time • Flow observations • Over the reach of a feature represented by the RIV Package

26. Observations 10 0 0 50 1.E+30 NH,MOBS,MAXM,IUHOBSV,HOBDRY 1. TOMULTH (below, lay,r,c,ts,roff,coff,toff,obs) hd01.ss 1 3 1 1 0.0 0.0 0.0 101.80 hd02.ss 1 4 4 1 0.0 0.0 0.0 128.12 hd03.ss 1 10 9 1 0.0 0.0 0.0 156.68 . . . hd10.ss 2 18 6 1 0.0 0.0 0.0 142.02 • Heads (HOB in name file) • River Gain (RVOB in name file) • Output file (“data 50 ex8._os” in name file) 1 18 1 50 NQxx,NQCxx,NQTx 1.00000E+00 TOMULTxx 1 18 NQOBxx,NQCLxx flow01.ss 1 0.0 -4.4 ts,toff,obs 1 1 1 1.00 lay,r,c,factor . . . 1 18 1 1.00 "SIMULATED EQUIVALENT" "OBSERVED VALUE" "OBSERVATION NAME" 100.209701538086 101.800003051758 hd01.ss 126.954444885254 128.119995117188 hd02.ss . . . -4.41627883911133 -4.40000009536743 flow01.ss

27. Constructing input files • In class we will either use 00-MFI2005.exe or the files will be constructed already. • Instructions for using 00_MFI2005.bat are provided in class.

28. Execute MODFLOW • Here, we will use 00-MFI2005.exe or already constructed batch files. • Detailed instructions are provided in the exercise instructions. • Basically, need to provide the name file filename on the same line (this is often done in a batch file) MODFLOW test.nam

29. Model results • Possible results for class problem (depends on options chosen) • Global budget (check for overall solution accuracy) • Heads at each active cell in the grid at each time step • Flows at each cell face • Simulated equivalents to observations • Often use software to visualize results. In class, use ModelViewer

30. Pathline Modeling • Advective transport. Used here as a first investigation of transport predictions • Requires: • Flow solution • Porosity – to determine velocity • Starting locations • The particle tracking is calculated using MODPATH, which uses results produced by MODFLOW.