MODFLOW
MODFLOW is a modular, three-dimensional finite-difference groundwater flow model developed by the United States Geological Survey (USGS) for simulating steady-state and transient groundwater conditions in aquifers, and it is widely regarded as an international standard for such hydrologic modeling.[1] First released publicly in 1984, MODFLOW has evolved through multiple versions to address increasingly complex groundwater systems, with its modular design allowing users to customize simulations by selecting and combining packages for various hydrologic stresses such as wells, rivers, recharge, and evapotranspiration.[1] Key historical updates include MODFLOW-2005, which enhanced solver capabilities and added support for nonlinear flow equations in unconfined aquifers via the MODFLOW-NWT variant released in 2011.[1] The latest core version, MODFLOW 6, introduced in 2017, represents a major overhaul as an object-oriented framework that supports multiple model types within a single simulation, including groundwater flow (GWF), solute transport (GWT), energy transport (GWE), and particle tracking (PRT).[2] This version improves upon predecessors by enabling unstructured grids (such as triangular or quadrilateral elements), nested grid capabilities, and seamless integration of surface-water and groundwater interactions, while including tools for converting older models to the new format.[2] MODFLOW's versatility has made it essential for applications in water resource management, environmental impact assessments, and scientific research worldwide, with ongoing development ensuring compatibility with emerging computational needs.[1]Overview
Description and Purpose
MODFLOW is a modular, three-dimensional finite-difference groundwater flow model developed by the United States Geological Survey (USGS) in Fortran, designed to solve the groundwater flow equation for simulating flow in saturated porous media.[3][4] The model simulates both steady-state and transient groundwater flow conditions in confined, unconfined, and mixed aquifers, enabling analyses of hydraulic heads, flow paths, and interactions with surface water features.[3] Its primary purposes include supporting groundwater resource management for water supply planning, assessing contaminant transport and plume migration, and predicting land subsidence due to aquifer compaction.[5][6] A hallmark of MODFLOW is its modular design, which allows users to selectively activate or deactivate simulation processes—such as river-aquifer interactions via the River (RIV) package or well pumping via the Well (WEL) package—without altering the core code, facilitating customized simulations for diverse hydrogeologic settings.[3][7] The software operates on both Windows and Unix-based systems and has been freely available in the public domain since its initial release, promoting accessibility for researchers, consultants, and agencies worldwide.[3] MODFLOW has become the de facto international standard for groundwater modeling, widely adopted by hydrogeologists and integrated into geographic information systems (GIS) and other hydrologic software tools like ModelMuse and FloPy for enhanced preprocessing, simulation, and visualization.[1] As of 2025, the latest iteration, MODFLOW 6, incorporates support for unstructured grids to improve resolution in complex geometries and provides an application programming interface (API) for custom integrations, such as coupling with energy transport models.[2][8]Development History
MODFLOW was conceived in 1981 by Michael G. McDonald and Allan W. Harbaugh at the U.S. Geological Survey (USGS) as a modular three-dimensional finite-difference groundwater flow model, with the initial version completed and released in 1984 as USGS Open-File Report 83-875.[9] Written in Fortran 66, the model drew inspiration from modular designs in petroleum reservoir simulation and earlier hydrologic models to enable flexible assembly of simulation components for diverse groundwater problems.[9] This development addressed the growing need for a standardized, user-friendly tool for groundwater modeling, spurred by 1980s regulatory demands for water resource management and environmental protection, including assessments related to hazardous waste sites under programs like Superfund.[1] Key early milestones included the 1988 release of MODFLOW-88, which updated the code to Fortran 77 and provided revised documentation in USGS Techniques and Methods, book 6, chapter A1, enhancing accessibility and public adoption.[9] In 1996, MODFLOW-96 introduced a name file for improved usability and significant Fortran code refinements, broadening its application within and beyond the USGS.[9] The 2000 version, MODFLOW-2000, incorporated the "process" concept to support solute transport and parameter estimation, along with sensitivity analysis capabilities, marking a shift toward more integrated hydrologic simulations.[9] This was followed in 2005 by MODFLOW-2005, which adopted Fortran 90 with modular data management, advanced solvers, and support for nested and refined grids to handle complex spatial variability.[9] A modern paradigm shift occurred with the 2017 launch of MODFLOW 6, an object-oriented framework in Fortran 2003 that supports unstructured grids, multiple models in a single simulation, and robust solvers like Newton-Raphson formulations, building on prior unstructured grid efforts in MODFLOW-USG.[9] In 2020, the MODFLOW One-Water Hydrologic Flow Model (MF-OWHM) version 2 was released, extending the MODFLOW-2005 framework to integrate surface water, groundwater, and agricultural processes for conjunctive-use analysis in one-water hydrology contexts.[10] Under continued USGS leadership, development has emphasized extensibility, with ongoing efforts as of 2025 focusing on API integrations based on the Basic Model Interface (BMI) for enhanced interoperability and climate-related modeling enhancements, such as the new Groundwater Energy Transport (GWE) model for heat simulation released in May 2024.[2][11]Theoretical Basis
Groundwater Flow Equation
The groundwater flow equation in MODFLOW is derived from the fundamental principles of Darcy's law and the continuity equation for fluid flow in porous media. Darcy's law states that the specific discharge \mathbf{q} is proportional to the hydraulic gradient, expressed as \mathbf{q} = - \mathbf{K} \nabla h, where \mathbf{K} is the hydraulic conductivity tensor and h is the hydraulic head.[12] The continuity equation ensures mass balance for incompressible flow, stating that the divergence of the specific discharge equals the rate of change in storage minus sources or sinks: \nabla \cdot \mathbf{q} = S_s \frac{\partial h}{\partial t} - W, where S_s is the specific storage and W represents volumetric fluid sources or sinks per unit volume.[9] Combining these yields the three-dimensional groundwater flow equation for saturated, anisotropic aquifers.[12] For transient flow in confined aquifers, the governing partial differential equation (PDE) is: \frac{\partial}{\partial x} \left( K_{xx} \frac{\partial h}{\partial x} \right) + \frac{\partial}{\partial y} \left( K_{yy} \frac{\partial h}{\partial y} \right) + \frac{\partial}{\partial z} \left( K_{zz} \frac{\partial h}{\partial z} \right) + W = S_s \frac{\partial h}{\partial t} Here, h is the hydraulic head (in units of length L); K_{xx}, K_{yy}, and K_{zz} are the principal components of hydraulic conductivity along the x, y, and z directions (L/T); W is the source/sink term (T^{-1}); S_s is the specific storage (L^{-1}); and t is time (T).[12][9] This form assumes constant saturated thickness and applies to heterogeneous and anisotropic conditions, with the coordinate axes aligned to the principal directions of anisotropy.[12] Variants of the equation address different aquifer types and flow regimes. In unconfined aquifers, the equation incorporates variable saturated thickness, using specific yield S_y (dimensionless) instead of S_s when the water table falls below the top of the model cell, to account for drainage under the water table; this introduces nonlinearity, often handled via approximations like the Newton-Raphson method.[9] For steady-state conditions, the transient term is removed (\frac{\partial h}{\partial t} = 0), simplifying to: \frac{\partial}{\partial x} \left( K_{xx} \frac{\partial h}{\partial x} \right) + \frac{\partial}{\partial y} \left( K_{yy} \frac{\partial h}{\partial y} \right) + \frac{\partial}{\partial z} \left( K_{zz} \frac{\partial h}{\partial z} \right) + W = 0 with storage terms omitted.[12][9] Boundary conditions define the domain edges and interactions. Dirichlet conditions impose fixed hydraulic head (h = constant), such as constant-head cells representing lakes or rivers with known water levels.[12] Neumann conditions specify flux, including no-flow boundaries (zero flux) or fixed influx like recharge.[9] Cauchy conditions handle head-dependent flux, such as Q = C (h - h_b), where Q is flow rate (L^3/T), C is conductance (L^3/T), h is aquifer head, and h_b is boundary head, applied to features like streams or drains.[12][9] The equation assumes saturated flow throughout the domain, with constant fluid density and viscosity, which limits its direct applicability to scenarios like saltwater intrusion.[12] These assumptions are relaxed in extensions such as SEAWAT, a MODFLOW-based program that couples variable-density flow with solute transport.[13]Numerical Methods
MODFLOW employs the finite-difference method to approximate solutions to the groundwater flow equation on a discretized domain. In early versions, such as the original 1984 model and MODFLOW-2005, the domain is divided into a block-centered grid consisting of rows, columns, and layers, where each cell is centered on the block and represents a control volume. Spatial derivatives are approximated using central differences; for example, the partial derivative of hydraulic head with respect to the x-direction is computed as \partial h / \partial x \approx (h_{i+1,j,k} - h_{i-1,j,k}) / (2 \Delta x), enabling the calculation of inter-cell flows based on Darcy's law.[14][7] Grid discretization has evolved to accommodate complex geometries. Early implementations relied on structured rectangular grids, which are efficient for regular domains but less flexible for irregular boundaries. In MODFLOW 6, introduced in 2017, the model supports unstructured grids, including Voronoi tessellations via the DISU package, allowing cells with arbitrary polygonal shapes and user-defined connectivity to better represent heterogeneous aquifers and intricate boundaries while maintaining mass conservation through a control-volume finite-difference (CVFD) formulation.[7][9] Time discretization in MODFLOW uses an implicit finite-difference scheme to ensure numerical stability, particularly for transient simulations. The backward Euler method is applied, where the change in hydraulic head over a time step \Delta t is approximated as (h^n - h^{n-1}) / \Delta t, with all terms evaluated at the current time level n, allowing larger time steps without stability restrictions.[7][9] The discretized equations are assembled into a linear system of the form A \mathbf{h} = \mathbf{Q}, where A is the conductance matrix representing hydraulic connections between cells, \mathbf{h} is the vector of hydraulic heads, and \mathbf{Q} includes source and sink terms such as recharge, pumping, and storage changes. Conductances are computed from aquifer properties like hydraulic conductivity and grid spacing, with iterative solvers addressing the sparse symmetric matrix.[7][9] For multi-aquifer systems, vertical flow between layers is handled through leakage terms based on vertical conductance. In block-centered formulations, the flow Q between adjacent layers is given by Q = C (h_{k+1} - h_k), where C is the vertical conductance derived from vertical hydraulic conductivity and layer thicknesses, facilitating simulation of confining beds and inter-aquifer exchange. In MODFLOW 6, this extends to unstructured grids with generalized connectivity.[7][9]Assumptions and Limitations
MODFLOW operates under several fundamental assumptions that simplify the representation of groundwater systems for numerical simulation. It assumes homogeneous (uniform within each cell) but potentially anisotropic hydraulic properties, such as hydraulic conductivity, within each finite-difference cell, allowing for uniform flow calculations across the cell despite potential subsurface heterogeneity.[9] The model is designed exclusively for saturated groundwater flow, neglecting unsaturated zone dynamics and capillary fringe effects, with water release or uptake from storage assumed to occur instantaneously.[9] Additionally, MODFLOW presumes constant fluid properties, including density and viscosity, which restricts its direct application to scenarios without these invariances.[9] It does not natively support multiphase or unsaturated flow processes, focusing instead on single-phase, saturated conditions.[9] These assumptions impose notable limitations on MODFLOW's applicability. The model adheres strictly to Darcy's law for flow between cells, rendering it unsuitable for turbulent or non-Darcian flows, such as those in high-velocity karst systems.[9] Its grid-based finite-difference discretization can limit resolution for sharp hydraulic gradients, potentially introducing errors in areas requiring fine-scale detail, as accuracy depends on cell size and alignment with principal conductivity directions.[9] Large-scale simulations are computationally intensive, particularly with nonlinearities, transient conditions, or complex packages, though later versions incorporate parallelization to mitigate this.[9] Common challenges include numerical oscillations in unconfined aquifers, which can arise during wetting and drying cycles or with Newton-Raphson formulations, often addressed through specialized solvers like the Newton formulation.[9] MODFLOW lacks built-in support for density-dependent flow, necessitating extensions such as SEAWAT for variable-density simulations.[13] To overcome some resolution limitations, workarounds like nested grids via the Local Grid Refinement (LGR) utility in MODFLOW-2005 allow embedding finer child grids within coarser parent grids, enabling targeted high-resolution modeling in areas of interest without excessive computational cost across the entire domain.[15] Extensions for advanced physics, such as those integrating solute transport, further expand capabilities beyond core saturated flow.[13] Model reliability is highly sensitive to parameterization, where small changes in inputs like hydraulic conductivity can significantly alter simulated heads and fluxes, underscoring the need for rigorous calibration against field data such as observed hydraulic heads and streamflows to ensure accuracy and reduce non-uniqueness.[16]Version History
Early Versions (1980s–1990s)
The development of MODFLOW began in the early 1980s at the U.S. Geological Survey (USGS), with initial work on a modular three-dimensional finite-difference groundwater flow model conducted between 1981 and 1983 by Michael G. McDonald and Arlen W. Harbaugh.[9] This prototype, coded primarily in Fortran 66, focused on basic simulation of saturated groundwater flow using block-centered finite differences on structured grids, without full modularity in its earliest form.[17] The model was first publicly documented in 1984 as USGS Open-File Report 83-875, establishing the foundational framework for solving the groundwater flow equation through finite-difference approximations.[17] In 1988, the USGS released MODFLOW-88, a significant rewrite in Fortran 77 that introduced the model's hallmark modular structure, allowing independent packages to handle specific hydrologic processes such as wells, rivers, and recharge. Core packages included the Basic (BAS) package for input/output management and the Block-Centered Flow (BCF) package for simulating internal flow between cells using conductance calculated via Darcy's law and the harmonic mean method. This version supported both steady-state and transient simulations, with capabilities for confined, unconfined, and convertible aquifers, though limited to a single instance of each package.[9] Incremental updates through the early 1990s, such as version 2.6, added improved error handling and numerical stability enhancements, including preconditioned conjugate-gradient solvers to address convergence issues in large models.[18] By the mid-1990s, MODFLOW had evolved further with the release of MODFLOW-96 in 1996, which refined input/output handling through a standardized name file and updated documentation to enhance user accessibility. This version incorporated early integrations like the MODPATH particle-tracking program, first developed in 1989 as a post-processor for pathline analysis in steady-state flow fields.[9] Additional packages emerged during this period, such as the Streamflow-Routing (STR1) package for stream-aquifer interactions introduced in 1989 and the Interbed-Drainage package for subsidence modeling in 1991, expanding the model's ability to simulate surface-groundwater exchanges and geomechanical effects.[9] Initial experiments with graphical user interfaces, like the MODFLOW-GUI developed for version 96, began to address the command-line limitations of earlier releases.[19] MODFLOW's early versions standardized block-centered finite-difference methods, enabling efficient discretization of heterogeneous aquifers on rectangular grids, which became a cornerstone for groundwater modeling. By the early 1990s, the model had gained widespread adoption within the USGS and external agencies, including early Environmental Protection Agency (EPA) assessments of contaminant plumes and aquifer management, due to its flexibility and open-source availability.[9] However, these releases faced limitations, such as reliance on symmetric solvers without parallel processing support, restrictions to structured grids, and challenges in handling wetting/drying in unconfined systems, which constrained simulations of complex transient behaviors.[9]MODFLOW 2000–2005
MODFLOW-2000, released in July 2000 by the U.S. Geological Survey (USGS), represented a significant evolution of the modular groundwater flow model, introducing an enhanced structure that facilitated the integration of advanced capabilities such as parameter estimation and solute transport simulation. This version maintained the finite-difference formulation of prior iterations but expanded the modular framework to include "processes" alongside traditional packages, allowing for more flexible model configuration and broader applicability to complex hydrologic systems. The initial release (version 1.0) was followed by iterative updates, with version 1.19.01 in March 2010 providing compiled executables for Microsoft Windows and incorporating bug fixes for stability in large simulations.[20] Key enhancements in MODFLOW-2000 focused on improving computational efficiency and model calibration. It incorporated advanced linear equation solvers, including the Preconditioned Conjugate-Gradient (PCG2) package for faster convergence in iterative solutions and the Geometric Multigrid (GMG) solver, which leveraged multigrid techniques to handle systems with millions of nodes more effectively. Additionally, an interface for external algebraic multigrid solvers like SAMG was enabled through the Link-AMG (LMG) package, further accelerating solutions for large-scale problems.[21] For parameter estimation, MODFLOW-2000 introduced built-in processes akin to PEST functionality, including the Sensitivity Process (SEN1), Observation Process (OBS1), and Parameter-Estimation Process (PES1), which allowed users to optimize model parameters against field data directly within the simulation workflow. The Link-MT3DMS (LMT6) package was also added to seamlessly couple MODFLOW-2000 with the MT3DMS solute transport model, enabling multi-species mass transport simulations essential for contaminant studies. MODFLOW-2005, released in 2006 and supported through 2017, built upon the MODFLOW-2000 framework with refinements aimed at handling heterogeneous aquifers and multi-scale simulations.[22] It introduced support for running multiple models within a single execution, particularly through integration with locally refined grid (LGR) capabilities, which allowed embedding high-resolution child grids within coarser parent grids to capture fine-scale features without excessive computational demand.[18] Version 1.12.00, released in 2017, enhanced the Hydrogeologic-Unit Flow (HUF) package for representing heterogeneous aquifer units with vertically varying properties, improving accuracy in layered systems.[23] For unconfined flow simulations, MODFLOW-2005 laid the groundwork for the Upstream Weighting (UPW) package, which was later formalized in extensions like MODFLOW-NWT to better manage nonlinearities in wetting and drying cells through continuous head-dependent functions. Among the key innovations during this period, the Observation package (OBS1 in MODFLOW-2000 and refined in 2005) enabled systematic comparison of simulated heads, fluxes, and advective transport against field observations, supporting robust model calibration and uncertainty analysis.[24] Basic parallel processing was also introduced, allowing MODFLOW-2000 and 2005 to be compiled with Message Passing Interface (MPI) for distributed computing on multiprocessor systems, which reduced runtime for large models by distributing matrix solving tasks.[25] These developments from 2000 to 2005 were driven by the growing demand for simulating regional basin-scale groundwater systems, where increased model sizes and integration with surface water processes required more efficient solvers and flexible parameterization to inform water resource management decisions.[18]MODFLOW 6 and Later Developments
MODFLOW 6, released in 2017 by the U.S. Geological Survey (USGS), represents a major redesign of the MODFLOW framework using an object-oriented Fortran structure to enable flexible simulation of groundwater flow and related processes.[26] This version incorporates capabilities from previous variants, supporting both structured and unstructured grids through the Discretization by Layers (DIS) package for regular grids and the Discretization by Vertices (DISV) package for irregular, vertex-based grids, allowing for more accurate representation of complex geologic features.[26] Additionally, MODFLOW 6 integrates the Groundwater Transport (GWT) model for simulating solute transport in three dimensions, facilitating coupled flow and transport analyses within a unified framework. Key enhancements in MODFLOW 6 versions include the addition of an application programming interface (API) in version 6.2.2, released in July 2021, which allows external programs to interact with the model without source code modifications, enhancing customization and integration with other software.[2] Subsequent releases, such as version 6.6.0 in December 2024 and 6.6.3 in 2025, have focused on bug fixes, performance optimizations, and expanded support for advanced solvers. As of November 2025, the latest version is 6.6.3.[27] MODFLOW 6 builds on precursor developments, including MODFLOW-NWT from 2011, which introduced a Newton-Raphson formulation to handle nonlinearities in unconfined flow and drying/rewetting cells more robustly.[28] This was followed by MODFLOW-USG in 2016, which provided the foundational unstructured grid capabilities later integrated into MODFLOW 6 for simulating tightly coupled processes on flexible meshes.[29] As of November 2025, MODFLOW 6 development continues actively on GitHub, with ongoing contributions from the USGS and community to refine core functionalities and add new packages.[1] Relatedly, the MODFLOW One-Water Hydrologic Flow Model (OWHM) version 2.3.1b-4, released in February 2025, bolsters coupled simulations of farm water use, river-aquifer interactions, and conjunctive management scenarios. Among its advantages, MODFLOW 6 supports multi-model simulations within a single run, enabling seamless coupling of flow, transport, and surface processes for comprehensive hydrologic assessments.[30] Python bindings via the FloPy library facilitate model setup, execution, and visualization, broadening accessibility for researchers and practitioners. The framework also maintains backward compatibility with earlier MODFLOW input formats, easing transitions from legacy models.[2]Model Components
Core Packages
The core packages in MODFLOW form the foundational components for simulating groundwater flow, defining the model domain, hydraulic properties, and essential boundary conditions. These packages are modular, allowing users to activate them selectively through input files to construct simulations tailored to specific hydrogeologic settings. All MODFLOW simulations require at least the Basic Package and a flow package, with boundary packages added as needed to represent external stresses. This modular design, introduced in the original MODFLOW version in 1984, enables flexible configuration while ensuring compatibility across versions.[31] The Basic Package (BAS in early versions, BAS6 in MODFLOW 6) establishes the fundamental structure of the model, including the spatial grid, temporal discretization, initial conditions, and output specifications. It defines the number of layers, rows, and columns; the lengths of stress periods and time steps within each; and the status of each cell (active, inactive, or constant-head). BAS6 also handles the calculation and reporting of hydraulic heads, drawdowns, and water budgets, providing essential outputs for model verification. Originally introduced in 1984 as part of the inaugural MODFLOW release, the Basic Package has evolved to support both structured and unstructured grids in MODFLOW 6 (released 2017), while maintaining backward compatibility for legacy models.[31][9] Flow packages govern the internal representation of groundwater movement through the aquifer, specifying properties such as hydraulic conductivity, storage coefficients, and wetting/drying behavior. The Block-Centered Flow Package (BCF or BCF6 in later adaptations) was the original flow package, introduced in 1984, which uses a block-centered finite-difference scheme to compute inter-cell conductances for confined, unconfined, or convertible layers. It calculates horizontal and vertical flow terms based on layer transmissivities, thicknesses, and anisotropy, supporting dewatering corrections for unconfined conditions. In MODFLOW 6, the Newtonian Flow Package (NPF) supersedes BCF and other legacy flow packages (e.g., Layer-Property Flow or LPF), providing enhanced capabilities for storage variations, hydraulic conductivity assignments (isotropic or anisotropic), and wetting thresholds to reactivate dry cells. NPF integrates optional Newton-Raphson linearization for handling nonlinear water-table dynamics and supports variable vertical conductance options like dewatered or perched conditions. These packages are activated once per model and are essential for all flow simulations.[31][9] Boundary packages simulate interactions between the groundwater system and external features, introducing head-dependent or fixed fluxes at model edges or internal points. The Well Package (WEL) represents extraction or injection at discrete locations, specifying flow rates independent of local heads, with options for automatic reduction if drawdown limits are exceeded. Rivers (RIV), introduced in 1984, model head-dependent leakage using river stage, bed conductance, and bottom elevation, ceasing flow when the water table drops below the riverbed. The Drain Package (DRN), also from 1984, simulates surface drains or tiles that remove water only when heads exceed the drain elevation, based on conductance. General-Head Boundaries (GHB), originating in 1984, allow flexible head-dependent fluxes to external reservoirs by specifying reference heads and conductances. Constant-Head boundaries (CHD), defined via the Basic Package's cell status array since 1984 (with dedicated input enhancements by 1991), enforce fixed potentials at specified cells, influencing adjacent flows without direct rate specification. Later additions include the Lake Package (LAK), introduced in 2005 for MODFLOW-2005 to simulate dynamic lake levels, storage, and interactions with aquifers via connected cells, and the Multi-Aquifer Well Package (MAW) in MODFLOW 6 (2017), which extends WEL to handle wells penetrating multiple layers with variable flow distribution. These packages are invoked as needed per stress period through input files, contributing terms to the finite-difference equations for budget tracking.[31][32][9]| Package | Key Function | Head Dependency | Introduced |
|---|---|---|---|
| WEL | Pumping/injection at points | None (specified rates) | 1984[31] |
| RIV | River-aquifer leakage | Yes (conductance-based) | 1984[31] |
| DRN | Drainage when head > elevation | Yes (conductance-based) | 1984[31] |
| GHB | Flux to external heads | Yes (conductance-based) | 1984[31] |
| CHD | Fixed potentials | None (enforced heads) | 1988 |
| LAK | Dynamic lake storage and flow | Yes (level-dependent) | 2005 |
| MAW | Multi-layer well flows | None (specified rates, distributed) | 2017[9] |