Water Availability and Use Science Program

Documentation for the MODFLOW 6

Groundwater Transport Model
Chapter 61 of
Section A, Groundwater
Book 6, Modeling Techniques

Techniques and Methods 6–A61

U.S. Department of the Interior

U.S. Geological Survey
Cover. Binary computer code illustration.
Documentation for the MODFLOW 6 Groundwater
Transport Model
By Christian D. Langevin, Alden M. Provost, Sorab Panday, and Joseph D. Hughes

Chapter 61 of
Section A, Groundwater
Book 6, Modeling Techniques

Techniques and Methods 6–A61

U.S. Department of the Interior

U.S. Geological Survey
U.S. Geological Survey, Reston, Virginia: 2022

For more information on the USGS—the Federal source for science about the Earth, its natural and living resources,
natural hazards, and the environment—visit [Link] or call 1–888–ASK–USGS.

For an overview of USGS information products, including maps, imagery, and publications, visit
[Link]

Any use of trade, frm, or product names is for descriptive purposes only and does not imply endorsement by the U.S.
Government.
Although this information product, for the most part, is in the public domain, it also may contain copyrighted materials
as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner.

Suggested citation:
Langevin, C.D., Provost, A.M., Panday, Sorab, and Hughes, J.D., 2022, Documentation for the MODFLOW 6
Groundwater Transport Model: U.S. Geological Survey Techniques and Methods, book 6, chap. A61, 56 p.,
[Link]

ISSN 2328-7055 (online)

 iii

Preface
The report describes the Groundwater Transport (GWT) Model for the U.S. Geological Survey (USGS) modular
hydrologic simulation program called MODFLOW 6. The program can be be downloaded from the USGS for
free. The performance and accuracy of the GWT Model has been tested in a variety of applications. Future
applications, however, might reveal errors that were not detected in the test simulations. Users are requested
to send notifcation of any errors found in this model documentation report or in the model program to the
MODFLOW contact listed on the Web page. Updates might be made to both the report and to the model pro-
gram. Users can check for updates on the MODFLOW Web page ([Link]

Acknowledgments
The authors are grateful for the constructive reviews provided by Vivek Bedekar of S.S. Papadopulos & Asso-
ciates, Inc. and Eric D. Morway of the U.S. Geological Survey. Editorial review of this report was conducted by
Angel Martin. The authors are also grateful for technical input provided by U.S. Geological Survey colleagues.
iv

Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Chapter 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–1
History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–2
Overview of the Groundwater Transport Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–3
Information for Existing Solute Transport Modelers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–6
Organization and Scope of This Report . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–8
Chapter 2. Formulation and Solution of the Control-Volume Finite-Difference Equation . . . . . . . . . . . . . . . . . . . . 2–1
Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–1
Control-Volume Finite-Difference Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–2
Characteristics of a Model Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–3
Control-Volume Finite-Difference Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–3
Numerical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–5
Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–6
Flow Model Interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–6
Flow Imbalance Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–6
Time Stepping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–7
Special Considerations for Dry Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–7
Chapter 3. Mobile Storage and Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–1
Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–1
Sorption . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–1
Linear Isotherm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–3
Freundlich Isotherm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–3
Langmuir Isotherm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–4
Considerations for Unconfned Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–4
Decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–5
First Order Decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–5
Zero-Order Decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–5
Decay of Sorbed Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–6
Linear Isotherm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–6
Nonlinear Freundlich and Langmuir Isotherms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3–7
Chapter 4. Advective and Dispersive Solute Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–1
Advection (ADV) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–1
Central-In-Space Weighting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–1
Upstream Weighting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–2
Total Variation Diminishing (TVD) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–2
Numerical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–3
 v

Dispersion (DSP) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–4

Mechanical Dispersion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–6
Simplifed Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–9
XT3D Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–12
Chapter 5. Sources and Sinks of Solute Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5–1
Constant-Concentration (CNC) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5–1
Flow Model Source and Sink Mixing (SSM) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5–1
Source Mass Loading (SRC) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5–3
Chapter 6. Transport for Advanced Stress Packages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–1
Storage Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–2
Advection Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–2
Water Mover Terms and the Mover Transport (MVT) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–3
Streamfow Transport (SFT) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–3
Lake Transport (LKT) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–4
Multi-Aquifer Well Transport (MWT) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–5
Unsaturated Zone Transport (UZT) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–5
Chapter 7. Immobile Domain Storage and Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7–1
References Cited . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . R–1

Figures
1–1. Diagram showing the structure of a MODFLOW 6 simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–4
1–2. Diagram showing the domains and packages for the Groundwater Transport Model. . . . . . . . . . . . 1–5
2–1. Diagram of a porous media model cell that is partially saturated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–4
4–1. Diagram showing, in two dimensions, the cell connections used by the XT3D method to
estimate the concentration gradient. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4–14
7–1. Diagram showing two different conceptualizations for mobile and immobile domains. . . . . . . . . . . 7–2

Tables
1–1. List of packages available for use with the Groundwater Transport Model . . . . . . . . . . . . . . . . . . . . . 1–6
5–1. List of Groundwater Flow Model Packages that can act as a solute source or sink for the
MODFLOW 6 Groundwater Transport Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5–2
6–1. Source and sink equations for Streamfow Transport Package of the MODFLOW 6 Ground-
water Transport Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–4
6–2. Source and sink equations for Lake Transport Package of the MODFLOW 6 Groundwater
Transport Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–5
6–3. Source and sink equations for Multi-Aquifer Well Transport Package of the MODFLOW 6
Groundwater Transport Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–5
vi

6–4. Source and sink equations for Unsaturated Zone Transport Package of the MODFLOW 6
Groundwater Transport Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6–6
Documentation for the MODFLOW 6 Groundwater
Transport Model
By Christian D. Langevin,1 Alden M. Provost,1 Sorab Panday,2 and Joseph D. Hughes1

Abstract
This report documents a new Groundwater Transport (GWT) Model for MODFLOW 6. The GWT Model
simulates three-dimensional transport of a single chemical species in fowing groundwater based on a gen-
eralized control-volume fnite-difference approach. Although each GWT Model is only able to represent a
single chemical species, multiple GWT Models may be invoked within a single MODFLOW 6 simulation to
represent solute transport of multiple non-interacting chemical species. The GWT Model is designed to work
with the Groundwater Flow (GWF) Model for MODFLOW 6, which simulates transient, three-dimensional
groundwater fow. The version of the GWT model documented here must use the same spatial discretization
used by the GWF Model; however, that spatial discretization can be represented by regular MODFLOW grids
consisting of layers, rows, and columns, or by more general unstructured grids. The GWT Model simulates
(1) advective transport, (2) the combined hydrodynamic dispersion processes of velocity-dependent mechan-
ical dispersion and molecular diffusion, (3) adsorption and absorption (collectively referred to as sorption) of
solutes by the aquifer matrix, (4) transfer between the mobile domain and one or more immobile domains, (5)
frst- or zero-order solute decay or production, (6) mixing from groundwater sources and sinks, and (7) direct
addition of solute mass. The GWT Model can also represent advective solute transport through advanced
package features, such as streams, lakes, multi-aquifer wells, and the unsaturated zone. If the GWF Model
application uses the Water Mover (MVR) Package to connect fow packages, then solute transport between
these packages can also be represented. The transport processes described in this report have been imple-
mented in a fully implicit manner and are solved in a system of equations using iterative numerical methods.
The present version of the GWT Model for MODFLOW 6 does not have an option to calculate steady-state
transport solutions; if a steady-state solution is required, then transient evolution of the solute must be repre-
sented using multiple time steps until no further changes in solute concentrations are detected.

1 U.S. Geological Survey.

2 GSI Environmental Inc.
 Chapter 1. Introduction 1–1

Chapter 1. Introduction
MODFLOW 6 is the latest core version of the MODFLOW software. It was released by the U.S. Geo-
logical Survey (USGS) in 2017 (Hughes and others, 2017). This new version of MODFLOW was redesigned
from scratch using an object-oriented design that allows for multiple models to be included in a single simu-
lation. The new MODFLOW 6 framework facilitates more than one model of the same type in a simulation.
For example, there may be a regional groundwater fow model and a locally refned inset model, and these
models can be tightly coupled along their interface. A unique feature of the MODFLOW 6 framework allows
these models to be coupled at the matrix solution level whereby a single system of equations is constructed
and solved simultaneously for both models.
In addition to supporting multiple models of the same type in a single simulation, the MODFLOW 6
framework was also designed to support multiple models of different types in the same simulation. The frst
type of model introduced in MODFLOW 6 was the Groundwater Flow (GWF) Model (Langevin and others,
2017), which simulates three-dimensional, transient, groundwater fow. The GWF Model synthesizes many of
the newer capabilities that were added to the MODFLOW-2005 program (Harbaugh, 2005). The following list
summarizes select features available in the GWF Model for MODFLOW 6:
‚ The GWF Model includes a Newton fow formulation for handling cell wetting and drying that can occur
in unconfned aquifers (Niswonger and others, 2011; Panday and others, 2013).
‚ The GWF Model can represent groundwater fow using a regular MODFLOW grid consisting of layers,
rows, and columns, but the GWF Model also supports unstructured grids, following the approach imple-
mented in Panday and others (2013).
‚ The GWF Model has a set of advanced stress packages for representing streams, lakes, multi-aquifer
wells, and fow in the unsaturated zone.
‚ The GWF Model includes a rule-based approach to transfer water between the different stress packages.
The approach is implemented in the Water Mover (MVR) Package (Morway and others, 2021) to transfer
water from providers to receivers.
‚ Because of its object-oriented design a single GWF Model can contain multiple packages of the same
type. This feature allows separate Well Package input fles to be created for each well feld, for example.
This report describes and documents a new Groundwater Transport (GWT) Model for MODFLOW 6. The
GWT Model simulates three-dimensional transport of a single solute species in fowing groundwater. Simula-
tion of changing solute concentrations requires the solution of a partial differential equation governing solute
transport. The GWT Model solves the solute transport equation using numerical methods and a generalized
control-volume fnite-difference approach, which can be used with regular MODFLOW grids or with unstruc-
tured grids. The GWT Model is designed to work with most of the new capabilities released with the GWF
Model, including the Newton fow formulation, unstructured grids, advanced packages, and the movement of
water between packages. The GWF and GWT Models operate simultaneously during a MODFLOW 6 simu-
lation to represent coupled groundwater fow and solute transport. The GWT Model can also run separately
from a GWF Model by reading the heads (groundwater levels) and fows saved by a previously run GWF
Model. The GWT model is also capable of working with the fows from another groundwater fow model,
provided the fows from that model are written in the correct form to fow and head fles.
The purpose of the GWT Model is to calculate changes in solute concentration in both space and time.
Solute concentrations within an aquifer can change in response to multiple solute transport processes. These
processes include (1) advective transport of solute with fowing groundwater, (2) the combined hydrody-
namic dispersion processes of velocity-dependent mechanical dispersion and molecular diffusion, (3) sorp-
tion of solutes by the aquifer matrix either by adsorption to individual solid grains or by absorption into solid
1–2 Documentation for the MODFLOW 6 Groundwater Transport Model

grains, (4) transfer of solute into low permeability aquifer material (called an immobile domain) where it can
be stored and later released, (5) frst- or zero-order solute decay or production in response to chemical or bio-
logical reactions, (6) mixing with fuids from groundwater sources and sinks, and (7) direct addition of solute
mass.
With the present GWT Model implementation, there can be multiple domains and multiple phases. There
is a single mobile domain, which normally consists of fowing groundwater, and there can be one or more
immobile domains. The GWT Model simulates the dissolved phase of chemical constituents in both the
mobile and immobile domains. The dissolved phase is also referred to in this report as the aqueous phase. If
sorption is represented, then the GWT Model also simulates the solid phase of the chemical constituent in both
the mobile and immobile domains. The dissolved and solid phases of the chemical constituent are tracked in
the different domains by the GWT Model and can be reported as output as requested by the user. There is no
provision in the version of the GWT Model described here to calculate and track a vapor phase.

History
Numerical models are often used to simulate and predict the fate and transport of a dissolved chemical
constituent in groundwater. These models range in complexity from simple particle tracking models to disper-
sive solute transport models. This section summarizes some of the popular particle tracking and solute trans-
port models that have been developed for MODFLOW.
A frst step in simulating and predicting the fate and transport of a dissolved chemical constituent is
development and calibration of a groundwater fow model, such as MODFLOW, for the area of interest. If
the objective of a transport investigation is merely to develop an understanding of groundwater fow direc-
tions and rates, then simple particle tracking methods can be used to simulate and predict predict groundwa-
ter travel times and paths. These groundwater travel times and paths are often a good frst approximation of
solute movement if the solute is conservative (it does not sorb or react with the aquifer) and transport is advec-
tion dominated such that hydrodynamic dispersion can be neglected. To this end, MODPATH offers a semi-
analytical particle tracking method designed to work with simulated fows from a MODFLOW groundwater
fow simulation. A recent version of MODPATH (Version 7), was designed to work with the simulated fows
from MODFLOW 6, provided the model grid used for the fow simulation consisted of rectangular cells (Pol-
lock, 2016) . Unstructured model grids consisting of quad-based refnement, such as quadpatch and quadtree
grids, meet this criterion, and can be used with MODPATH Version 7. MODPATH Version 7 has several limi-
tations that restrict its use with the GWF Model, including (1) simulated particle paths are not consistent with
the Newton fow formulation if cells are dry, (2) the model grid used for fow must have rectangular cells, and
(3) particles cannot be tracked through advanced model fow packages.
The USGS distributes a MODFLOW-based solute transport model, called the Groundwater Transport
(GWT) Process. The GWT Process (not to be confused with the GWT Model for MODFLOW 6 described
in this report) is implemented directly in MODFLOW-2000, and is referred to as MODFLOW-GWT.
MODFLOW-GWT is based on the MOC2D (Konikow and Bredehoeft, 1978) and MOC3D (Konikow and
others, 1996) programs and simulates advection, hydrodynamic dispersion, and simple reactions for a single
chemical species. Winston and others (2018) document a volume-weighted particle tracking method for the
most recent release of this solute transport model. MODFLOW-GWT has many advanced capabilities, includ-
ing solute routing through lakes (Merritt and Konikow, 2000), streams (Prudic and others, 2004), and multi-
node wells (Hornberger and Konikow, 2006). These advanced capabilities provided the impetus for many of
the design features described in this report. MODFLOW-GWT does not work with some of the newer capabil-
ities developed for MODFLOW, such as the Newton fow formulation and unstructured grids.
MT3D is another popular numerical model developed to simulate solute transport using fows from a
MODFLOW simulation. MT3D is not implemented directly inside of MODFLOW. Rather, it reads heads and
fows from a separate fle that MODFLOW creates while the model is running. The fow and transport link fle
 Chapter 1. Introduction 1–3

has been updated over time to work with MODFLOW-2000 (Zheng and others, 2001) and newer MODFLOW
versions. The original version of MT3D (Zheng, 1990) was a single species solute transport model, but MT3D
was later extended to represent dual-domain transport and multiple chemical species under the MT3DMS
name (Zheng and Wang, 1999) (where “MS” was added to indicate multiple species). Additional capabili-
ties were added to MT3DMS to support more MODFLOW packages, simulate zero-order reactions, and solve
for steady-state conditions, culminating in MT3DMS Version 5.3 (Zheng, 2010). MT3DMS has been widely
used by consultants, academics, and government agencies to simulate complex groundwater transport prob-
lems. MT3DMS is used in other codes, such as SEAWAT, to solve the transport equation for variable-density
fow problems (Guo and Langevin, 2002; Langevin and others, 2003, 2008), and RT3D and PHT3D to simu-
late reactive transport (Clement, 1997; Prommer and others, 2003). MT3DMS does not work with some of the
newer capabilities developed for MODFLOW, such as the Newton fow formulation and unstructured grids.
Bedekar and others (2016) extended MT3DMS to work with the Newton fow formulation and transport in
streams, lakes, and the unsaturated zone. Many other changes were also made as part of this extension, which
resulted in a new version of the program called MT3D-USGS. At this time (2022), MT3D-USGS is continu-
ing to receive updates and fxes in response to user requests. MT3D-USGS can be used with heads and fows
from a MODFLOW 6 simulation with a GWF Model provided a regular MODFLOW grid is used for the fow
simulation. MT3D-USGS cannot be used to simulate transport within the fow domains represented by the
advanced stress packages of the GWF Model.
MODFLOW-USG is a popular unstructured grid version of MODFLOW (Panday and others, 2013) that
can be used with regular MODFLOW grids or fexible unstructured grids. MODFLOW-USG has a robust
Newton fow formulation for solving diffcult unconfned aquifer problems, and it has the Connected Lin-
ear Network (CLN) Process for simulating fow in boreholes, fractures, and conduits. Since its release by the
USGS, MODFLOW-USG has been updated to include many new capabilities, including the addition of a gen-
eralized control-volume solute transport model (Panday, 2020). The solute transport model in MODFLOW-
USG runs concurrently with the fow model to simulate the fate and transport of multiple chemical species in
response to advection, hydrodynamic dispersion, sorption, and zero- or frst-order growth and decay. Many
of the features and concepts developed for MODFLOW-USG have provided the foundation for the GWF and
GWT Models in MODFLOW 6.

Overview of the Groundwater Transport Model

The GWT Model documented here builds on the successes of previous solute transport models and takes
advantage of design and solver advances implemented in the MODFLOW 6 framework. Following the mod-
ular concepts developed by earlier MODFLOW developers, the GWT model is highly extensible so that new
capabilities and packages can be added in the future. Because the GWT Model is part of the broader MOD-
FLOW 6 framework, it conforms to a consistent design that should make it look, feel, and operate like the
GWF Model and other models being developed for the framework. As the MODFLOW 6 framework evolves
and new enhancements are made, these enhancements will also apply to the GWT Model.
The GWT Model works with the heads and fows generated by the GWF Model. These simulated heads
and fows can be passed through memory from the GWF Model to a GWT Model running concurrently as part
of a single MODFLOW 6 simulation. Alternatively, a GWF Model can be run frst to save simulated heads
and fows to binary fles. Then, a second MODFLOW 6 simulation with a GWT Model can be used to repre-
sent solute transport, using the saved heads and fows from the frst simulation with the GWF Model. There
are advantages and disadvantages to both approaches.
For the single-simulation use case, a GWT Model will run within a MODFLOW 6 simulation in conjunc-
tion with a GWF Model (see fgure 1–1a). The GWF Model will calculate heads and fows, which are made
available to the GWT Model through the GWF-GWT Exchange. Time stepping is controlled by the Tempo-
ral Discretization (TDIS) Package. Heads and fows are calculated for the GWF Model by the Iterative Model
1–4 Documentation for the MODFLOW 6 Groundwater Transport Model

(a) MODFLOW 6
Simulation

Temporal Groundwater Groundwater Flow-Transport Iterative Model Iterative Model

Discretization Flow (GWF) Transport Exchange Solution (IMS) Solution (IMS)
(TDIS) Model (GWT) Model (GWF-GWT) for Flow for Transport

(b)
MODFLOW 6 MODFLOW 6
Simulation Simulation

Temporal Groundwater Temporal Groundwater

Iterative Model Iterative Model
Discretization Flow (GWF) Discretization Transport
Solution (IMS) Solution (IMS)
(TDIS) Model (TDIS) (GWT) Model

Figure 1–1. Structure of a MODFLOW 6 simulation. (a) Flow and transport models are part of the same simulation. (b) Flow and
transport models are in separate simulations. In this case with separate simulations, GWF Model fows are saved to a binary
fle, which is read as input to the transport model.

Solution (IMS). A second IMS is used to solve for concentration and solute fuxes of the GWT Model. Alter-
natively, the GWF and GWT Models can be run as separate simulations (fgure 1–1b). In this case, the user
runs the GWF Model frst and saves all heads and fows for every time step to binary fles. The user then runs
the GWT Model, which reads the heads and fows as input.
The GWT Model described in this report is divided into “packages.” A package is a part of the model
that deals with a single aspect of simulation. For example, the Advection Package simulates the transport pro-
cess of advection, and the Dispersion Package simulates the transport processes of mechanical dispersion and
molecular diffusion. Some packages are always required for a simulation, whereas other packages are only
activated if their capabilities are needed for a particular application.
The GWT Model is comprised of the packages shown in fgure 1–2. The packages shown on the left of
fgure 1–2 are used to provide data to the model, such as discretization information, initial concentrations, the
frequency and type of output to save, locations and types of observations to save, and information on how the
GWT Model interfaces with the GWF Model. These data input packages do not represent transport processes,
but are needed to provide information for the GWT Model. The remaining packages shown in fgure 1–2 are
separated into mobile domain and immobile domain packages. The mobile domain represents the “fast” part
of the groundwater system in which a dissolved constituent is transported through an aquifer with fowing
groundwater. The immobile domain represents the “slow” part of the system in which groundwater movement
can be considered negligible. Solute mass can move between the mobile and immobile domains in response to
a transfer coeffcient and the concentration difference between the mobile and immobile domains. A unique
aspect of the GWT Model described here is that there can be any number of immobile domains, each with
their own transfer coeffcients, immobile domain porosity values, sorption and decay parameters, and calcu-
lated concentrations. The mobile domain has packages for representing solute sources and sinks, the direct
addition of solute mass, and the specifcation of constant concentration conditions. Both the mobile and immo-
bile domains have packages for entering properties for internal transport (such as porosity, and sorption and
decay parameters). Specifcation of one or more immobile domains is optional; however, most simulations
will require specifcation of mobile domain packages.
The various packages of the GWT Model documented in this report, the shortened character abbreviation
used for each package, and the package category are listed in table 1–1.
 Chapter 1. Introduction 1–5

Groundwater
Transport
(GWT) Model

Immobile
Mobile Domain
Domain

Discretization Source and Mobile Storage Streamflow Immobile Immobile Immobile

Advection Sink Mixing
(DIS, DISV, and Transfer Transport Storage and Storage and Storage and
(ADV) (SSM)
DISU) (MST) (SFT) Transfer (IST) Transfer (IST) Transfer (IST)

Initial Dispersion Mass Source Lake Transport

Conditions (IC) (DSP) Loading (SRC) (LKT)

Constant Multi-Aquifer
Output Control Well Transport
Concentration
(OC) (MWT)
(CNC)

Unsaturated
Observations Zone Transport
(OBS) (UZT)

Mover
Flow Model Transport
Interface (FMI) (MVT)

Data Input Cell-to-Cell External Within-Cell Advanced Within-Cell Within-Cell Within-Cell

Process Stress Process Package Process Process Process

Figure 1–2. Domains and packages for the MODFLOW 6 Groundwater Transport (GWT) Model. There are three Immobile Stor-
age and Transfer (IST) Packages shown here, however, the user can include as many IST Packages as needed for a model sim-
ulation.
1–6 Documentation for the MODFLOW 6 Groundwater Transport Model

Table 1–1. List of packages available for use with the Groundwater Transport Model.

Package Name Abbreviation Package Category

Spatial Discretization DIS, DISV, or DISU Data Input
Initial Conditions IC Data Input
Model Observations OBS Data Input
Output Control OC Data Input
Flow Model Interface FMI Data Input
Advection ADV Cell-to-Cell Process
Dispersion DSP Cell-to-Cell Process
Source-Sink Mixing SSM External Stress
Mass Source Loading SRC External Stress
Constant Concentration CNC External Stress
Mobile Storage and Transfer MST Within-Cell Process
Immobile Storage and Transfer IST Within-Cell Process
Streamfow Transport SFT Advanced Package
Lake Transport LKT Advanced Package
Multi-Aquifer Well Transport MWT Advanced Package
Unsaturated Zone Transport UZT Advanced Package
Mover Transport MVT Advanced Package

Information for Existing Solute Transport Modelers

The MODFLOW 6 GWT Model contains most of the functionality of MODFLOW-GWT, MT3DMS,
MT3D-USGS and MODFLOW-USG. The following list summarizes major differences between the GWT
Model in MODFLOW 6 and previous MODFLOW-based solute transport programs.
1. The GWT Model simulates transport of a single chemical species; however, because MODFLOW 6
allows for multiple models of the same type to be included in a single simulation, multiple species can
be represented by using multiple GWT Models. A clear advantage of this design is that each chemical
species has a dedicated GWT Model. The fexibility allows all transport input parameters to be specifed
separately for each species, including initial conditions, boundary conditions, output frequency, and so
forth. A disadvantage of this approach is that the present implementation requires separate memory for
each GWT Model, which may not be effcient.
2. There is no specialized fow and transport link fle (Zheng and others, 2001) used to pass the simulated
groundwater fows to the transport model, as is required by MT3D. Instead, simulated fows from the
GWF Model are passed in memory to the GWT Model while the program is running. Alternatively, the
GWT Model can read the standard binary head and fow fles saved by the GWF Model during a pre-
ceding fow simulation. If the user intends to simulate transport through the advanced stress packages
and Water Mover Package, then fows from these advanced packages must also be saved to binary fles.
Names for these binary fles are provided as input to the Flow Model Interface (FMI) Package.
3. The GWT Model is based on a generalized control-volume fnite-difference method, which means that
solute transport can be simulated using regular MODFLOW grids (specifed with the DIS Package),
which consist of layers, rows, and columns, or unstructured MODFLOW grids (specifed with the DISV
or DISU Packages; table 1–1).
 Chapter 1. Introduction 1–7

4. Advection can be simulated using central-in-space weighting, upstream weighting, or an implicit second-
order Total Variation Diminishing (TVD) scheme. The GWT Model does not have the Method of Char-
acteristics (particle-based approaches) or an explicit TVD scheme. Consequently, the GWT Model may
require a higher level of spatial discretization than other transport models that use higher order terms for
advection dominated groundwater systems. This can be an important limitation for some problems, which
require the preservation of sharp solute fronts.

5. Variable-density fow and transport can be simulated by including a GWF Model and a GWT Model in
the same MODFLOW 6 simulation. The Buoyancy Package should be activated for the GWF Model if
density variations are expected to affect groundwater fow so that fuid density is calculated as a function
of simulated concentration. If more than one chemical species is represented then the Buoyancy Package
allows the simulated concentration for each species to be used in the density equation of state. Langevin
and others (2020) describe the hydraulic-head formation that is implemented in the Buoyancy Package
for variable-density groundwater fow and present results from MODFLOW 6 variable-density simu-
lations. The variable-density capabilities available in MODFLOW 6 replicate and extend the capabili-
ties available in SEAWAT to include, for example, the Newton fow formulation for unconfned aquifers,
transport through advanced packages, and transport using unstructured grids.

6. The GWT model includes the MST and IST Packages (fg. 1–2). These two packages collectively com-
prise the capabilities of the MT3DMS Reaction Package.

7. The MST Package supports the linear isotherm for representation of sorption as well as the nonlinear
Freundlich and Langmuir isotherms. The MST Package described in this report does not support the
nonequilibrium sorption model that is available in MT3DMS. The IST Package supports only linear sorp-
tion.

8. The GWT Model was designed so that the user can specify as many immobile domains as necessary
to represent observed contaminant transport patterns and solute breakthrough curves. The effects of an
immobile domain are represented using the IST Package, and the user can specify as many IST Packages
as necessary.

9. The GWT Model documented in this report does not support kinetic reactions between species, a feature
that is available in MT3D-USGS. Interactions between species in separate GWT Models is an option that
may be developed in the future through a GWT-GWT Exchange.

10. There is no option to automatically run the GWT Model to steady state using a single time step. This is
an option available in MT3DMS (Zheng, 2010). Steady-state conditions must be determined by running
the transport model under transient conditions until solute concentrations stabilize.

11. The GWT Model described in this report is capable of simulating solute transport in the advanced stress
packages of MODFLOW 6, including the Lake, Streamfow Routing, Multi-Aquifer Well and Unsatu-
rated Zone Transport Packages (Langevin and others, 2017). Solute transport between these advanced
packages is also supported, such as the transport of solute from a stream into a lake. The present imple-
mentation simulates solute advection between package features, such as between two stream reaches, but
dispersive transport between package features is not represented. Similarly, solute transport between the
advanced packages and the aquifer is only through advection.

12. There are many other differences between the MODFLOW 6 GWT Model and other solute transport
models that work with MODFLOW, especially with regards to program design and input and output.
1–8 Documentation for the MODFLOW 6 Groundwater Transport Model

MODFLOW 6 input and output are described in a separate user guide, which is included with the distri-
bution. A full suite of test and example problems for the GWT Model is also included with the software
distribution available on the internet for download.

Organization and Scope of This Report

This report documents the GWT Model for MODFLOW 6 and contains the underlying mathematical
equations and fundamental concepts. The report focuses on foundational equations and concepts that are
expected to persist even as new versions of the model are released. Equations and concepts that may change
as the program evolves are mentioned only briefy here, and users should consult the release notes and supple-
mentary information that is distributed with MODFLOW 6. Likewise, instructions for preparing input fles,
running the model, and post-processing the results are not included in this report. Instead, user input instruc-
tions are provided as a separate document with the program distribution. Distributing the input instructions
separately allows them to evolve and change with the program as new options are added. In some instances,
users may need to refer to additional reports that describe new capabilities added to MODFLOW 6, including
capabilities to simulate additional solute transport processes. The MODFLOW 6 distribution also includes a
rich set of example problems for the GWT Model. These example problems contain problem descriptions and
model input fles so that users can run the examples and evaluate the results.
This report is organized as follows. Chapter 2 describes the derivation of the control-volume fnite-
difference (CVFD) equation upon which the GWT Model is based. Chapter 3 describes the Mobile Storage
and Transfer (MST) Package, which calculates changes in solute storage, sorption, and solute growth and
decay processes. Chapter 4 describes the equations and approaches for simulating the transport of solute due
to advection and dispersion. Chapter 5 describes the representation of simple solute sources and sinks, which
can be represented with the Constant Concentration (CNC) Package, the Source and Sink Mixing (SSM) Pack-
age for sources and sinks associated with the groundwater fow model, and the mass source loading (SRC)
Package, which represents the direct addition of solute mass to a model cell. Chapter 6 describes the advanced
package transport capabilities, which can be used to simulate solute transport in streams, lakes, multi-aquifer
wells, and in the unsaturated zone. Chapter 7 describes the Immobile Storage and Transfer (IST) Package,
which can be used to represent one or more immobile domains embedded within the aquifer. These domains
may correspond to silt or clay lenses, for example, that can store or release solute mass from or to the mobile
domain. The GWT Model packages documented in this report are those listed in table 1–1.
 Chapter 2. Formulation and Solution of the Control-Volume Finite-Difference Equation 2–1

Chapter 2. Formulation and Solution of the Control-Volume Finite-

Difference Equation
The groundwater transport equation in MODFLOW 6 is discretized using a control-volume fnite-
difference (CVFD) method. This chapter describes the mathematical equations discretized in the GWT model,
discretization options, and the general forms of the fnite-difference equations used to simulate solute transport
within an aquifer.

Mathematical Model
Transport of a solute dissolved in groundwater is described mathematically by a partial differential equa-
tion that represents the conservation of solute mass (Konikow and Grove, 1977; Zheng and Bennett, 2002). At
any location, the accumulation of solute mass is equal to the difference between mass entering and mass leav-
ing a specifed volume of aquifer. Such an equation can be written in a variety of different forms, including the
following form (eq. 2–1), which includes the transport mechanisms represented by the GWT Model:

B pSw θCq
“ ´∇ ¨ pqCq ` ∇ ¨ pSw θD∇Cq ` qs1 Cs ` Ms ´ λ1 θSw C ´ γ1 θSw
Bt
` ˘ nim (2–1)
B Sw C ÿ
´fm ρb ´ λ2 fm ρb Sw C ´ γ2 fm ρb Sw ´ ζim Sw pC ´ Cim q ,
Bt im“1

where Sw is the water saturation (dimensionless) defned as the volume of water per volume of voids, θ is the
effective porosity of the mobile domain (dimensionless), defned as volume of voids participating in mobile
transport per unit volume of aquifer, C is volumetric concentration of the mobile domain expressed as mass of
dissolved solute per unit volume of fuid (M {L3 ), t is time (T ), q is the vector of specifc discharge (L{T ), D
is the second-order tensor of hydrodynamic dispersion coeffcients (L2 {T ), qs1 is the volumetric fow rate per
unit volume of aquifer (defned as positive for fow into the aquifer) for mass sources and sinks (1{T ), Cs is
the volumetric solute concentration of the source or sink fuid (M {L3 ), Ms is rate of solute mass loading per
unit volume of aquifer (M {L3 T ), λ1 is the frst-order decay rate coeffcient for the liquid phase (1{T ), γ1 is
the zero-order decay rate coeffcient for the liquid phase (M {L3 T ), fm is the fraction of aquifer solid material
available for sorptive exchange with the mobile phase under fully saturated conditions, ρb is the bulk density
of the aquifer material (M {L3 ), C is the sorbed concentration of solute mass in the mobile domain (M {M ),
λ2 is the frst-order decay rate coeffcient (1{T ) for the sorbed phase of the mobile domain, γ2 is the zero-order
decay rate coeffcient for the sorbed phase of the mobile domain (M {M T ), nim is the number of immobile
domains, ζim is the rate coeffcient for the transfer of mass between the mobile domain and immobile domain
im (1{T ), and Cim is the solute concentration for immobile domain im (M {L3 ).
Equation 2–1 can be rewritten in the following manner to correspond to the design of the GWT Model (see
table 1–1 for package names and abbreviations).

` ˘
B pSw θCq B Sw C
´ ´ f m ρb ´ λ1 θSw C ´ γ1 θSw ´ λ2 fm ρb Sw C ´ γ2 fm ρb Sw
Bt Bt
loooooooooooooooooooooooooooooooooooooooooooooooooooomoooooooooooooooooooooooooooooooooooooooooooooooooooon
M ST
nim
ÿ (2–2)
1
´ ∇ ¨ pqCq `∇ ¨ pS w θD∇Cq
looooomooooon looooooooomooooooooon `q C `M
s s loomoosn
loomoon ´ ζim Sw pC ´ Cim q “ 0.
im“1
SRC loooooooooooooomoooooooooooooon
ADV DSP SSM
IST
2–2 Documentation for the MODFLOW 6 Groundwater Transport Model

In this form (eq. 2–2), the terms have been grouped and labeled according to the GWT transport package that
represents their effects. The effects of solute storage, sorption, and decay on the mobile domain are repre-
sented with the Mobile Storage and Transfer (MST) Package. The effects of advection are represented with the
Advection (ADV) Package. The effects of hydrodynamic dispersion, which includes mechanical dispersion
and molecular diffusion, are represented by the Dispersion (DSP) Package. Groundwater infow and outfow
from stress packages in the fow model are represented with the Source and Sink Mixing (SSM) Package. The
direct loading of solute mass can be represented with the Mass Source Loading (SRC) Package (table 1–1).
The effects of a diffusive exchange between the mobile domain and an immobile domain is represented with
the Immobile Storage and Transfer (IST) Package or multiple instances of the IST Package if the intent is to
represent multiple immobile domains.
Separating the transport equation in the manner described above is a common way of breaking the prob-
lem into individual terms that can be addressed independently of one another. Zheng (1990) and Zheng and
Wang (1999) were among the frst to write the transport equation in this manner, which allowed them to write
“packages” to solve the individual terms. Their approach built on the MODFLOW concept of a package,
which makes it relatively easy for a user to learn about individual processes, turn packages on and off, and
develop new packages. A simplifed form of equation 2–2 is

f M ST ` f ADV ` f DSP ` f SSM ` f SRC ` f IST “ 0, (2–3)

where

` ˘
M ST B pSw θCq B Sw C
f “´ ´ fm ρb ´ λ1 θSw C ´ γ1 θSw ´ λ2 fm ρb Sw C ´ γ2 fm ρb Sw
Bt Bt
f ADV “ ´∇ ¨ pqCq
f DSP “ ∇ ¨ pSw θD∇Cq
f SSM “ qs1 Cs (2–4)

f SRC “ Ms
nim
ÿ
f IST “ ´ ζim Sw pC ´ Cim q .
im“1

Control-Volume Finite-Difference Method

Equation 2–1, together with initial conditions and any relevant boundary conditions, represents mathemat-
ically the solute mass balance at any point in the model domain. In certain simple cases, equation 2–1 can be
solved analytically to obtain a mathematical expression for the distribution of solute concentration through-
out a model domain. For models of real-world feld sites, which tend to be too complex to solve analytically,
numerical solutions are often sought. Using the CVFD method, MODFLOW 6 discretizes the model domain
into cells. The balance of solute mass is formulated for each model cell, taking into account the fows of
solute to and from neighboring cells by advection and dispersion, as well as external sources and sinks. Taken
together, the solute mass balance equations for all the cells form a system of linear equations that is solved
iteratively using a linear matrix solver. Details of the CVFD implementation in the GWT Model are described
below.
 Chapter 2. Formulation and Solution of the Control-Volume Finite-Difference Equation 2–3

Characteristics of a Model Cell

As described by Langevin and others (2017) a MODFLOW 6 groundwater model cell is a prism with ver-
tical sides and a horizontal top and bottom. An example of a groundwater model cell, with the size of individ-
ual sediment grains exaggerated to illustrate porosity concepts, is shown in fgure 2–1. The area of cell n in
plan view is defned as An . The top and bottom cell elevations are defned as T OPn and BOTn , respectively,
and the head in the cell is defned as hn . The total volume of the cell, including the solid and void space, is
defned as

Vcell “ An pT OPn ´ BOTn q . (2–5)

Equation 2–5 is used to calculate cell volume in MODFLOW 6 instead of one based on cell widths and depths,
because although cells are prisms, they can have shapes other than rectangles in plan view. The volume of
void space (Vvoid ) and the volume of space occupied by solids (Vsolid ) within a cell is related to porosity by

Vvoid Vvoid
θ“ “ . (2–6)
Vcell pVvoid ` Vsolid q

In this context, porosity is intended to represent the pores available for mobile transport, also referred to as the
effective porosity. This effective porosity is different from the immobile domain porosity, which is defned
in chapter 7 “Immobile Domain Storage and Transfer”. Unless stated otherwise, porosity refers to effective
porosity.
The volume of water is related to the cell saturation, Sw , which is calculated by the fow model. When a
cell is fully saturated, it has a saturation of one, which means all of the pore space is flled with water. When a
cell is completely dry, it has a saturation of zero. When a cell is partially saturated, the saturation is assumed to
depend on the head value within the cell as

h ´ BOT
Sw “ . (2–7)
pT OP ´ BOT q

If the Newton-Raphson formulation is used in the fow model, then the equation for cell saturation is slightly
different in that the transitions from fully saturated to partially saturated and from partially saturated to dry are
smoothed (Langevin and others, 2017).
From these equations, the volume of water in a cell, Vw , is expressed as

Vw “ θVcell Sw . (2–8)

The volume of water in a cell may change during the simulation in response to infows and outfows, as calcu-
lated by the fow model. An inherent assumption in the approximation of the volume of water in a cell using
this approach is that the aquifer porosity does not change with time. As noted by Goode (1990), this is not
strictly correct, because changes in head and the corresponding change in the volume of water in storage is
normally associated with a change in porosity. The GWT Model, however, does not explicitly account for the
change in aquifer porosity that results from water being added to or released from storage.

Control-Volume Finite-Difference Equation

A general solute mass balance equation for a model cell n can be written as
 2–4

M9 nADV “

mulations for Fn,m

ÿ

mPηn
considered to be unsaturated.

Fn,m
ADV
,

ADV are discussed in chapter 4.

Documentation for the MODFLOW 6 Groundwater Transport Model

groundwater storage and the sorbed phase and by decay or production.

where ηn is the set of neighbors of (cells connected to) cell n, and Fn,m
advective fows of solute mass between cell n and each of its neighbors:
sha1_base64="HlbaWngJl2+9LR1NU+PhGZjgz4I=">AAAB7HicbVBNS8NAEJ3Ur1q/qh69LBahp5KIoMeiF29WaNpCG8pmO2mXbjZhdyOU0t/gxYMiXv1B3vw3btsctPXBwOO9GWbmhang2rjut1PY2Nza3inulvb2Dw6PyscnLZ1kiqHPEpGoTkg1Ci7RN9wI7KQKaRwKbIfju7nffkKleSKbZpJiENOh5BFn1FjJv31o9mW/XHFr7gJknXg5qUCORr/81RskLItRGiao1l3PTU0wpcpwJnBW6mUaU8rGdIhdSyWNUQfTxbEzcmGVAYkSZUsaslB/T0xprPUkDm1nTM1Ir3pz8T+vm5noJphymWYGJVsuijJBTELmn5MBV8iMmFhCmeL2VsJGVFFmbD4lG4K3+vI6aV3WPLfmPV5V6tU8jiKcwTlUwYNrqMM9NMAHBhye4RXeHOm8OO/Ox7K14OQzp/AHzucPURSORA==</latexit>
sha1_base64="HlbaWngJl2+9LR1NU+PhGZjgz4I=">AAAB7HicbVBNS8NAEJ3Ur1q/qh69LBahp5KIoMeiF29WaNpCG8pmO2mXbjZhdyOU0t/gxYMiXv1B3vw3btsctPXBwOO9GWbmhang2rjut1PY2Nza3inulvb2Dw6PyscnLZ1kiqHPEpGoTkg1Ci7RN9wI7KQKaRwKbIfju7nffkKleSKbZpJiENOh5BFn1FjJv31o9mW/XHFr7gJknXg5qUCORr/81RskLItRGiao1l3PTU0wpcpwJnBW6mUaU8rGdIhdSyWNUQfTxbEzcmGVAYkSZUsaslB/T0xprPUkDm1nTM1Ir3pz8T+vm5noJphymWYGJVsuijJBTELmn5MBV8iMmFhCmeL2VsJGVFFmbD4lG4K3+vI6aV3WPLfmPV5V6tU8jiKcwTlUwYNrqMM9NMAHBhye4RXeHOm8OO/Ox7K14OQzp/AHzucPURSORA==</latexit><latexit
sha1_base64="HlbaWngJl2+9LR1NU+PhGZjgz4I=">AAAB7HicbVBNS8NAEJ3Ur1q/qh69LBahp5KIoMeiF29WaNpCG8pmO2mXbjZhdyOU0t/gxYMiXv1B3vw3btsctPXBwOO9GWbmhang2rjut1PY2Nza3inulvb2Dw6PyscnLZ1kiqHPEpGoTkg1Ci7RN9wI7KQKaRwKbIfju7nffkKleSKbZpJiENOh5BFn1FjJv31o9mW/XHFr7gJknXg5qUCORr/81RskLItRGiao1l3PTU0wpcpwJnBW6mUaU8rGdIhdSyWNUQfTxbEzcmGVAYkSZUsaslB/T0xprPUkDm1nTM1Ir3pz8T+vm5noJphymWYGJVsuijJBTELmn5MBV8iMmFhCmeL2VsJGVFFmbD4lG4K3+vI6aV3WPLfmPV5V6tU8jiKcwTlUwYNrqMM9NMAHBhye4RXeHOm8OO/Ox7K14OQzp/AHzucPURSORA==</latexit><latexit
sha1_base64="HlbaWngJl2+9LR1NU+PhGZjgz4I=">AAAB7HicbVBNS8NAEJ3Ur1q/qh69LBahp5KIoMeiF29WaNpCG8pmO2mXbjZhdyOU0t/gxYMiXv1B3vw3btsctPXBwOO9GWbmhang2rjut1PY2Nza3inulvb2Dw6PyscnLZ1kiqHPEpGoTkg1Ci7RN9wI7KQKaRwKbIfju7nffkKleSKbZpJiENOh5BFn1FjJv31o9mW/XHFr7gJknXg5qUCORr/81RskLItRGiao1l3PTU0wpcpwJnBW6mUaU8rGdIhdSyWNUQfTxbEzcmGVAYkSZUsaslB/T0xprPUkDm1nTM1Ir3pz8T+vm5noJphymWYGJVsuijJBTELmn5MBV8iMmFhCmeL2VsJGVFFmbD4lG4K3+vI6aV3WPLfmPV5V6tU8jiKcwTlUwYNrqMM9NMAHBhye4RXeHOm8OO/Ox7K14OQzp/AHzucPURSORA==</latexit><latexit
<latexit
sha1_base64="7drEBFSbqF4fju5mWGHm87YOfYI=">AAAB7HicbVBNS8NAEJ34WetX1aOXxSL0VBIR9Fjw4s0KTVtoQ9lsJ+3SzSbsboQS+hu8eFDEqz/Im//GbZuDtj4YeLw3w8y8MBVcG9f9djY2t7Z3dkt75f2Dw6PjyslpWyeZYuizRCSqG1KNgkv0DTcCu6lCGocCO+Hkbu53nlBpnsiWmaYYxHQkecQZNVbyWw/NgRxUqm7dXYCsE68gVSjQHFS++sOEZTFKwwTVuue5qQlyqgxnAmflfqYxpWxCR9izVNIYdZAvjp2RS6sMSZQoW9KQhfp7Iqex1tM4tJ0xNWO96s3F/7xeZqLbIOcyzQxKtlwUZYKYhMw/J0OukBkxtYQyxe2thI2poszYfMo2BG/15XXSvqp7bt17vK42akUcJTiHC6iBBzfQgHtogg8MODzDK7w50nlx3p2PZeuGU8ycwR84nz9mjI5S</latexit>
sha1_base64="7drEBFSbqF4fju5mWGHm87YOfYI=">AAAB7HicbVBNS8NAEJ34WetX1aOXxSL0VBIR9Fjw4s0KTVtoQ9lsJ+3SzSbsboQS+hu8eFDEqz/Im//GbZuDtj4YeLw3w8y8MBVcG9f9djY2t7Z3dkt75f2Dw6PjyslpWyeZYuizRCSqG1KNgkv0DTcCu6lCGocCO+Hkbu53nlBpnsiWmaYYxHQkecQZNVbyWw/NgRxUqm7dXYCsE68gVSjQHFS++sOEZTFKwwTVuue5qQlyqgxnAmflfqYxpWxCR9izVNIYdZAvjp2RS6sMSZQoW9KQhfp7Iqex1tM4tJ0xNWO96s3F/7xeZqLbIOcyzQxKtlwUZYKYhMw/J0OukBkxtYQyxe2thI2poszYfMo2BG/15XXSvqp7bt17vK42akUcJTiHC6iBBzfQgHtogg8MODzDK7w50nlx3p2PZeuGU8ycwR84nz9mjI5S</latexit><latexit
sha1_base64="7drEBFSbqF4fju5mWGHm87YOfYI=">AAAB7HicbVBNS8NAEJ34WetX1aOXxSL0VBIR9Fjw4s0KTVtoQ9lsJ+3SzSbsboQS+hu8eFDEqz/Im//GbZuDtj4YeLw3w8y8MBVcG9f9djY2t7Z3dkt75f2Dw6PjyslpWyeZYuizRCSqG1KNgkv0DTcCu6lCGocCO+Hkbu53nlBpnsiWmaYYxHQkecQZNVbyWw/NgRxUqm7dXYCsE68gVSjQHFS++sOEZTFKwwTVuue5qQlyqgxnAmflfqYxpWxCR9izVNIYdZAvjp2RS6sMSZQoW9KQhfp7Iqex1tM4tJ0xNWO96s3F/7xeZqLbIOcyzQxKtlwUZYKYhMw/J0OukBkxtYQyxe2thI2poszYfMo2BG/15XXSvqp7bt17vK42akUcJTiHC6iBBzfQgHtogg8MODzDK7w50nlx3p2PZeuGU8ycwR84nz9mjI5S</latexit><latexit
sha1_base64="7drEBFSbqF4fju5mWGHm87YOfYI=">AAAB7HicbVBNS8NAEJ34WetX1aOXxSL0VBIR9Fjw4s0KTVtoQ9lsJ+3SzSbsboQS+hu8eFDEqz/Im//GbZuDtj4YeLw3w8y8MBVcG9f9djY2t7Z3dkt75f2Dw6PjyslpWyeZYuizRCSqG1KNgkv0DTcCu6lCGocCO+Hkbu53nlBpnsiWmaYYxHQkecQZNVbyWw/NgRxUqm7dXYCsE68gVSjQHFS++sOEZTFKwwTVuue5qQlyqgxnAmflfqYxpWxCR9izVNIYdZAvjp2RS6sMSZQoW9KQhfp7Iqex1tM4tJ0xNWO96s3F/7xeZqLbIOcyzQxKtlwUZYKYhMw/J0OukBkxtYQyxe2thI2poszYfMo2BG/15XXSvqp7bt17vK42akUcJTiHC6iBBzfQgHtogg8MODzDK7w50nlx3p2PZeuGU8ycwR84nz9mjI5S</latexit><latexit
<latexit
sha1_base64="sH7bVKd8y7W60UxbzdAv02sCNcg=">AAAB6nicbVBNS8NAEJ3Ur1q/qh69LBahp5KIoMeCF48VTVtoQ9lsJ+3SzSbsboQS+hO8eFDEq7/Im//GbZuDtj4YeLw3w8y8MBVcG9f9dkobm1vbO+Xdyt7+weFR9fikrZNMMfRZIhLVDalGwSX6hhuB3VQhjUOBnXByO/c7T6g0T+SjmaYYxHQkecQZNVZ6GA/koFpzG+4CZJ14BalBgdag+tUfJiyLURomqNY9z01NkFNlOBM4q/QzjSllEzrCnqWSxqiDfHHqjFxYZUiiRNmShizU3xM5jbWexqHtjKkZ61VvLv7n9TIT3QQ5l2lmULLloigTxCRk/jcZcoXMiKkllClubyVsTBVlxqZTsSF4qy+vk/Zlw3Mb3v1VrVkv4ijDGZxDHTy4hibcQQt8YDCCZ3iFN0c4L86787FsLTnFzCn8gfP5A0fyjbM=</latexit>
sha1_base64="sH7bVKd8y7W60UxbzdAv02sCNcg=">AAAB6nicbVBNS8NAEJ3Ur1q/qh69LBahp5KIoMeCF48VTVtoQ9lsJ+3SzSbsboQS+hO8eFDEq7/Im//GbZuDtj4YeLw3w8y8MBVcG9f9dkobm1vbO+Xdyt7+weFR9fikrZNMMfRZIhLVDalGwSX6hhuB3VQhjUOBnXByO/c7T6g0T+SjmaYYxHQkecQZNVZ6GA/koFpzG+4CZJ14BalBgdag+tUfJiyLURomqNY9z01NkFNlOBM4q/QzjSllEzrCnqWSxqiDfHHqjFxYZUiiRNmShizU3xM5jbWexqHtjKkZ61VvLv7n9TIT3QQ5l2lmULLloigTxCRk/jcZcoXMiKkllClubyVsTBVlxqZTsSF4qy+vk/Zlw3Mb3v1VrVkv4ijDGZxDHTy4hibcQQt8YDCCZ3iFN0c4L86787FsLTnFzCn8gfP5A0fyjbM=</latexit><latexit
sha1_base64="sH7bVKd8y7W60UxbzdAv02sCNcg=">AAAB6nicbVBNS8NAEJ3Ur1q/qh69LBahp5KIoMeCF48VTVtoQ9lsJ+3SzSbsboQS+hO8eFDEq7/Im//GbZuDtj4YeLw3w8y8MBVcG9f9dkobm1vbO+Xdyt7+weFR9fikrZNMMfRZIhLVDalGwSX6hhuB3VQhjUOBnXByO/c7T6g0T+SjmaYYxHQkecQZNVZ6GA/koFpzG+4CZJ14BalBgdag+tUfJiyLURomqNY9z01NkFNlOBM4q/QzjSllEzrCnqWSxqiDfHHqjFxYZUiiRNmShizU3xM5jbWexqHtjKkZ61VvLv7n9TIT3QQ5l2lmULLloigTxCRk/jcZcoXMiKkllClubyVsTBVlxqZTsSF4qy+vk/Zlw3Mb3v1VrVkv4ijDGZxDHTy4hibcQQt8YDCCZ3iFN0c4L86787FsLTnFzCn8gfP5A0fyjbM=</latexit><latexit
sha1_base64="sH7bVKd8y7W60UxbzdAv02sCNcg=">AAAB6nicbVBNS8NAEJ3Ur1q/qh69LBahp5KIoMeCF48VTVtoQ9lsJ+3SzSbsboQS+hO8eFDEq7/Im//GbZuDtj4YeLw3w8y8MBVcG9f9dkobm1vbO+Xdyt7+weFR9fikrZNMMfRZIhLVDalGwSX6hhuB3VQhjUOBnXByO/c7T6g0T+SjmaYYxHQkecQZNVZ6GA/koFpzG+4CZJ14BalBgdag+tUfJiyLURomqNY9z01NkFNlOBM4q/QzjSllEzrCnqWSxqiDfHHqjFxYZUiiRNmShizU3xM5jbWexqHtjKkZ61VvLv7n9TIT3QQ5l2lmULLloigTxCRk/jcZcoXMiKkllClubyVsTBVlxqZTsSF4qy+vk/Zlw3Mb3v1VrVkv4ijDGZxDHTy4hibcQQt8YDCCZ3iFN0c4L86787FsLTnFzCn8gfP5A0fyjbM=</latexit><latexit
<latexit
hn

BOTn
T OPn

M9 nM ST ` M9 nADV ` M9 nDSP ` M9 nSSM ` M9 nSRC ` M9 nIST ` M9 nF M I ` M9 nAP T “ 0,

n and m. Fn,m has dimensions of M {T and is positive when fow is from cell m and into cell n.
where M9 nM ST is the rate of change of solute mass in the cell due to storage, sorption, and decay; M9 nADV ,

For the advection term, the net rate at which solute mass is entering or leaving cell n is the sum of the
a simple “change in storage is equal to infow minus outfow” equation such that the addition of solute to a
the terms in equation 2–9 have dimensions of M {T . The sign convention is chosen for equation 2–9 using
elevations are T OPn and BOTn , respectively; and hn is the head (water level) in the cell. Above the water level, the cell is

cell is positive and the subtraction (removal) of solute from a cell is negative. When expressed in this form,
an immobile domain. Equation 2–9 contains two additional terms, M9 nF M I and M9 nAP T , which are not repre-

a cell to account for errors in the fow solution and exchange with an advanced package, respectively. All of

the M9 nM ST term is the net rate at which solute mass enters the cell from groundwater storage and the sorbed

ADV is the advective fow rate of solute

sented in equations 2–1 and 2–2. These terms represent the addition or subtraction of solute mass to or from

mass (M {T ) into cell n from cell m (the rates are positive for fow into cell n from cell m). Alternative for-
Figure 2–1. Porous media model cell that is partially saturated. The area of cell n in plan view is An ; the cell top and bottom

The advection and dispersion terms in equation 2–9 involve the transfer of solute mass between adjacent
and mixing from external groundwater fuid sources and sinks, respectively; M9 nSRC is the rate of solute mass

phase and due to decay or production. A negative value for M9 nM ST indicates a net uptake of solute mass into

model cells. These terms require an expression for the fow of solute mass, Fn,m , between two adjacent cells,
M9 nDSP , M9 nSSM are the net rates at which solute mass fows into or out of the cell due to advection, dispersion,

loading added directly to a cell; and M9 nIST is the rate of change of solute mass in the cell due to exchange with

(2–10)
(2–9)
 Chapter 2. Formulation and Solution of the Control-Volume Finite-Difference Equation 2–5

The M9 nDSP term of equation 2–9 is the net rate at which solute mass is entering or leaving cell n due to
hydrodynamic dispersion. This term is the sum of the hydrodynamic dispersive fows of solute mass between
cell n and each of its surrounding cells:

ÿ
M9 nDSP “ DSP
Fn,m , (2–11)
mPηn

where Fn,m
DSP is the dispersive fow rate of solute mass (M {T ) into cell n from cell m. Alternative formula-

tions for Fn,m

DSP are discussed in chapter 4.

All of the terms in equation 2–9 are discussed in more detail in subsequent chapters of this report.

Numerical Solution
Solution of the system of solute transport equations relies on the MODFLOW 6 framework described by
Hughes and others (2017), which is customized to solve generalized CVFD equations in which a cell can be
connected to any number of surrounding cells. The CVFD equation for solute transport can be written for
model cell n as

ÿ
An,n Cn ` An,m Cm “ bn , (2–12)
mPηn

where An,n is the coeffcient (L3 {T ) for the concentration in cell n, An,m is the coeffcient (L3 {T ) for the con-
centration in cell m, a neighbor of cell n, Cn and Cm are the concentrations (M {L3 ) in cells n and m, respec-
tively, and bn is the right-hand-side value of the balance equation (M {T ). The summation term in equation 2–
12 is written in a general way to indicate that the balance equation for cell n may depend on the concentra-
tions in any number of surrounding cells. The set of cells surrounding cell n is denoted by ηn . The terms in
equation 2–12 are assembled by the GWT Model piece-by-piece as each package adds contributions to the
coeffcients and right-hand-side terms. Once all of the assembly routines are complete, the concentration coef-
fcients and the right-hand-side term in equation 2–12 are the sum of the contributions from the different pack-
ages. For example, the An,m term may contain contributions from both the advection and dispersion packages.
These contributions are described in detail in subsequent chapters on the different transport packages.
The solute balance equation in the GWT Model is written using an implicit formulation in which the con-
centrations in equation 2–12 represent values at the end of the time step. This approach means that the concen-
tration values in equation 2–12 are unknown and must be solved simultaneously. The implicit formulation is
often preferred over an explicit formulation, because it is generally stable and allows for relatively large time
steps. Use of an explicit formulation would require a relatively small time step and fow expressions that use
known concentrations from the end of the previous time step so that equation 2–12 could be solved directly.
Equation 2–12 is the solute balance equation for a single model cell. Application of this balance equation
to every cell in the model grid results in a system of equations that can be expressed in matrix form as

AC “ b. (2–13)

In equation 2–13, A is a sparse square matrix with the number of rows and columns equal to the number of
cells. C is a vector of cell concentrations with the number of entries equal to the number of cells. b is the
right-hand-side vector also with the length equal to the number of cells. The number of cells here refers to the
number of cells in the model grid plus the number of features represented with advanced packages. Depending
2–6 Documentation for the MODFLOW 6 Groundwater Transport Model

on the options selected by the user, the A matrix may be symmetric or asymmetric, which affects the type of
linear solution methods that can be used. For some GWT Model applications, equation 2–13 is linear in that A
and b are independent of the dependent variable C. However, in many applications A and b depend on con-
centration, and the balance equation is, therefore, nonlinear. In such cases, solution of equation 2–13 must be
repeated multiple times, each time with updated coeffcients and right-hand-side values, until solution conver-
gence is achieved.

Initial Conditions
Starting concentrations are required as input for all cells in the model grid. These starting concentrations
represent the initial solute condition for the aquifer system. If one or more immobile domains are included
in the model, then starting concentrations for each immobile domain are also required. Starting concentra-
tions for the mobile domain are entered in the Initial Conditions (IC) Package. Starting concentrations for each
immobile domain are specifed in each IST Package.

Flow Model Interface

For most applications, the GWT Model is intended to work with a corresponding GWF Model. The GWF
Model calculates (1) heads for all active model cells, (2) fows between all connected model cells, (3) fows
between model cells and stress features, and (4) rates of storage changes in model cells. Flows calculated by
the GWF model represent volumetric fow rates with dimensions of L3 {T . The FMI Package is used within
the GWT Model to provide calculated groundwater fows. These calculated groundwater fows may come
from a corresponding GWF Model that is run concurrently in the same simulation as the GWT Model, or the
fows may be read from a fle that was created from a previous GWF simulation or from some other program
that is capable of writing fles in the correct format.

Flow Imbalance Correction

The GWT Model of MODFLOW 6 uses numerical methods to solve for the water balance of a model cell.
In most cases, there will be a small residual fow error Qe for a cell, which can be controlled through the toler-
ances specifed for the numerical solution. When Qe is zero, the infows and outfows for the cell are perfectly
balanced. If Qe is nonzero and fow into a model cell is greater than the outfow for the cell, then the calcu-
lated concentration may increase erroneously. Likewise, if the simulated fow out of a cell is larger than the
fow into a cell, then calculated concentrations may erroneously decrease. In most cases, these small errors are
inconsequential and go unnoticed. These errors can become more apparent in the results when simulated con-
centrations are larger or smaller than surrounding cell and boundary concentrations. In some situations, these
errors can grow over time and complicate interpretation of calculated concentrations.
Panday and others (2018) developed a simple approach to help minimize the effects of these fow imbal-
ance errors. The approach works by simply adding or removing solute mass at a rate equal to the product of
the residual fow error Qe and the calculated cell concentration,

M9 nF M I “ Qne Cn . (2–14)

In effect, the residual fow error is treated as a source or sink with the concentration equal to the calculated
cell concentration. To represent this addition or subtraction of solute mass in the system of equations, the A
coeffcient matrix is updated as
 Chapter 2. Formulation and Solution of the Control-Volume Finite-Difference Equation 2–7

An,n Ð An,n ´ Qen . (2–15)

For most solute transport applications, the fow model is solved with suffcient accuracy that the residual fow
error Qe for a cell is so small as to not have a noticeable effect on the calculated concentration. In some cases,
however, this optional correction available in the FMI Package can improve the accuracy of calculated con-
centrations, improve model stability, and minimize the presence of calculated concentrations that are above or
below expected values (Panday and others, 2018).

Time Stepping
For the present implementation of the GWT Model, all terms in the solute transport equation are solved
implicitly. With the implicit approach applied to the transport equation, it is possible to take relatively large
time steps and effciently obtain a stable solution. If the time steps are too large, however, accuracy of the
model results will decline, so there is usually some compromise required between the desired level of accu-
racy and length of the time step. If an explicit method were to be added to the GWT Model in the future, such
as the Method of Characteristics approach for solving the advection term, then time-step constraints would
likely be required to ensure numerical stability.
In MODFLOW 6, time step lengths are controlled by the user and specifed in the Temporal Discretiza-
tion (TDIS) input fle. When the fow model and transport model are included in the same simulation, then the
length of the time step specifed in TDIS is used for both models. If the GWT Model runs in a separate simula-
tion from the GWF Model, then the time steps used for the transport model can be different, and likely shorter,
than the time steps used for the fow solution. Instructions for specifying time steps are included in the input
and output guide that is distributed with the MODFLOW 6 software.

Special Considerations for Dry Cells

There are special considerations for handling dry model cells. The MODFLOW 6 GWF Model has two
alternatives for handling dry cells. With the traditional MODFLOW approach, a cell becomes dry when the
calculated head is below the cell bottom. When this drying happens, the cell is marked as dry and is removed
from fow calculations. In MODFLOW 6, a dry cell is assigned a special head value by the GWF Model. The
GWT Model checks for this dry cell value and, if present, removes the cell from the transport calculations.
Dry cells can be reactivated by the GWF Model if the rewetting option is selected by the user and if calcu-
lated heads in surrounding cells meet certain criteria (McDonald and others, 1992). The GWT Model checks
to see if a dry cell has rewet, and will reactivate that cell for transport calculations if that situation is encoun-
tered. This procedure for handling dry model cells is consistent with the approach implemented in MT3DMS
(Zheng, 1990; Zheng and Wang, 1999) and MT3D-USGS (Bedekar and others, 2016).
MODFLOW also has a Newton approach for handling dry cells (Painter and others, 2008), which
has become the preferred approach for simulating aquifers with unconfned conditions. With the Newton
approach, cells remain active, even if they are dry. Dry cells have calculated heads that are at or beneath the
cell bottom elevation. Although cells with heads below the cell bottom are dry, they can still transmit water to
underlying cells that are not dry. This vertical transmission happens instantaneously, as there is no water stor-
age within a dry cell. The Newton approach is available in MODFLOW-NWT (Niswonger and others, 2011),
MODFLOW-USG (Panday and others, 2013), and MODFLOW 6 (Langevin and others, 2017).
Bedekar and others (2016) implemented a method in MT3D-USGS to route solute mass through cells
with a calculated head below the cell bottom. The method works by inactivating these dry cells in the trans-
port model and accumulating mass fuxes in all overlying cells and then passing this mass directly to the high-
2–8 Documentation for the MODFLOW 6 Groundwater Transport Model

est cell that is not dry. The method has been shown to work for complicated groundwater modeling problems
involving many dry cells (Bedekar and others, 2016).
The MODFLOW 6 GWT Model handles transport through dry cells in a different manner than MT3D-
USGS. Instead of deactivating dry cells for solute transport, these cells remain active and a balance equation
is solved. Because there is no water in the cell, and thus, there is no solute mass in the cell, the balance equa-
tion reduces to a steady state form in which all solute mass coming into the cell instantaneously exits the cell.
Thus, for Newton models with dry cells, solute mass entering a dry cell is instantaneously transmitted down to
the uppermost cell that is not dry. Because the calculated head is below the cell bottom for these cells, the cal-
culated saturation is zero. Thus, whereas there may be an advective fux of solute through these model cells,
there is no dispersive fux as the dispersive fux equation contains a saturation term (eqn. 4–17). The advan-
tage of this approach is that solute transmission through dry cells is calculated as part of the matrix equations,
rather than requiring a separate accumulation step. A disadvantage of this approach is that the calculated con-
centrations in these dry cells are meaningless, and should be removed as part of a post-processing step. This
approach, as implemented in MODFLOW 6, has proven stable and robust for a variety of problems, including
transport through an aquifer system with perched conditions (Keating and Zyvoloski, 2009), and a complicated
water-table fuctuation problem reported by Langevin and others (2020).
 Chapter 3. Mobile Storage and Transfer 3–1

Chapter 3. Mobile Storage and Transfer

The Mobile Storage and Transfer (MST) Package of the GWT Model for MODFLOW 6 represents solute
mass storage, sorption, and frst- or zero-order decay. The mathematical expression for these terms is grouped
in equation 2–4 as f M ST . To isolate and describe these processes, the f M ST term is expanded as

f M ST “ f storage ` f sorption ` f decay ` f decaysorbed . (3–1)

This chapter describes the mathematical expressions for these individual terms and shows how the terms are
added to the system of equations. Because this chapter does not include equations describing fow between
cells, for example, between cells n and m, the n subscript is not included for cell variables.

Storage
The f storage term in equation 3–1 describes the rate of change of dissolved solute mass as

B pSw θCq
f storage “ ´ . (3–2)
Bt
The backward-in-time fnite-difference approximation of equation 3–2 for a model cell is

C t`∆t V t`∆t ´ C t Vwt

w
M9 storage “ ´ , (3–3)
∆t
where t is time (T ), ∆t is the length of the transport time step (T ), C is the solute concentration (M {L3 ), and
Vw is the volume of water in the cell (L3 ). Superscripts indicate the discrete time level at which each quan-
tity is evaluated, and the subscript n, which indicates that quantities M , C, and V pertain to cell n, has been
omitted for clarity. Equation 3–3 is derived by multiplying the f storage term in equation 3–2 by the cell vol-
ume (Vcell ), then using equation 2–8 to replace Sw θVcell with the volume of water in the cell, Vw . The expres-
sions in equations 3–2 and 3–3 are preceded by negative signs because the M9 storage term is on the left side of
the cell balance equation (eq. 2–9). Thus, if more mass is added to a cell than is lost (and the volume of water
remains the same), the concentration in the cell must increase.
The equation for aqueous solute storage is added to the system of equations by updating the diagonal posi-
tion of the A matrix for row n with the coeffcient of the C t`∆t term in equation 3–3 as

Vwt`∆t
An,n Ð An,n ´ , (3–4)
∆t
C t Vwt
and by moving the remaining ∆t term to the right-hand side to update bn as

C t Vwt
bn Ð bn ´ . (3–5)
∆t

Sorption
As a dissolved solute moves through an aquifer, some of the solute mass can bind to the surface of the
solid material, a process known as adsorption, or penetrate into the solid material, a process known as absorp-
3–2 Documentation for the MODFLOW 6 Groundwater Transport Model

tion (Zheng and Bennett, 2002). Adsorption and absorption, collectively called sorption, can slow the move-
ment of chemicals dissolved in groundwater. This apparent slowing of solute movement relative to water fow
is also known as retardation. Sorption can be modeled as a transfer of solute mass from a dissolved state in the
aqueous phase to a sorbed state in the solid phase, where it can no longer be transported by advection and dis-
persion. Desorption refers to the reverse process, whereby sorbed mass is released from the solid material and
reenters the fow system.
The sorption capability of the MST Package simulates sorption and desorption by transferring solute mass
between the groundwater and the solid material of the aquifer following the general approach implemented
by Zheng and Wang (1999) in the RCT Package of the MT3DMS program. The sorption capability described
here is similar to that of the MT3DMS RCT Package in that the MST Package also includes the linear equi-
librium sorption model and the Freundlich and Langmuir isotherms, but it differs in that it does not presently
include the nonequilibrium sorption model.
The concentration of sorbed mass, C, is expressed as the mass of sorbed solute per mass of solid aquifer
material available for sorption (M {M ). When all of the solid aquifer material is available for sorption, the
sorbed mass per volume of aquifer is ρb C, where ρb is the bulk density of the solid in M {L3 , the mass of
solid aquifer material per volume of aquifer. However, sorptive exchange between the mobile water and the
solid aquifer material requires contact between the mobile water and the pore walls. When one or more immo-
bile domains are present, only some fraction of the solid is available for sorptive exchange with the mobile
phase. Under fully saturated conditions, this fraction of porosity that is mobile is defned as fm , and the sorbed
mass per volume of aquifer is fm ρb C. In the absence of immobile domains, fm is internally assigned a value
of one as there is no other domain for sorption to occur. If one or more immobile domains are present, then
fm is calculated internally by the program as one minus the sum of the immobile domain porosities. Under
partially saturated conditions, the walls of desaturated pores are assumed to be effectively dry; any residual
flm of water that may remain on desaturated pore walls is assumed not to contribute appreciably to sorptive
exchange. The fraction of solid aquifer material available for sorption is assumed to be proportional to the
mobile water saturation, Sw , and the most general expression for sorbed mass per volume of aquifer is given
as fm ρb Sw C.
The rate of change of the sorbed mass per volume of aquifer is

` ˘
sorption B fm ρb Sw C
f “´ . (3–6)
Bt
Assuming fm and the bulk density do not change with time, equation 3–6 can be simplifed to

` ˘
sorption B CSw
f “ ´fm ρb . (3–7)
Bt
Equation 3–7 is a convenient expression in some cases, such as when a simple linear expression can be
used to relate the sorbed concentration to the aqueous concentration. For more complicated relations, however,
it is benefcial to use the product rule to expand equation 3–7 into

ˆ ˙
sorption BC BC BSw
f “ ´fm ρb Sw `C . (3–8)
BC Bt Bt

Equations 3–7 and 3–8 are discretized and used with the linear, Freundlich, and Langmuir isotherms to
develop fnite-difference approximations that include the effect of sorption.
 Chapter 3. Mobile Storage and Transfer 3–3

Linear Isotherm
There are several different conceptual models for relating the concentration of sorbed mass, C, to the con-
centration of dissolved solute, C. The most common approach is to assume that the dissolved solute is in equi-
librium with the sorbed solute, and that the concentration of sorbed mass is proportional to the concentration
of dissolved solute mass. This equilibrium-controlled linear approach can be represented mathematically by
the equation

C “ Kd C, (3–9)

where Kd is the linear distribution coeffcient (L3 {M ), often referred to as the partition coeffcient or the
adsorption ratio. This equilibrium-controlled linear sorption model is one of three options available in the
GWT Model for MODFLOW 6.
The sorption process is included in the system of equations by adding terms to the coeffcient matrix A
and the right-hand-side vector b. For the linear isotherm, the transfer of mass to and from the solid phase can
be approximated from the substitution of equation 3–9 into 3–7 and by multiplication by the volume of the
model cell to give

fm ρb Vcell fm ρb Vcell
M9 sorption “ ´ Kd pSw Cqt`∆t ` Kd pSw Cqt , (3–10)
∆t ∆t
where Kd is assumed to be constant. The subscript n, which indicates that the quantities M , fm , Vcell , Kd ,
Sw , and C pertain to cell n, has been omitted for clarity. The sorption mass transfer M9 sorption is part of the
M9 M ST term on the left side of equation 2–9. Thus, to place these terms into the matrix equation, the terms
on the right side of equation 3–10 that are coeffcients of C t`∆t are added to the diagonal position of the A
matrix for row n as

fm ρb Vcell
An,n Ð An,n ´ Kd Swt`∆t . (3–11)
∆t
The remaining terms must be moved to the right-hand-side vector, b, which requires that the signs of the terms
be changed:

fm ρb Vcell
bn Ð bn ´ Kd pSw Cqt . (3–12)
∆t

Freundlich Isotherm
The Freundlich isotherm is nonlinear with respect to the aqueous concentration and is written as

C “ Kf C a , (3–13)

where Kf is the Freundlich constant pL3 {M qa and a is the dimensionless Freundlich exponent. The Fre-
undlich isotherm is implemented in the MST Package by writing a fnite-difference expression for equation 3–
BC . Equation 3–13 can be differentiated with respect to the aqueous
8, which also requires an expression for BC
concentration to give
3–4 Documentation for the MODFLOW 6 Groundwater Transport Model

BC
“ aKf C a´1 . (3–14)
BC

The fnite-difference approximation to equation 3–8 is

ˆ ˙t` 12 ∆t
sorption fm ρb Vcell BC ` t`∆t
´ Ct
˘
M9 “´ Sw C
∆t BC
fm ρb Vcell ` ˘t` 12 ∆t ` t`∆t
´ Swt .
˘
´ C Sw (3–15)
∆t

To evaluate the t ` 12 ∆t terms an average value over the time step is used. For C and BC BC
, an average aqueous
concentration is calculated using the concentration at the start of the time step, and the most recent iterate of
the concentration at the end of the time step. This average aqueous concentration is then substituted into equa-
tions 3–13 and 3–14. Equation 3–15 is included in the system of equations by adding terms to the coeffcient
matrix A and the right-hand-side vector b.

Langmuir Isotherm

The Langmuir isotherm is also nonlinear with respect to the aqueous concentration and written as

Kl SC
C“ , (3–16)
1 ` Kl C

where Kl is the Langmuir constant pL3 {M q and S is the total concentration of sorption sites available
pM {M q. The Langmuir isotherm is implemented in a manner similar to the Freundlich isotherm by using the
fnite-difference approximation shown in equation 3–15. For the Langmuir isotherm, the derivative term is

BC Kl S
“ . (3–17)
BC p1 ` Kl Cq2

Considerations for Unconfned Conditions

The sorption formulation implemented in the GWT Model is mass conservative. A consequence of this
sorption formulation described here, however, is that simulated concentrations may be higher than or lower
than expected concentrations under unconfned conditions. When the water table decreases in a model cell,
sorbed mass within the dewatered part of the cell is instantaneously applied to the saturated part of the cell at
the end of the time step. Thus, in the absence of other solute infows and outfows, solute is distributed over a
smaller volume at the end of the time step than the volume at the start of the time step. In extreme cases this
may cause the solute concentration to increase above expected concentrations simply due to the reduction
in the volume of water in the cell. MT3D-USGS has options for treating this condition (Bedekar and others,
2016); however, these options are not available in the version of the GWT Model described in this report.
 Chapter 3. Mobile Storage and Transfer 3–5

Decay
The decay capability of the MST Package simulates the effects of frst- or zero-order decay and produc-
tion of the dissolved aqueous phase. Decay and production can also be represented for the sorbed phase, as
described in a following section (Decay of Sorbed Mass), and for the dissolved aqueous and sorbed phases in
an immobile domain as described in chapter 7, “Immobile Domain Storage and Transfer”. The mathematical
expression for decay,

f decay “ ´λ1 θSw C ´ γ1 θSw , (3–18)

can be written to include frst-order decay and zero-order decay, although only one can be active in the present
GWT Model implementation. λ1 is the frst-order decay rate coeffcient for the mobile domain (T ´1 ), and
γ1 is the zero-order decay rate coeffcient for the mobile domain (M L´3 T ´1 ). Implementation of decay is
handled differently depending on whether frst-order or zero-order is selected.

First Order Decay

First-order decay is typically used to represent radioactive decay and biodegradation or chemical decom-
position. When frst-order decay is represented, f decay can be written as

f decay “ ´λ1 θSw C. (3–19)

For frst-order decay defned by equation 3–19, the rate of solute-mass decay in a cell of volume Vcell is

M9 decay “ ´λ1 θVcell pSw Cqt`∆t . (3–20)

To include the effect of frst-order decay in the matrix equations, the coeffcient term on the right side of equa-
tion 3–20 (that is a coeffcient of the dependent variable C t`∆t ) is added to the diagonal position of the A
matrix for row n as

An,n Ð An,n ´ λ1 θVcell Swt`∆t . (3–21)

First-order decay is commonly written in terms of a half life t1{2 , which is the length of time for the solute
concentration to decrease by half. The half life is related to the frst-order decay rate coeffcient by

ln 2
t1{2 “ . (3–22)
λ1

Zero-Order Decay
Under some transport conditions, zero-order decay may be used to mathematically represent the process
of biodegradation in which the rate of biodegradation does not depend on concentration. Zero-order decay can
also be used to model groundwater age mathematically as a solute “concentration,” as described by Goode
(1996). When groundwater age is simulated in this way, a zero-order (constant) growth rate of one (a decay
rate of minus one) causes the groundwater age to increase by one unit of time for every unit of time simulated.
Thus, during each time step the groundwater age increases by the length of the time step.
3–6 Documentation for the MODFLOW 6 Groundwater Transport Model

When zero-order decay is represented, f decay is

f decay “ ´γ1 θSw , (3–23)

which is expressed in fnite-difference form as

M9 decay “ ´γ1 θVcell Swt`∆t . (3–24)

Unlike in the case of frst-order decay (eq. 3–20), the dependent variable C does not appear in the expres-
sion for zero-order decay (eq. 3–24), so the zero-order term is moved to the right-hand-side vector, b, which
requires that the sign of the term be changed:

bn Ð bn ` γ1 θVcell Swt`∆t . (3–25)

Zero-order decay can result in negative concentrations unless the decay rate is automatically reduced by
the program based on simulated solute concentrations. A simple method was implemented to reduce the user-
specifed zero-order decay rate if that rate would result in negative concentrations. This reduced rate is esti-
mated by comparing the user-specifed decay rate with the solute concentration divided by the length of the
time step. If the user-specifed decay rate is larger than the solute concentration divided by the time step, then
the program replaces γ1 with the solute concentration divided by the time step. Because this approach requires
the current concentration, which is estimated as the previous concentration iterate, additional outer iterations
may be required for solution convergence.

Decay of Sorbed Mass

When sorption is represented in the GWT Model, separate terms are included in the mass balance equation
to represent decay of the sorbed mass. Decay of the sorbed mass can be expressed with the following mathe-
matical expression (Toride and others, 1993),

f sorbeddecay “ ´λ2 fm ρb Sw C ´ γ2 fm ρb Sw . (3–26)

Linear Isotherm
For frst-order decay with the linear sorption isotherm, the rate of decay of sorbed mass in a cell is

M9 sorbeddecay “ ´λ2 fm ρb Sw Kd Vcell C t`∆t , (3–27)

and for zero-order decay, the rate of decay of sorbed mass in a cell is

M9 sorbeddecay “ ´γ2 fm ρb Sw Vcell . (3–28)

The decay of sorbed mass is added to the system of equations differently depending on whether the decay
is frst or zero order. For frst-order decay, in which the mass transfer rate depends on concentration, the coeff-
cient of C t`∆t in equation 3–27 is added to the diagonal position of the A matrix for row n as
 Chapter 3. Mobile Storage and Transfer 3–7

An,n Ð An,n ´ λ2 fm ρb Sw Kd Vcell . (3–29)

For zero-order decay, in which the mass transfer rate is independent of solute concentration, the mass transfer
rate is moved to the right-hand-side vector, b, as

bn Ð bn ` γ2 fm ρb Sw Vcell . (3–30)

Additional precautions are also necessary, as described in the previous section, to ensure that zero-order decay
does not result in negative concentrations.

Nonlinear Freundlich and Langmuir Isotherms

For the nonlinear Freundlich and Langmuir isotherms, zero-order decay is also represented using equa-
tions 3–28 and 3–30 as the rate is not a function of the sorbed concentration (except when the rate must be
reduced to prevent negative concentrations). For frst-order decay, however, sorbed mass is expressed using C
as

t`∆t
M9 sorbeddecay “ ´λ2 fm ρb Sw Vcell C , (3–31)
t`∆t
where C is approximated using the isotherm equation (eq. 3–13 for the Freundlich isotherm or 3–16
for the Langmuir isotherm) and the most recent iterate for the aqueous concentration. Unlike for the linear
isotherm with frst-order decay, which can be added to the A matrix, the term resulting from frst-order decay
with the nonlinear isotherms must be added to the right-hand-side vector. Additional outer iterations are often
required for solution convergence when the Freundlich isotherm or the Langmuir isotherm is used to represent
sorption.
 Chapter 4. Advective and Dispersive Solute Transport 4–1

Chapter 4. Advective and Dispersive Solute Transport

The Advection (ADV) and Dispersion (DSP) Packages add terms to the system of equations to represent
advective transport and the combined effects of mechanical dispersion and molecular diffusion. These two
transport processes are represented in the GWT Model by transferring mass between adjacent model cells
according to groundwater fow rates and solute concentration gradients.

Advection (ADV) Package

Advection is the movement of a dissolved solute as it is transported through an aquifer at the average
linear velocity of the groundwater fow. The advective fux f ADV (M {L2 T ) of a solute of concentration C
(M {L3 ) transported by a specifc discharge of groundwater q (L{T ) is

f ADV “ qC. (4–1)

This advective fux term comprises the advective solute fux that appears in parentheses in the frst term on the
right-hand side of equation 2–1. As indicated in equation 2–10, the discrete analog of equation 4–1 used by
the CVFD method in the GWT Model expresses the total advective fow rate of solute mass across the inter-
face between cells n and m, Fn,mADV (M {T ), as the product of the volumetric groundwater fow rate across the

interface, Qn,m (L3 {T ), and a representative solute concentration of groundwater crossing the interface, Cn,m
(M {L3 ):

ADV
Fn,m “ Qn,m Cn,m , (4–2)

where Fn,m
ADV and Q
n,m are defned to be positive for fow into cell n from cell m.
Three options are available for calculating the concentration at the cell face, Cn,m : central-in-space
weighting, upstream weighting, and a total variation diminishing (TVD) scheme. In each case, the solute con-
centration at the cell face can be expressed as a weighted average of the concentrations in the two cells:

Cn,m “ ωn,m Cn ` p1 ´ ωn,m q Cm , (4–3)

where ωn,m is a dimensionless weighting factor.

Central-In-Space Weighting
The central-in-space weighting scheme is based on a simple distance-weighted linear interpolation
between the center of cell n and the center of cell m to calculate solute concentration at the shared face
between cell n and cell m. Although “central-in-space” is a misnomer for grids without equal spacing between
connected cells, it is retained here for consistency with nomenclature used by other MODFLOW-based trans-
port programs, such as MT3D. The value for ωn,m is a distance-weighted interpolation factor calculated based
on cell dimensions:

Lm,n
ωn,m “ , (4–4)
Ln,m ` Lm,n
4–2 Documentation for the MODFLOW 6 Groundwater Transport Model

where Ln,m and Lm,n are the distance (L) from the center of cell n to its shared face with cell m and the dis-
tance from the center of cell m to its shared face with cell n, respectively. Central-in-space weighting is not
often used because it can result in spurious oscillations in the simulated concentrations. It is included as an
option in the ADV Package, however, because it may be useful for model testing and comparison purposes.

Upstream Weighting
Upstream weighting is a commonly used approach for calculating the interface concentration. For
upstream weighting, the weighting factor, ωn,m , depends on the sign of the fow between cell n and m accord-
ing to

$
&0 for Qn,m ą 0
ωn,m “ . (4–5)
%1 for Qn,m ĺ 0

Total Variation Diminishing (TVD)

A TVD expansion of Cn,m can be expressed as

ups TV D
Cn,m “ Cn,m ` Cn,m , (4–6)

with

TV D σn,m ´ dwn ups

¯
Cn,m “ Cn,m ´ Cn,m , (4–7)
2
ups
where Cn,m is the concentration of the upstream cell, Cn,m
dwn is the concentration of the downstream cell,
T V D is an additional fux-limiting concentration term, and σ
Cn,m n,m is the fux limiter. The van Leer fux lim-
iter (Forsyth and others, 1998) can be written as

$
&0 for rn,m ĺ 0
σn,m “ 2rn,m
(4–8)
%
1`rn,m for rn,m ą 0.

rn,m is the smoothness sensor defned as

ups 2up
Cn,m ´ Cn,m Ln,m ` Lm,n
rn,m “ ¨ dwn ups , (4–9)
Lups,2up ` L2up,ups Cn,m ´ Cn,m

2up
where Cn,m is the concentration of the second upstream cell, and Lups,2up and L2up,ups are the distance (L)
from the center of the upstream cell to its shared face with the second upstream cell and the distance from the
center of the second upstream cell to its shared face with the upstream cell, respectively. The second upstream
cell is determined by identifying the cell with the largest fow into the upstream cell. If a second upstream cell
cannot be identifed because of fow conditions or because the upstream cell is on the edge of the grid, then the
TVD adjustment is not applied.
 Chapter 4. Advective and Dispersive Solute Transport 4–3

As shown in equation 4–6 the interface concentration Cn,m consists of two terms for the TVD expan-
ups
sion. The frst term Cn,m is simply an upstream-weighted concentration calculated using equations 4–3 and
4–5. Thus, the upstream weighting factor ωn,m (eq. 4–5) enters the advective fow expression, equation 4–
ups
2, through the Cn,m term in equation 4–6. The contribution to the advective fow associated with Cn,m
T V D is

treated separately, as described below in the discussion of the numerical solution procedure. However, for the
purpose of comparison with other weighting schemes, the TVD scheme can be expressed entirely in terms of
its own weighting factor, ωn,mT V D , which depends on the fux limiter, σ
n,m :

TV D TV D
` ˘
Cn,m “ ωn,m Cn ` 1 ´ ωn,m Cm , (4–10)

where

$
& σn,m for Qn,m ą 0
TV D
ωn,m “ 2 . (4–11)
σn,m
%1 ´
2 for Qn,m ĺ 0

Numerical Solution
When the central-in-space or upstream weighting scheme is used, the advective-transport term for the face
shared by cells n and m is incorporated into the system of linear equations, equation 2–13, as follows. Substi-
tution of equation 4–3 into equation 4–2 gives

ADV
Fn,m “ ωn,m Qn,m Cn ` p1 ´ ωn,m q Qn,m Cm . (4–12)

The matrix coeffcient in the diagonal position for row n is updated by adding the term ωn,m Qn,m as

An,n Ð An,n ` ωn,m Qn,m . (4–13)

The matrix coeffcient in row n that corresponds to the connection between cell n and neighboring cell m is
updated by adding the coeffcient part of the term p1 ´ ωn,m q Qn,m Cm as

An,m Ð An,m ` p1 ´ ωn,m q Qn,m . (4–14)

When the TVD option is used, the expression for Fn,m

ADV includes the additional TVD term Q
n,m Cn,m as
TV D

ADV TV D
Fn,m “ ωn,m Qn,m Cn ` p1 ´ ωn,m q Qn,m Cm ` Qn,m Cn,m , (4–15)

where the weighting factor ωn,m is evaluated as for upstream weighting (equation 4–5), and Cn,m
T V D is the fux-

limiting concentration introduced in equation 4–6. The terms in equation 4–15 that involve ωn,m are incorpo-
rated in the coeffcient matrix just as they would be for upstream weighting. The TVD term is added to bn , the
row-n entry in the right-hand-side vector, b:

TV D
bn Ð bn ´ Qn,m Cn,m . (4–16)
4–4 Documentation for the MODFLOW 6 Groundwater Transport Model

When upstream weighting or the TVD scheme is used, the dependence of ωn,m on the fow direction
causes the coeffcient matrix, A, to be nonsymmetric, that is, An,m ‰ Am,n . Solution of the resulting matrix
problem requires use of a linear solver that can accommodate nonsymmetric matrices.

Dispersion (DSP) Package

Over time, a mass of solute dissolved in groundwater spreads by molecular diffusion. If groundwater is
fowing, solute spreads also by mechanical dispersion, which is caused by “[v]ariations in local velocity, both
in magnitude and direction, along the tortuous paths and between adjacent fow paths” within the pores of the
porous medium (Bear, 1972). The combination of molecular diffusion and mechanical dispersion is called
hydrodynamic dispersion.
Although mechanical dispersion originates from pore-scale processes, it is generally not feasible to sim-
ulate its effects at the pore scale. Mathematical representation of mechanical dispersion at the macroscopic
scale, which is of interest in practical applications, has been the subject of active research for decades, and
various “macrodispersion” models have been developed (for example, Scheidegger (1961), Burnett and Frind
(1987), and Lichtner and others (2002)).
The simple “Fickian” model of hydrodynamic dispersion conceptualizes mechanical dispersion as being
analogous to molecular diffusion. In the Fickian model, the fux of solute mass due to hydrodynamic disper-
sion, f DSP (M {L2 T ), is related to the solute concentration gradient, ∇C (M {L4 ), through a tensorial “con-
stant of proportionality,” the hydrodynamic dispersion tensor, D (L2 {T ), which is the sum of contributions
from mechanical dispersion and molecular diffusion:

f DSP “ ´Sw θD∇C, (4–17)

with

D “ Dmech ` Dmol I, (4–18)

where Dmech (L2 {T ) is the mechanical dispersion tensor, which may be anisotropic, Dmol (L2 {T ) is an effec-
tive molecular diffusion coeffcient that takes into account the effect of porous medium tortuosity, I is the
identity tensor (dimensionless), and θ is effective porosity (dimensionless). Equation 4–17 represents the dis-
persive solute fux that appears in parentheses in the second term on the right-hand side of equation 2–1. The
second term on the right-hand side of equation 4–18 places the contribution from the molecular diffusion coef-
fcient on the diagonal of the hydrodynamic dispersion tensor, which corresponds to assuming that molecular
diffusion is isotropic. Following Zheng and Wang (1999), θ is included explicitly in equation 4–17 to account
for the fact that hydrodynamic dispersion occurs only within the pore space, not within the solid matrix. Sim-
ilarly, inclusion of the saturation, Sw , in equation 4–17 accounts for the fact that hydrodynamic dispersion
occurs only within the portion of the porespace that contains water.
Equation 4–17 shows the dispersive fux as proportional to the gradient of the solute concentration. This
relation is consistent with the dispersion model presented by Zheng and Bennett (2002) and implemented in
popular solute transport modeling codes such as MODFLOW-GWT and MT3DMS. Other literature suggests,
however, that the dispersive fux should be proportional to the gradient of solute mass fraction (Bird and oth-
ers, 2006), which is the dispersion model implemented in the SUTRA program (Voss and Provost, 2010). For
most transport problems with slight density variations, the difference between the two approaches is expected
to be negligible.
The mechanical dispersion tensor, Dmech , is typically assumed to be characterized by three mutually per-
pendicular “principal” directions of spreading, which implies that the tensor is defned by a real, symmetric
 Chapter 4. Advective and Dispersive Solute Transport 4–5

matrix. In a homogeneous aquifer, a solute mass that is initially spherical will assume an ellipsoidal shape as
it advects with the groundwater fow and spreads due to mechanical dispersion, and the principal axes of the
ellipsoid will coincide with the principal directions of the dispersion tensor. The entries in the matrix, or “dis-
persion coeffcients,” control the rates and directions of spreading. The dispersion coeffcients can vary with
fow direction, and different models for the dispersion coeffcients lead to different symmetries in the solute
spreading pattern. For example, “isotropic” dispersion is controlled by two dispersion coeffcients that do not
vary with fow direction, and the rates of spreading along and perpendicular to the fow direction are inde-
pendent of the fow direction. The model of Scheidegger (1961), which in its most general form admits non-
symmetric matrices and is defned by 61 dispersion coeffcients, encompasses a wide variety of symmetries.
In practice, however, the mechanical dispersion model is typically simplifed considerably by assuming that
one of the principal directions of the dispersion tensor, called the “longitudinal” direction, is always aligned
with the fow direction. The remaining two principal directions, called the “transverse” directions, are then
perpendicular to the fow direction. Special cases of this “fow-aligned” dispersion tensor have been devel-
oped by incorporating additional simplifying assumptions. For example, the popular model of Burnett and
Frind (1987), which is defned by two longitudinal and three transverse dispersion coeffcients, allows differ-
ent spreading behavior for vertical fow than for fow within the horizontal plane and includes isotropic disper-
sion as a special case. For fow within the horizontal plane, the rates of longitudinal and transverse dispersion
are independent of fow direction, but the rate of horizontal transverse dispersion can be different from the rate
of vertical transverse dispersion. For vertical fow, transverse dispersion is horizontal and is axially symmetric
about the vertical direction. For intermediate fow directions, dispersion coeffcients are interpolated a specifc
way between their horizontal-fow and vertical-fow values.
The Fickian model described by equation 4–17 tends to produce unrealistic back dispersion, and the dis-
persion coeffcients are scale-dependent, particularly in applications involving estimation of effective dis-
persion coeffcients in heterogeneous aquifers (Konikow, 2010). Nevertheless, conceptual simplicity, ease
of integration into conventional groundwater fow and transport models, and the ability to simulate varying
degrees and directions of solute spreading make the Fickian model a popular tool for representing mechanical
dispersion in practical applications. The GWT Model of MODFLOW 6 offers a Fickian dispersion model that
includes isotropic dispersion and the mechanical dispersion models of Burnett and Frind (1987) and Lichtner
and others (2002) as special cases.
The CVFD method used in the GWT Model is based on a solute mass balance over each cell (equation 2–
9), which includes dispersive fows of solute between each cell and its surrounding cells (equation 2–11). The
GWT Model offers two choices for formulating Fn,m DSP , the dispersive fow rate of solute mass (M {T ) into cell

n from cell m across their shared face, which is a discrete analog of equation 4–17. The “simplifed” formu-
lation is mathematically analogous to the “conductance-based” formulation of groundwater fow across cell
faces in MODFLOW 6. Similar to the conductance-based groundwater fow formulation, the simplifed solute-
mass fow formulation works well when the model grid and governing tensor—in this case the dispersion ten-
sor—satisfy certain requirements, as discussed in Langevin and others (2017) and summarized below. When
these requirements are not met use of the simplifed formulation introduces numerical error (in addition to the
usual discretization error), which may or may not be signifcant in a given application. In such cases the alter-
native XT3D formulation of the solute-mass fow, which is mathematically analogous to the XT3D formula-
tion of groundwater fow in MODFLOW 6 (Provost and others, 2017), can be helpful because it automatically
accounts for irregularities in the grid and tensor anisotropy, albeit at the expense of longer simulation times.
The XT3D formulation is used by default for solute-mass fow in MODFLOW 6, but the simplifed fow for-
mulation can be activated by the user to evaluate model performance and accuracy.
When the GWT Model is used with a GWF Model that uses the Horizontal Flow Barrier (HFB) Pack-
age, transport results should be evaluated with caution as noted by Hornberger and others (2002), especially
if the intent of the fow barrier is to impede solute movement. Flow barriers are not explicitly represented in
the GWF Model. Instead, their effects are implicitly incorporated in the fow equations by adjusting the con-
4–6 Documentation for the MODFLOW 6 Groundwater Transport Model

ductance between two model cells to account for lower permeability material. A consequence of the implicit
incorporation of a fow barrier is that the simulated dispersive fux through the barrier will be too large as it
does not have any thickness or storage capacity. Therefore, the effectiveness of the barrier to contain a solute
plume or impede solute movement will be underrepresented by the GWT Model in this situation.

Mechanical Dispersion Model

In the MODFLOW 6 GWT Model, the mechanical dispersion tensor, Dmech , is assumed to be character-
ized by three mutually perpendicular “principal” directions of spreading, one of which, called the “longitudi-
nal” direction, is always aligned with the direction of groundwater fow. The remaining two principal direc-
tions, called the “transverse” directions, are then perpendicular to the fow direction and to each other. One of
the transverse directions lies within the (x, y) (horizontal) plane. Expressed in coordinates (xL , xT 1 , xT 2 ) that
align with the longitudinal and two transverse directions, respectively, the mathematical form of the mechani-
cal dispersion tensor used in the GWT Model is

¨ ˛
αL v 0 0
Dmech αT 1 v 0 ‚, (4–19)
˚ ‹
“˝ 0
0 0 αT 2 v

where αL , αT 1 , and αT 2 are called the longitudinal and frst and second transverse dispersivities (L), and v
is the groundwater fow velocity (L{T ) or “seepage velocity” calculated as the specifc discharge divided by
the porosity. Although the equations presented here are based on groundwater velocity, implementation of
these equations in the Dispersion Package is based on specifc discharge (instead of velocity) by including the
porosity term in equation 4–17 in the dispersion tensor and by replacing vθ with q.
The longitudinal and transverse dispersivities vary with the groundwater fow direction as follows:

vz2 vz2
ˆ˙
2 2
αL “ αLH cos θ2 ` αLV sin θ2 “ αLH 1 ´ 2 ` αLV
v v2
vz2 vz2
ˆ ˙
2 2
αT 1 “ αT H1 cos θ2 ` αT V sin θ2 “ αT H1 1 ´ 2 ` αT V , (4–20)
v v2
vz2 vz2
ˆ ˙
2 2
αT 2 “ αT H2 cos θ2 ` αT V sin θ2 “ αT H2 1 ´ 2 ` αT V
v v2

where θ2 is the angle at which the groundwater velocity vector, v, is inclined upward from the (x, y) (horizon-
tal) plane, and vz is the z (vertical) component of v. The model has fve parameters that can be specifed by
the user: two longitudinal dispersivities, αLH and αLV , and three transverse dispersivities, αT H1 , αT H2 , and
αT V , all with dimensions of (L).
Note that equation 4–19 expresses Dmech in (xL , xT 1 , xT 2 ) coordinates, but velocity is expressed relative
to (x, y, z) model coordinates in equation 4–20. The “horizontal plane” relative to which the dispersion tensor
is defned is the (x, y) model coordinate plane, and the “vertical” direction is the model coordinate z direction.
Allowing the “horizontal” and “vertical” directions to be rotated relative to the model coordinates would allow
the defnition of the dispersion tensor to align with dipping beds; for example, in much the same way that three
user-specifed angles allow the conductivity tensor to be defned with respect to coordinates rotated relative to
(x, y, z) model coordinates (Langevin and others, 2017). However, an option to set the reference coordinates
for the dispersion tensor to be other than the model coordinates is not available in the current implemention.
 Chapter 4. Advective and Dispersive Solute Transport 4–7

The behavior of the mechanical dispersion model can be understood by considering the values assumed by
the longitudinal and transverse dispersivities when the groundwater fow is either horizontal or vertical. When
fow is in the horizontal plane (θ2 “ 00 ; vz “ 0), the longitudinal dispersivity (αL ) is αLH , the transverse dis-
persivity for horizontal spreading (αT 1 ) is αT H1 , and the transverse dispersivity for vertical spreading (αT 2 )
is αT H2 . When fow is in the vertical direction (θ2 “ 900 ; vx “ vy “ 0, vz “ v), the longitudinal dispersiv-
ity (αL ) is αLV , and the transverse dispersivities (αT 1 and αT H1 ), which represent horizontal spreading, are
both αT V . For groundwater fow directions between horizontal and vertical, the longitudinal and transverse
dispersivities vary smoothly with fow direction according to equation 4–20.
When transformed entirely into (x, y, z) model coordinates, the mechanical dispersion tensor has the fol-
lowing form:

vx2 vy2 v2
mech
Dxx “ αL ` αT H1 ` αT x z
v v v
2 v 2 2
v y v
mech
Dyy “ αT H1 x ` αL ` αT y z
v v v
2 2
v v y v2
mech
Dzz “ αT 2 x ` αT 2 ` αL z (4–21)
v v v
mech vx vy
Dxy “ pαL ´ αT 2 ` αT H2 ´ αT H1 q
v
mech vx vz
Dxz “ pαL ´ αT 2 q
v
mech vy vz
Dyz “ pαL ´ αT 2 q
v

where

v2 vx2
ˆ ˙
αT x “ αT H2 x2 ` αT V 1´ 2
v v
˜ ¸. (4–22)
vy2 vy2
αT y “ αT H2 ` αT V 1´ 2
v2 v

The generalized mechanical dispersion tensor, written in model coordinates (eq. 4–21), can be recast to
correspond to other popular dispersion models, such as the models of Burnett and Frind (1987), Lichtner and
others (2002), and the isotropic dispersion model. Setting

αLH “ αL
αLV “ αL
αT H1 “ αTH , (4–23)
αT H2 “ αTV
αT V “ αTV

where αL , αTH , and αTV are constants introduced for simplifcation purposes, gives the model of Burnett and
Frind (1987):
4–8 Documentation for the MODFLOW 6 Groundwater Transport Model

vx2 vy2 v2
mech
Dxx “ αL ` αTH ` αTV z
v v v
2 2
v vy v2
mech
Dyy “ αTH x ` αL ` αTV z
v v v
2 v 2 2
v y v
mech
Dzz “ αTV x ` αTV ` αL z . (4–24)
v v v
mech H v x vy
` ˘
Dxy “ αL ´ αT
v
mech V v x vz
` ˘
Dxz “ αL ´ αT
v
mech V v y vz
` ˘
Dyz “ αL ´ αT
v

Setting

H
αLH “ αL
V
αLV “ αL
αT H1 “ αTH , (4–25)
αT H2 “ αTV
αT V “ αTH

where αL
H , αV , αH , and αV are constants, gives a form of the model of Lichtner and others (2002):
L T T

vx2 vy2 v2
mech
Dxx “ αL ` αTH ` αT x z
v v v
2 2
v vy v2
mech
Dyy “ αTH x ` αL ` αT y z
v v v
2 v 2 2
v y v
mech
Dzz “ αT 2 x ` αT 2 ` αL z , (4–26)
v v v
mech H vx vy
` ˘
Dxy “ αL ´ αT 2 ` αT H2 ´ αT
v
mech vx vz
Dxz “ pαL ´ αT 2 q
v
mech vy vz
Dyz “ pαL ´ αT 2 q
v

with
 Chapter 4. Advective and Dispersive Solute Transport 4–9

vz2 2
ˆ ˙
H V vz
αL “ αL 1 ´ 2 ` αL
v v2
2 vx2
ˆ ˙
V vx H
αT x “ αT 2 ` αT 1 ´ 2
v v
2
˜ ¸. (4–27)
v
V y H
vy2
αT y “ αT 2 ` αT 1 ´ 2
v v
v2 v2
ˆ ˙
αT 2 “ αTV 1 ´ z2 ` αTH z2
v v

Setting

αLH “ αL
αLV “ αL
αT H1 “ αT , (4–28)
αT H2 “ αT
αT V “ α T

where αL and αT are constants, gives an isotropic mechanical dispersion tensor, which is a special case of
both the Burnett and Frind (1987) model and the Lichtner and others (2002) model:

vx2 vy2 v2
mech
Dxx “ αL ` αT ` αT z
v v v
2 2
v vy v2
mech
Dyy “ αT x ` αL ` αT z
v v v
2 v 2 2
v y v
mech
Dzz “ αT x ` αT ` αL z . (4–29)
v v v
mech vx vy
Dxy “ pαL ´ αT q
v
mech vx vz
Dxz “ pαL ´ αT q
v
mech vy vz
Dyz “ pαL ´ αT q
v
With the isotropic dispersion model, the rate of longitudinal spreading (along the fow direction) is controlled
by αL , and the rate of transverse spreading (symmetrically about an axis oriented with the fow direction) is
controlled by αT . The rates of longitudinal and transverse spreading are independent of the fow direction.
In general, an initially spherical patch of solute in a uniform fow feld spreads into a diffuse ellipsoid that is
symmetric about an axis oriented with the fow direction. In the special case αT “ αL , an initially spherical
patch of solute in a uniform fow feld spreads into a diffuse sphere.

Simplifed Formulation
Hydrodynamic dispersion results in solute spreading from areas of higher concentration to areas of lower
concentration. In an anisotropic porous medium, however, dispersive transport is not necessarily in the direc-
4–10 Documentation for the MODFLOW 6 Groundwater Transport Model

tion of the concentration gradient; that is, the solute mass fux vector is not necessarily aligned with the con-
centration gradient vector. In the Fickian dispersion model (equation 4–17) in three dimensions (3D), each of
the three components of the solute mass fux vector is related, by way of the hydrodynamic dispersion tensor,
to each of the three components of the concentration gradient:

ˆ ˙
BC BC BC
fxDSP “ ´Sw θ Dxx ` Dxy ` Dxz
Bx By Bz
ˆ ˙
DSP BC BC BC
fy “ ´Sw θ Dxy ` Dyy ` Dyz , (4–30)
Bx By Bz
ˆ ˙
BC BC BC
fzDSP “ ´Sw θ Dxz ` Dyz ` Dzz
Bx By Bz

where Dxx , Dxy , Dxz , Dyy , Dyz , and Dzz pL{T q are the elements (dispersion coeffcients) of the hydrody-
namic dispersion tensor D, which is assumed to be symmetric (Dyx “ Dxy , Dzx “ Dxz , and Dzy “ Dyz ).
Thus, in the most general case, the solute mass fux component along a particular direction cannot be com-
puted solely based on the component of the concentration gradient along that one direction; three independent
components of the gradient are required. However, in certain special cases that are discussed later in this sec-
tion, fux components reduce to the simplifed form

BC
fsDSP “ ´Sw θDss , (4–31)
Bs

where s represents the direction along which the concentration gradient is evaluated and the fux component is
computed. Direction s may represent one of the coordinate directions, x, y, or z, or another direction, depend-
ing on the special case.
The “simplifed formulation” of the dispersive solute mass fow rate between adjacent cells n and m is a
discrete analog of equation 4–31 and has the form

DSP
Fn,m ˜ n,m pCm ´ Cn q .
“D (4–32)

The coeffcient D̃n,m (L3 {T ) is analogous to a “conductance” for fow of solute mass driven by a concentra-
tion difference and incorporates the effects of the dispersion coeffcients and porosities in the two model cells,
the interfacial area over which the dispersive fux occurs, and the distance between the nodes at which the cell
concentrations are calculated.
For the simplifed formulation, a dispersion conductance is defned as a coeffcient, that when multiplied
by a concentration difference, will result in the dispersive mass fux. The dispersion conductance is calculated
based on the harmonic mean of two half-cell dispersion conductances as

d˜n d˜m
D̃n,m “ , (4–33)
d˜n ` d˜m

where d˜n is the calculated half-cell dispersion conductance for cell n in the direction of cell m and d˜m is the
calculated half-cell dispersion conductance for cell m in the direction of cell n. The half-cell dispersion con-
ductance is calculated for cell n (in the direction of cell m) as
 Chapter 4. Advective and Dispersive Solute Transport 4–11

Dn,m An,m
d˜n “ , (4–34)
Ln,m

and for cell m in the direction of cell n as

Dm,n Am,n
d˜m “ . (4–35)
Lm,n

In equations 4–34 and 4–35, the effective dispersion coeffcients, Dn,m and Dm,n , are interpolated from
the principal fow-aligned dispersion components D11 , D22 , and D33 . D11 is aligned with the fow direction,
and D22 and D33 are aligned with the two orthogonal transverse directions. An,m and Am,n are equal and are
calculated as the area for fow between cells n and m. For a horizontal connection, the fow area is a function
of cell saturation. The distance between cell n and its shared face with cell m is denoted by Ln,m . Likewise,
the distance between cell m and its shared face with cell n is denoted by Lm,n .
With the simplifed formulation, the value for the effective dispersion coeffcient Dn,m in the n-m direc-
tion is calculated from D11 , D22 , and D33 . For the NPF Package (Langevin and others, 2017) the simplifed
approach interpolates an effective hydraulic conductivity value from a hydraulic conductivity ellipsoid. For
the dispersion coeffcient, however, the following simple linear equation is used to calculate an effective value
as

2 2 2
Dn,m “ νD1 D11 ` νD2 D22 ` νD3 D33 , (4–36)

where νD1 , νD2 , and νD3 defne the three components of a unit vector pointing in the n-m direction and refer-
enced in the local xL , xT 1 , xT 2 fow-aligned coordinate system for cell n. A separate calculation is made for
Dm,n , which is the effective dispersion coeffcient for model cell m in the n direction.
The simplifed formulation defned by equation 4–32 can provide an accurate estimate of the solute mass
fow if (1) all the fux expressions that need to be evaluated reduce to the simplifed form in equation 4–31, and
(2) the model grid satisfes certain geometric requirements. The MODFLOW 6 GWF Model documentation
(Langevin and others, 2017) discusses these “CVFD requirements” in detail and summarizes them as follows:

“For accurate solutions, the standard CVFD formulation requires that a line drawn between the cen-
ters of two connected cells should intersect the shared face at a right angle ... . Furthermore, the
intersection point should coincide with an appropriate mean position on the shared face (Narasimhan
and Witherspoon, 1976). ... Although this CVFD requirement is met for a simple grid of regular
polygons, equilateral triangles, and rectangles, it is violated for nested grids and may be violated for
grids with nonregular polygon-shaped cells. ... The smaller the deviation from this CVFD require-
ment, the smaller the loss of accuracy in the groundwater fow solution. In addition, the errors gener-
ally decrease as resolution increases, but they are diffcult to quantify.”

The constraints on the fux expressions and grid geometry described above imply that the simplifed formula-
tion defned by equation 4–32 can provide an accurate estimate of the solute mass fow in the following special
cases:

1. The MODFLOW grid is regular and aligned with the model coordinates, and fow is unidirectional
along one of the grid directions. In this case, every fux component calculation is performed along a
grid direction and depends only on the concentration gradient component along that same grid direction.
4–12 Documentation for the MODFLOW 6 Groundwater Transport Model

For example, if fow is uniformly in the x direction in a three-dimensional transport simulation, the y and
z components of velocity are zero, the hydrodynamic dispersion tensor is diagonal (Dxy “ Dxz “ Dyz “
0), and the fux components in equation 4–30 simplify to

BC
fxDSP “ ´Sw θDxx
Bx
DSP BC
fy “ ´Sw θDyy (4–37)
By
BC
fzDSP “ ´Sw θDzz .
Bz

Each of the three fux expressions in equation 4–37 is of the form given in equation 4–30, with s set to x,
y, or z. Furthermore, because the grid is regular, it satisfes the CVFD requirements.

2. The MODFLOW grid satisfes the CVFD requirements, and the longitudinal and relevant trans-
verse dispersivities are always equal to each other, regardless of the fow direction. Satisfaction of
the CVFD requirements is stipulated, so it remains only to discuss the ramifcations of equal longitudi-
nal and transverse dispersivities, which are defned in equation 4–20. For a three-dimensional simulation,
both transverse dispersivities are relevant, and the requirement is that αL “ αT 1 “ αT 2 for all fow direc-
tions, in which case the hydrodynamic dispersion tensor is diagonal (Dxy “ Dxz “ Dyz “ 0) and the
fux components in equation 4–30 simplify to the forms in equation 4–37 (with Dxx “ Dyy “ Dzz ). This
is accomplished by setting αLH “ αT H1 “ αT H2 and αLV “ αT V . For a two-dimensional simulation
within the (x, y) (horizontal) plane, αT 2 is irrelevant, and the requirement that αL “ αT 1 for all fow
directions is satisfed by setting αLH “ αT H1 ; the values of αLV , αT V , and αT H2 have no effect. For
a two-dimensional simulation within either the (x, z) or (y, z) (vertical) plane, αT 1 is irrelevant, and the
requirement that αL “ αT 2 for all fow directions is satisfed by setting αLH “ αT H2 and αLV “ αT V ;
the value of αT H1 has no effect. For a one-dimensional simulation in either the x, y, or z direction, nei-
ther transverse dispersivity is relevant, and the constraint on the dispersivities is satisfed by default.

The special cases listed above may apply to an entire grid or to portions of a grid. In cases in which the con-
ditions described above are not met, the XT3D option discussed below can provide a more accurate estimate
of the solute mass fow rate. Although it may be possible to implement some type of ghost-node correction,
as implemented for fow by Panday and others (2013) and Langevin and others (2017), the XT3D option pro-
vides an easier way to improve the fux calculation and is the only method implemented as an alternative to the
simple fux calculation.

XT3D Formulation
The XT3D formulation overcomes the limitations of the simplifed formulation described above by relat-
ing each component of solute mass fux to all three components of the solute concentration gradient vector
(equation 4–30). In doing so, XT3D accounts for anisotropy of the dispersion tensor and irregularities in the
model grid. The solute concentration gradient vector is estimated by spatial interpolation of concentration
values at the nodes of model cells. The implementation of the XT3D formulation for solute mass fow in the
GWT Model is mathematically analogous to its implementation for groundwater fow in the GWF Model. The
following description of the XT3D formulation is adapted and summarized from the detailed discussion in
Provost and others (2017).
 Chapter 4. Advective and Dispersive Solute Transport 4–13

The XT3D method produces an expression for the solute mass fow between two model cells, n and m, as
a function of the solute concentrations in those two cells and their surrounding cells. Conceptually, the method
is the result of three main mathematical steps listed below:

1. On each side of the interface between cells n and m, construction of an expression for the concentration-
gradient vector. The expression for the “cell n” side is a function of the concentrations in cell n and its
neighbors, and an unknown concentration at the interface. The expression for the “cell m” side is a func-
tion of the solute concentrations in cell m and its neighbors, and the unknown concentration at the inter-
face.

2. On each side of the interface between cells n and m, application of the Fickian dispersion equation (equa-
tion 4–17) and calculation of an expression for the component of the solute mass fux normal to the inter-
face in terms of the concentrations mentioned above. D can be anisotropic and different on each side of
the interface.

3. Application of the continuity principle, which requires that the solute mass fow crossing the interface
between cells n and m be the same on each side of the interface. This allows the unknown concentration
at the interface to be solved for and leads to a single expression for the solute mass fow across the inter-
face in terms of the concentrations in cells n and m and their neighbors.

Implementation of the steps enumerated above is based on the following concepts:

‚ the use of weighted averaging to incorporate gradient information from neighboring connections,

‚ dependence of weights on distances and orientations, and

‚ combining information from both sides of the interface.

Figure 4–1 shows the connections used by the XT3D method to estimate the concentration gradient at the
interface between two cells, n and its neighbor m, on an unstructured grid. (Although the grid shown in fg-
ure 4–1 consists of regular hexagons, the XT3D method is applicable to more general unstructured grids, as
well.) In this context, the connection between cells n and m is the “primary connection,” and the correspond-
ing interface is the “primary interface.” The gradient is estimated at the point at which the primary connec-
tion intersects the primary interface, which is marked with an “X.” The component of the gradient along the
primary connection is estimated simply by differencing along the primary connection. An obvious approach
would be to difference over the entire length of the primary connection; that is, between the nodes n and m.
Instead, an “unknown” concentration value is temporarily assigned to the point on the primary interface, and
differencing along the primary connection is performed between node n and the interface, and between node
m and the interface. The unknown solute concentration is eventually eliminated as a variable by enforcing
continuity of the normal fux at the primary interface.
In the simplifed formulation, the solute mass fow from cell m into cell n (equation 4–32) involves con-
centrations only at nodes n and m. In contrast, the spatial interpolation (weighted averaging) performed by
XT3D results in a solute mass fow expression that involves concentration values from all neighbors of cells n
and m:

ÿ ÿ
DSP
Fn,m “ Dn,m,pn,mq pCm ´ Cn q ` Dn,p,pn,mq pCp ´ Cn q ´ Dm,q,pn,mq pCq ´ Cm q , (4–38)
pPηn qPηm
p‰m q‰n

where the frst summation is over neighbors of cell n, excluding cell m, the second summation is over neigh-
bors of cell m, excluding cell n. The “D” coeffcients are analogous to “conductances” for fow of solute mass
4–14 Documentation for the MODFLOW 6 Groundwater Transport Model

q5
p2 q6
p3

q4
q1 p1
p4
n X m

q3
p5 p6 q2

Figure 4–1. The connections used by the XT3D method, in two dimensions, to estimate the concentration gradient at a point
(“X”) on the interface between cells n and m (the “primary” interface). A separate estimate of the concentration gradient at
point X is formulated using information from each side of the primary interface. On the “node n” side of the primary interface,
the component of the gradient along the primary connection (black line) is estimated by fnite differencing the concentrations
at node n and point X, and the component of the gradient perpendicular to the primary connection is estimated using gradient-
component information from connections between node n and its neighbors p2 , ..., p6 (orange lines). An analogous procedure
involving node m and its connections with its neighbors q2 , ..., q6 (blue lines) is used to formulate an estimate of the concentra-
tion gradient for the “node m” side of the primary interface. Modifed from Provost and others (2017).
 Chapter 4. Advective and Dispersive Solute Transport 4–15

driven by concentration differences and incorporate the effects of the dispersion coeffcients in cells n and m,
the saturation-dependent interfacial area over which the dispersive fux occurs, and other geometric informa-
tion. Formulation of the “D” coeffcients for dispersive transport is directly analogous to the formulation of
the XT3D conductance-like “C” coeffcients for groundwater fow described in Provost and others (2017).
Subscripts on the “D” coeffcients indicate the cell-cell connection from which the coeffcient derives and
the cell-cell interface to which it applies. For example, Dn,p,pn,mq is the coeffcient that derives from the con-
nection between cell n and neighboring cell p and applies to calculating the solute mass fow at the interface
between cells n and m.
The XT3D option is applicable to both regular and irregular model grids, whether the dispersion tensor
is isotropic or anisotropic. The XT3D option is on by default for transport (XT3D is not on by default for the
GWF model). The simplifed formulation can be activated by turning off the XT3D option. The XT3D option
tends to be more computationally intensive than using the simplifed formulation. Before deciding whether
to use the simplifed formulation or XT3D option for production runs, the user should consider whether the
simplifed formulation alone can provide acceptable accuracy for the particular problem being solved. Trial
runs that compare solution accuracy and run times for different formulations can be helpful in this regard.
When the concentration gradient is uniform in the vicinity of a cell interface and its neighboring connec-
tions, the XT3D estimate of solute mass fow across the interface is exact. When the concentration gradient
is nonuniform, as is typically the case in practice, and anisotropy of the dispersion tensor is not aligned with
the model-coordinate axes, XT3D handles the gradient nonuniformity by weighted averaging of gradient-
component information. Such averaging is conceptually similar to the averaging done in standard fnite dif-
ferencing on a rectangular grid. If accuracy is a concern when simulating dispersion with anisotropy that is
not aligned with the model coordinates and driven by substantially nonuniform gradients, such as the gradients
associated with strong sources and sinks of solute mass, grid refnement can be used to estimate the discretiza-
tion error.
If a cell is inactive or the head is below the cell bottom, as is possible with the Newton fow formulation,
the XT3D method for the GWT Model is not used and the dispersive fux is set to zero.
Numerical simulations of highly anisotropic fow or transport based on CVFD and fnite-element dis-
cretizations can exhibit “spurious oscillations” in the solution (Pal and Edwards, 2011). Although it can be
diffcult to distinguish spurious oscillations from legitimate variations in concentration in complex fow and
transport systems, solutions calculated using XT3D for highly anisotropic systems should be evaluated criti-
cally for evidence of unrealistic patterns in concentration or solute mass fow. For example, in a steady-state
groundwater fow and transport problem, the concentration solution should not exhibit a local maximum or
minimum within the interior of the model domain unless there is a corresponding source or sink of solute mass
at that location or solute decay or production that can explain the pattern.
 Chapter 5. Sources and Sinks of Solute Mass 5–1

Chapter 5. Sources and Sinks of Solute Mass

There are three different packages in the MODFLOW 6 GWT Model that can be used to represent sources
and sinks of solute mass. These packages include the Constant-Concentration (CNC), Source-Sink Mixing
(SSM), and the Source Mass Loading (SRC) Packages, which are described in this chapter.

Constant-Concentration (CNC) Package

The CNC Package for the MODFLOW 6 GWT Model allows a cell to be marked as a constant-
concentration cell in which the solute concentration is fxed at a user-specifed value for the duration of a
time step. Users may change the value of this fxed concentration for each time step, and so it is important
to emphasize that “constant” means that the concentration is fxed for the time step. Furthermore, the user
may change the cell status from a constant-concentration condition to an active cell or from an active cell to
a constant-concentration condition.
A constant-concentration cell is handled differently from other types of sources and sinks in that the solute
concentration for the cell is not solved. Instead, the cell is marked to be handled in a special way by the solu-
tion routine. When a cell is marked as a constant concentration cell, the balance equation is manipulated such
that

An,n “ 1, (5–1)

bn “ Cs , (5–2)

and

An,m‰n “ 0, (5–3)

where Cs is the concentration specifed by the user. With this approach, the concentration for the cell is cal-
culated (Cn “ bn {An,n ) to be the concentration specifed by the user. Adjacent active cells connected to
a constant-concentration cell n will likely have non-zero A coeffcients in their balance equations, because
of advection and dispersion terms added by the ADV and DSP Packages, respectively. Upon solution of the
system of equations, those active cells will show a non-zero mass fow to or from the connected constant-
concentration cell. These mass fows are tabulated and provided to the user as model output.
Although constant-concentration conditions are easy to implement and specify for a solute transport
model, they are not often used in practice, because they can result in unrealistically large mass fuxes. For
most applications, users are encouraged to use the SSM or SRC Packages to represent solute sources and
sinks.

Flow Model Source and Sink Mixing (SSM) Package

The Source and Sink Mixing (SSM) Package simulates the change in the solute mass caused by groundwa-
ter infows and outfows from GWF stress and advanced stress packages. Stress boundaries can be assigned a
user-specifed solute concentration so that a solute mass fux can be calculated for infow to a model cell. The
addition of water from a stress boundary may increase or decrease the concentration in the cell, depending on
5–2 Documentation for the MODFLOW 6 Groundwater Transport Model

whether the user-specifed concentration is greater than or less than the cell concentration. Flow packages that
can act as sources or sinks within the SSM Package are shown in table 5–1.

Table 5–1. List of Groundwater Flow Model Packages that can act as a solute source or sink for the MODFLOW 6 Groundwater
Transport Model.

Package Name Abbreviation Package Category

Specifed Head CHD Stress
Well WEL Stress
Recharge RCH Stress
River RIV Stress
General-Head Boundary GHB Stress
Drain DRN Stress
Evapotranspiration EVT Stress
Stream-Flow Routing SFR Advanced Stress
Lake LAK Advanced Stress
Multi-Aquifer Well MAW Advanced Stress
Unsaturated Zone Flow UZF Advanced Stress

There are two alternative ways in which the advanced stress packages (SFR, LAK, MAW, and UZF;
table 5–1) can be represented in the GWT Model. The simplest way is for the user to assign a concentra-
tion for any infow into the GWT Model domain from the advanced stress package. This simpler approach
is part of the SSM Package described here. Alternatively, solute transport can be explicitly simulated within
these advanced stress packages using a more sophisticated approach, as described in the next chapter (chap-
ter 6). When this more sophisticated approach is used, solute concentrations are calculated for the individual
advanced stress package features, and these calculated concentrations are assigned to fows into the connected
model cells.
The effect of groundwater sources and sinks, as shown in the mathematical equations for solute transport
(eqs. 2–1, 2–2, and 2–4), is expressed as

f SSM “ qs1 Cs . (5–4)

In the GWT Model, the user is allowed to assign multiple solute sources and sinks to a single model cell, and
the concentrations for these individual sources and sinks also may vary. The fourth term on the left-hand side
of equation 2–9, the net rate at which solute mass is entering or leaving cell n due to external sources and
sinks associated with stress packages, is simply the sum of all nssm sources and sinks, defned as

nssm
ÿ
M9 nSSM “ Qn,issm Cn,issm (5–5)
issm“1

where Qn,issm is the volumetric groundwater fow rate (L3 {T ) into or out of cell n by way of a stress pack-
age source and sink. The variable issm is included here to indicate that the formulation works for multiple
groundwater sources and sinks assigned to the same cell. For infow to cell n, Cn,issm is the solute concentra-
tion of the stress package source. For outfow from cell n into a stress package feature, Cn,issm is the calcu-
lated solute concentration in cell n.
 Chapter 5. Sources and Sinks of Solute Mass 5–3

The effects of the mass infow or outfow due to stress package sources and sinks are added to the system
of equations differently, depending on the sign of each Qn,issm term. For fow into model cell n (Qn,issm >0),
the right-hand-side term is updated as

bn Ð bn ´ Qn,issm Cn,issm . (5–6)

If fow is out of the model cell n (Qn,issm <0), then the SSM Package has two different options for how
to determine the concentration that is assigned to the outfow. The default option is to remove the water at the
same concentration as the cell. For this condition, the A matrix is updated as

An,n Ð An,n ` Qn,issm . (5–7)

The other option is for the user to specify a concentration for the outfow as Cn,issm . If this user-specifed con-
centration is less than the concentration of the cell, then the right-hand-side term is updated according to equa-
tion 5–6. If the user-specifed concentration is greater than the cell concentration, then water is withdrawn at
the cell concentration and equation 5–7 is used.

Source Mass Loading (SRC) Package

Solute mass infow or outfow not associated with the infow or outfow from a GWF Model stress package
can be specifed by the user through the M9 nSRC term. This term allows users to add or remove solute mass
directly to or from a model cell without an associated addition or withdrawal of water. The SRC Package is
fexible in that multiple source and sink terms can be specifed, including multiple terms for a single cell, as

nsrc
ÿ
M9 nSRC “ SRC
Mn,isrc , (5–8)
isrc“1

where nsrc is the number of source terms for cell n. Values for Mn,isrc
SRC are specifed by the user, and can be

negative to remove solute mass or positive to add solute mass. If a negative value is specifed by the user
for Mn,isrc
SRC , then it is possible that the calculated solute concentration could be also negative. This cannot be

determined before the start of the simulation, however.

The effect of the mass source term is included in the system of equations by simply updating the right-
hand side as

SRC
bn Ð bn ´ Mn,isrc . (5–9)

This update is simply the addition of the user-specifed mass infow or removal rate to the balance equation,
which is written in terms of M {T .
Use of this source and sink rate can be problematic in some instances. If M9 nSRC is specifed as negative,
then it is possible that the simulated concentrations will be negative, and there is no mechanism to prevent this
other than the user adjusting the removal rate. Likewise, large source rates can result in concentrations that
may not be reasonable. It is important that the user ensures that realistic rates are used.
 Chapter 6. Transport for Advanced Stress Packages 6–1

Chapter 6. Transport for Advanced Stress Packages

The GWT Model has two different options for representing the effects of advanced stress packages (SFR,
LAK, MAW, UZF, and MVR; table 5–1) on solute transport in an aquifer. The frst option is to treat features in
the advanced stress packages as sources and sinks with user-specifed concentrations and use the SSM Package
described in the previous chapter. With this SSM approach, the user can assign a solute concentration to fea-
tures in the advanced stress packages. If water fows into a model cell from the feature, then that water enters
the cell at the concentration specifed by the user. If water leaves the model cell, then water exits the cell at the
calculated concentration of the cell. With this approach, solute is not transported or tracked once it leaves the
GWT Model cell and enters an advanced package feature.
In the second option, solute transport can be explicitly simulated for the advanced stress packages of the
GWF Model (Langevin and others, 2017) by solving solute balance equations for each feature. For example,
if the Streamfow Routing (SFR) Package is used to simulate fow in a stream network, then the GWT Model
can simulate solute transport within the stream network itself, and transport between the stream and underlying
aquifer. The approach presented here tightly couples solute transport in the aquifer with solute transport in the
advanced package features. Tight coupling is implemented by adding additional rows to the matrix equations
to solve a solute concentration for each advanced package feature. Specialized transport packages were writ-
ten for each of the advanced stress packages and the Water Mover Package. These transport packages include
the Streamfow Transport (SFT), Lake Transport (LKT), Multi-Aquifer Well Transport (MWT), Unsaturated
Zone Transport (UZT), and Mover Transport (MVT) Packages. For the SFR Package, a feature is an individ-
ual stream reach. For the LAK Package a feature is a single lake. For the MAW Package, a feature is a single
multi-aquifer well, and for the UZF Package, a feature is a single unsaturated zone fow cell.
The advanced package transport routines formulate a separate control-volume fnite-difference equation
for each feature in a way that is similar to what is done for a single model cell. The advanced package trans-
port routines follow the same general pattern, which involves solving a form of the following equation,

M9 nstorage ` M9 nadvection ` M9 nto´mover ` M9 nf rom´mover ` M9 nsinks{sources “ 0, (6–1)

where the n subscript indicates that the equation is for feature n, M9 nstorage is the rate of change of solute mass
in feature n per time (a positive value for M9 nstorage indicates that solute mass in feature n is decreasing),
M9 nadvection is the sum of all advective infow and outfow, M9 nto´mover is the rate of mass transfer to the water
sinks{sources
mover, M9 nf rom´mover is the rate of solute mass added to the feature from the water mover, and M9 n
is the rate of solute mass added to or removed from the feature due to sources and sinks, which are defned
separately for each advanced transport package (SFT, LKT, MWT, and UZT).
In the version of the GWT Model described here, the advanced package transport routines do not rep-
resent sorption, decay, or a dispersive fux between adjacent features or a feature and adjacent GWT Model
cells. Although these terms could be added to these advanced transport packages in the future, the focus of the
present implementation is on advective transport of a conservative dissolved constituent.
The terms in equation 6–1 are described in the remainder of this chapter. The storage and advection terms
are solved using generalized routines that work with all of the advanced packages, including SFT, LKT, MWT,
and UZT. The source and sink term in equation 6–1 has different implementations depending on the advanced
package. For example, SFT and LKT can both be affected by rainfall, which may have an associated solute
concentration, whereas rainfall is not a source term for MWT. Accordingly, the sink and source term, which
may contain multiple types for an advanced package, is described in subsequent sections for those packages.
The mover terms in equation 6–1 are described in the “Water Mover Terms and the Mover Transport (MVT)
Package” section of this chapter.
6–2 Documentation for the MODFLOW 6 Groundwater Transport Model

Storage Term
The storage term in equation 6–1 is expressed as

1 `
M9 storage “ ´C t`∆t Vwt`∆t ` C t Vwt ,
˘
(6–2)
∆t
where the solute concentration (C) and volume (Vw ) terms represent the concentration and water volume of
feature n at the start (t) and end (t ` ∆t) of the time step.
The equation for aqueous solute storage is added to the system of equations by updating the diagonal posi-
tion of the A matrix for row n, corresponding to a specifc feature within an advanced package, with the coef-
fcient of the C t`∆t term in equation 6–2 as

Vwt`∆t
An,n Ð An,n ´ , (6–3)
∆t
C t Vwt
and by moving the remaining ∆t term to the right-hand side to update bn as

C t Vwt
bn Ð bn ´ . (6–4)
∆t
Implementation of this generic storage term allows the solute concentration to be calculated within each fea-
ture based on changes in water volume and the solute balance for each feature.

Advection Term
The advection term in equation 6–1 includes fows to and from adjacent features in an advanced package
as well as fows to and from connected GWT Model cells. Following the weighted advective fux equation 4–
12, this advection term is expressed as

ÿ
M9 advection “ ωn,m Qn,m Cn ` p1 ´ ωn,m q Qn,m Cm , (6–5)
mPηn

where m is an adjacent advanced package feature or connected GWT Model cell, ηn is a list of all connected
features and GWT model cells, ωn,m is the upstream weighting factor determined using equation 4–5, and
Qn,m is the volumetric fow rate between n and m defned as positive into n.
As shown in equations 4–13 and 4–14, the matrix coeffcient in the diagonal position for row n is updated
by adding the term ωn,m Qn,m as

An,n Ð An,n ` ωn,m Qn,m . (6–6)

The matrix coeffcient in row n that corresponds to the connection between cell n and neighboring cell m is
updated by adding the term p1 ´ ωn,m q Qn,m as

An,m Ð An,m ` p1 ´ ωn,m q Qn,m . (6–7)

 Chapter 6. Transport for Advanced Stress Packages 6–3

Water Mover Terms and the Mover Transport (MVT) Package

The Water Mover (MVR) Package of the MODFLOW 6 GWF Model can be used to dynamically transfer
water from a provider package to a receiver package, based on rules described by (Langevin and others, 2017).
Providers can include regular stress package features, such as a drain, a general-head boundary, or a river cell,
for example. Providers may also include a feature in an advanced stress package, such as a stream, lake, multi-
aquifer well, or an unsaturated zone cell. Receivers are restricted to a feature in an advanced stress package;
there is no way to transfer water to one of the regular stress packages, because they do not solve a continuity
equation. The M9 nto´mover and M9 nf rom´mover terms in equation 6–1 represent the solute mass fux removed
and added by the mover, respectively. The present mover transport implementation does not include dispersive
transport as part of this mass transfer; only advective transport is represented.
When a feature in a regular or advanced stress package acts as a provider to the water mover, solute leaves
the feature at a rate equal to the product of the water mover fow rate and the solute concentration of the fea-
ture:

M9 nto´mover “ Qnto´mover Cn , (6–8)

where Qto´mover
n is the volumetric rate of water transferred to the water mover (which is negative in sign as
water is being removed from feature n), and Cn is the simulated concentration of the feature in the stress pack-
age. The effect of the transferred water is included in the system of equations by updating the diagonal posi-
tion of the A matrix for row n as

An,n Ð An,n ` Qto´mover

n . (6–9)

When a feature in an advanced stress package acts as a receiver, then the water that it receives can have an
associated solute concentration. The Water Mover Transport (MVT) Package was implemented for the MOD-
FLOW 6 GWT Model to facilitate the transfer of solute from providers to receivers. Because the GWF Model
MVR Package is generalized and allows multiple providers to transfer water to a single receiver, the GWT
Model MVT Package accumulates these solute mass fuxes into a M9 nf rom´mover term (which is positive in
sign) for each feature of an advanced package as the system of equations is formulated. The advanced package
then adds the accumulated solute mass fux as a source of mass from the mover by updating the right-hand-
side vector as

bn Ð bn ´ M9 nf rom´mover . (6–10)

Thus, the M9 nto´mover terms are handled in a fully implicit manner and added to the left-hand side, whereas the
M9 nf rom´mover terms are added to the right-hand side. A fully implicit approach could be added in the future
for the M9 nf rom´mover terms as a way to improve convergence of the numerical solution; however, this addition
of an implicit approach would require adding matrix-level connections for all provider and receiver combina-
tions.

Streamfow Transport (SFT) Package

The Streamfow Transport (SFT) Package simulates solute concentrations in each stream reach based on a
numerical solution of the solute transport balance equation (eq. 2–12). Because the corresponding SFR Pack-
age simulates one-dimensional vertically averaged fow within a stream channel, the SFT Package simulates a
6–4 Documentation for the MODFLOW 6 Groundwater Transport Model

single solute concentration for each reach. This approach assumes that all of the water entering a reach instan-
taneously mixes with water in the stream reach. Thus, there is no way to represent stratifcation or variations
of solute concentration with depth using the current SFT Package implementation.

The generalized balance equation 6–1 for an advanced package feature, which is a reach for SFT, includes
sinks{sources
a source and sink term, M9 n . Source and sink terms for the SFT Package include rainfall, evapora-
tion, external infow, and external outfow. These source and sink terms are included in the balance equation
by updating the diagonal position of the A matrix and the right-hand-side vector element for row n according
to the equations in table 6–1.

Table 6–1. Source and sink equations for Streamfow Transport Package of the MODFLOW 6 Groundwater Transport Model.
[The An,n Ð column represents the term that is added to the diagonal position of the A matrix; bn Ð column represents the term that is added to the right-hand-side vector; Cs is the user spec-
ifed concentration for this source or sink term; Qs is the volumetric fow rate provided by the corresponding Groundwater Flow Model stress package for this source or sink term; Cn is the solute
concentration of the stream reach; ω is a weighting factor used to shift between the left-hand-side and right-hand-side implementations based on solute concentration]

sinks{sources
Source or Sink M9 n An,n Ð bn Ð Note
Rainfall (source) Cs Qs zero ´Cs Qs - $
&1 if Cn ă Cs
Evaporation (sink) ωQs Cn ` p1 ´ ωq Qs Cs ωQs ´ p1 ´ ωq Qs Cs ω“
%0 if Cn ě Cs
External infow (source) Cs Qs zero ´Cs Qs -
External outfow (sink) Cn Qs Qs zero -

Lake Transport (LKT) Package

The Lake Transport (LKT) Package simulates solute concentrations in each lake based on a numerical
solution of the solute transport balance equation (eq. 2–12). Because the corresponding LAK Package simu-
lates a single stage for the entire lake and does not represent fow within the lake, the LKT Package simulates a
single concentration for each lake. This approach assumes that all of the water entering a lake instantaneously
and thoroughly mixes with lake water. Thus, there is no way to represent stratifcation or variations of solute
concentration with depth using the current LKT Package implementation.

The generalized balance equation 6–1 for an advanced package feature, which is a lake for LKT, includes
sinks{sources
a source and sink term, M9 n . Source and sink terms for the LKT Package include rainfall, evap-
oration, external infow, external outfow, and withdrawal. These source and sink terms are included in the
balance equation by updating the diagonal position of the A matrix and the right-hand-side vector element for
row n according to the equations in table 6–2.
 Chapter 6. Transport for Advanced Stress Packages 6–5

Table 6–2. Source and sink equations for Lake Transport Package of the MODFLOW 6 Groundwater Transport Model.
[The An,n Ð column represents the term that is added to the diagonal position of the A matrix; bn Ð column represents the term that is added to the right-hand-side vector; Cs is the user spec-
ifed concentration for this source or sink term; Qs is the volumetric fow rate provided by the corresponding Groundwater Flow Model stress package for this source or sink term; Cn is the solute

concentration of the lake; ω is a weighting factor used to shift between the left-hand-side and right-hand-side implementations based on solute concentration]

sinks{sources
Source or Sink M9 n An,n Ð bn Ð Note
Rainfall (source) Cs Qs zero ´Cs Qs - $
&1 if C ă C
n s
Evaporation (sink) ωQs Cn ` p1 ´ ωq Qs Cs ωQs ´ p1 ´ ωq Qs Cs ω“
%0 if Cn ě Cs
External infow (source) Cs Qs zero ´Cs Qs -
External outfow (sink) Cn Qs Qs zero -
Withdrawal (sink) Cn Qs Qs zero -

Multi-Aquifer Well Transport (MWT) Package

The Multi-Aquifer Well Transport (MWT) Package simulates solute concentrations in each well based
on a numerical solution of the solute transport balance equation (eq. 2–12). The MWT Package simulates a
single concentration for each well. This approach assumes that all of the water entering a multi-aquifer well is
instantaneously and thoroughly mixed with well water. There is no option at present to represent stratifcation
or variations of solute concentration with depth in a well.
The generalized balance equation 6–1 for an advanced package feature, which is a multi-aquifer well for
sinks{sources
MWT, includes a source and sink term, M9 n . Source and sink terms for the MWT Package include
withdrawal and fowing well rate (rate of free fowing discharge due to artesian conditions). These source
and sink terms are included in the balance equation by updating the diagonal position of the A matrix and the
right-hand-side vector element for row n according to the equations in table 6–3.

Table 6–3. Source and sink equations for Multi-Aquifer Well Transport Package of the MODFLOW 6 Groundwater Transport
Model.
[The An,n Ð column represents the term that is added to the diagonal position of the A matrix; bn Ð column represents the term that is added to the right-hand-side vector; Cs is the user spec-
ifed concentration for this source or sink term; Qs is the volumetric fow rate provided by the corresponding Groundwater Flow Model stress package for this source or sink term; Cn is the solute
concentration of the multi-aquifer well; ω is a weighting factor used to shift between the left-hand-side and right-hand-side implementations based on the sign of the well fow rate]

sinks{sources
Source or Sink M9 n An,n Ð bn Ð Note$
&1 if Q ă 0
s
Well rate (source or sink) ωQs Cn ` p1 ´ ωq Qs Cs ωQs ´ p1 ´ ωq Qs Cs ω“
%0 if Qs ě 0
Flowing well rate (sink) Cn Qs Qs zero -

Unsaturated Zone Transport (UZT) Package

The Unsaturated Zone Transport (UZT) Package simulates solute concentrations in each UZT cell based
on a numerical solution of the solute transport balance equation (eq. 2–12). Although the corresponding UZF
Package simulates individual wetting fronts, the UZT Package calculates an average concentration for the
entire UZT cell. Individual wetting front concentrations are not calculated or tracked, although this feature
could be implemented in the future.
6–6 Documentation for the MODFLOW 6 Groundwater Transport Model

The generalized balance equation 6–1 for an advanced package feature, which is an unsaturated zone cell
sinks{sources
for UZT, includes a source and sink term, M9 n . Source and sink terms for the UZT Package include
infltration, rejected infltration, and unsaturated zone evapotranspiration. These source and sink terms are
included in the balance equation by updating the diagonal position of the A matrix and the right-hand-side
vector element for row n according to the equations in table 6–4.

Table 6–4. Source and sink equations for Unsaturated Zone Transport Package of the MODFLOW 6 Groundwater Transport
Model.
[The An,n Ð column represents the term that is added to the diagonal position of the A matrix; bn Ð column represents the term that is added to the right-hand-side vector; Cs is the user spec-
ifed concentration for this source or sink term; Qs is the volumetric fow rate provided by the corresponding Groundwater Flow Model stress package for this source or sink term; Cn is the solute
concentration of the unsaturated zone feature; ω is a weighting factor used to shift between the left-hand-side and right-hand-side implementations based on solute concentration]

sinks{sources
Source or Sink M9 n An,n Ð bn Ð Note
Infltration (source) Cs Qs zero ´Cs Qs -
Rejected infltration (sink) Cs Qs zero ´Cs Qs - $
&1 if C ă C
n s
Evapotranspiration (sink) ωQs Cn ` p1 ´ ωq Qs Cs ωQs ´ p1 ´ ωq Qs Cs ω“
%0 if Cn ě Cs
 Chapter 7. Immobile Domain Storage and Transfer 7–1

Chapter 7. Immobile Domain Storage and Transfer

Models of solute transport in heterogeneous groundwater systems often rely on conceptual models that
involve multiple domains. For example, if a system can be conceptualized as having a fast (mobile) domain
and a slow (immobile) domain, then advection and dispersion can be represented in the mobile domain and
assumed to be negligible in the immobile domain. Under this conceptual model, the exchange between the
mobile and immobile domains can be represented as a diffusive process with a simple linear exchange coef-
fcient. Solute moves through the system in the mobile domain and can diffuse into or out of the immobile
domain based on concentration differences. If the concentration in the mobile domain is higher than the con-
centration in the immobile domain, then solute mass will diffuse into the immobile domain. Alternatively, if
the concentration in the mobile domain is less than the concentration in the immobile domain, then solute mass
will diffuse into the mobile domain. This conceptual model has been implemented in many solute transport
modeling programs, including MOC3D (Goode, 1999) and MT3DMS (Zheng and Wang, 1999).
Conceptual models for mobile and immobile domains are shown in fgure 7–1. In fgure 7–1A, an ide-
alized sand aquifer is shown to contain evenly distributed and isolated clay and silt nodules. If the size of
these nodules is similar or greater than the size of a model cell, then these low permeability features should
be explicitly represented in the fow model as cells with reduced hydraulic conductivity compared to aquifer
cells. If these nodules are much smaller than a model cell and are evenly distributed, then their effects on
solute transport can be represented using immobile domain concepts. Likewise, fractured aquifers (fg. 7–1B)
that can be approximated as an equivalent porous medium for groundwater fow may beneft from immobile
domain concepts for solute transport. The fractures represent the fast moving mobile domain, and the less per-
meable rock matrix acts as an immobile domain that is capable of storing and releasing solute mass.
With the immobile domain conceptual model, groundwater does not fow through these nodules, and thus,
their pore water is immobile. This conceptual model dictates that mobile domain porosity be calculated as the
volume of voids in the mobile domain divided by the total volume. Likewise, the immobile domain porosity
θim is defned as the volume of immobile domain pores divided by the total volume. Concentration within
the mobile domain is tracked separately from the concentration in the immobile domain. For a model cell,
an average concentration is calculated for the mobile domain, and an average concentration is calculated for
each immobile domain. The fux of solute mass between the mobile domain and immobile domain im is rep-
resented as a diffusive process that can be expressed in the form

f IST “ ´ζim Sw pC ´ Cim q , (7–1)

where Sw is the mobile-domain saturation, and the immobile domain is assumed to be fully saturated. The
transfer rate is calculated based on the concentration difference between the mobile and immobile domains and
a mass transfer coeffcient, ζim (1{T ), and the area available for mass transfer is assumed to be proportional
to Sw . Under idealized conditions, values for ζim can be calculated based on the surface area of the immobile
domain in contact with the mobile domain at full saturation and other geometric properties of the immobile
domain, but in practice, ζim and other immobile domain properties are estimated through model calibration.
The Immobile domain Storage and Transfer (IST) Package in MODFLOW 6 simulates the effects of mass
transfer between the mobile domain and an immobile domain. In the present implementation, there can be as
many immobile domains as necessary, with each immobile domain being represented by a separate instance of
an IST Package. For example, one immobile domain package can represent the clay nodules in fgure 7–1A,
and another immobile domain package can be used to represent the silt nodules. Each immobile domain can
be assigned different properties and exchange coeffcients, and these values can vary throughout the model
grid. In this case, a separate concentration is also calculated for each immobile domain. Immobile domains do
not interact with other immobile domains; an immobile domain can only interact with the mobile domain.
7–2 Documentation for the MODFLOW 6 Groundwater Transport Model

A. B.

silt

clay

rock

sand
fracture

Figure 7–1. Two conceptualizations for mobile and immobile domains: (A) sand aquifer containing clay and silt nodules, and
(B) fractured aquifer.

The immobile domain terms in equations 2–1, 2–2, and 2–4 contribute to the mass balance equation for
solute in the mobile domain. In order to represent the dual domain mass transfer process, a separate governing
equation must be written for the immobile domain (Zheng and Wang, 1999; Zheng and Bennett, 2002). If writ-
ten to include the effects of sorption and decay in the immobile domain, as well as the transfer with the mobile
domain, the governing equation for the immobile domain is

BCim BC im
θim ` fim ρb “ ´λ1,im θim Cim ´ λ2,im fim ρb C im
Bt Bt
´γ1,im θim ´ γ2,im fim ρb ` ζim Sw pC ´ Cim q , (7–2)

where θim is the volume of the immobile pores per volume of aquifer, fim is the fraction of aquifer solid mate-
rial available for sorptive exchange with the immobile domain under fully saturated conditions, C im is the
sorbed concentration of the immobile domain, expressed as the mass of the sorbed chemical per mass of solid,
λ1,im is the frst-order reaction rate coeffcient for the liquid phase of the immobile domain (1{T ), λ2,im is the
frst-order reaction rate coeffcient for the sorbed phase of the immobile domain (1{T ), γ1,im is the zero-order
reaction rate coeffcient for the liquid phase of the immobile domain (M L´3 T ´1 ), and γ2,im is the zero-order
reaction rate coeffcient for the sorbed phase of the immobile domain (M M ´1 T ´1 ). As noted by Zheng and
Bennett (2002) fim , can be approximated as fim “ θim {θ. This approximation is the default approach imple-
mented in MODFLOW 6 and can be shown to be restricted to cases where the immobile domain porosity is
the same as the mobile domain porosity and half of the domain is mobile and the other half is immobile, for
example. The MODFLOW 6 program may be modifed in the future so that users can specify values for fim
based on the fractions and porosities of the domains.
As noted by Ma and Zheng (2011) different forms of ζim have been reported in the literature. Ma and
Zheng (2011) defne classic and alternate forms, which differ by a factor of θim . With the approach devel-
oped here, based on equations 7–1 and 7–2, the defnition of ζim in this report corresponds to the classic form,
which is the form implemented in the various versions of MT3D.
Care should be taken when using the IST Package for simulations in which cells can become dry. When a
cell becomes dry, the cell is removed from the solution and no transport calculations are made. Thus, if there
is solute mass in the immobile domain, the mass may become trapped, and unable to reenter the groundwater
 Chapter 7. Immobile Domain Storage and Transfer 7–3

fow system. Any mass trapped in the immobile domain for a dry cell also will not be subjected to frst-order
or zero-order decay or production, though this may be supported in future versions. These issues do not occur
if the fow model uses the Newton formulation. With the Newton formulation, model cells are not removed
from the solution; instead, these cells remain active and subject to decay and production processes.
Following the solution approach outlined in Zheng and Bennett (2002), a discretized form of equation 7–1
can be written for a model cell as

t`∆t
M9 nIST “ ´Vcell Swt`∆t ζim C t`∆t ` Vcell Swt`∆t ζim Cim , (7–3)

where a positive value indicates solute transfer from the immobile domain into the mobile domain. The sec-
ond term on the right side of equation 7–3, which contains the immobile domain concentration, is recast as
t`∆t
a function of the mobile domain concentration. In order to express Cim as a function of mobile domain
solute concentrations and other known terms, an implicit fnite-difference approximation can be applied to
equation 7–2, the balance equation for the immobile domain. The implicit fnite-difference equation for the
immobile domain is

θim Vcell t`∆t θim Vcell t fim ρb Vcell t`∆t fim ρb Vcell t
Cim ´ Cim ` Kd Cim ´ Kd Cim
∆t ∆t ∆t ∆t
“
t`∆t t`∆t
´λ1,im θim Vcell Cim ´ λ2,im fim Vcell ρb Kd Cim
´γ1,im θim Vcell ´ γ2,im fim ρb Vcell
t`∆t
`Vcell Sw ζim C t`∆t ´ Vcell Swt`∆t ζim Cim
t`∆t
. (7–4)

t`∆t
From this equation, an expression for Cim as a function of C t`∆t and other known terms can be written as

θim Vcell fim ρb Vcell

t`∆t Vcell Swt`∆t ζim t`∆t ∆t ` ∆t Kd t γ1,im θim Vcell γ2,im fim ρb Vcell
Cim “ C ` Cim ´ ´ ,
F F F F
(7–5)

where F is defned as

θim Vcell fim ρb Vcell

F “ ` Kd ` λ1,im θim Vcell ` λ2,im fim Vcell ρb Kd ` Vcell Swt`∆t ζim . (7–6)
∆t ∆t
t`∆t
This expression for Cim in equation 7–5 can be substituted into the last term on the right side of equation 7–
3. The diagonal position of the A matrix and the right-hand-side vector element for row n are then updated to
include the transfer with the immobile domain as

` ˘2
Vcell Swt`∆t ζim
An,n Ð An,n ´ Vcell Swt`∆t ζim ` , (7–7)
F
and
7–4 Documentation for the MODFLOW 6 Groundwater Transport Model

Vcell Swt`∆t ζim

„ˆ ˙ ȷ
θim Vcell fim ρb Vcell t
bn Ð bn ´ ` Kd Cim ´ γ1,im θim Vcell ´ γ2,im fim ρb Vcell .
F ∆t ∆t
(7–8)
 References Cited R–1

References Cited

Bear, J., 1972, Dynamics of fuids in porous media: New York, New York, Dover Publications, Inc., 757 p.
Bedekar, V., Morway, E.D., Langevin, C.D., and Tonkin, M.J., 2016, MT3D-USGS version 1: A U.S. Geolog-
ical Survey release of MT3DMS updated with new and expanded transport capabilities for use with MOD-
FLOW: U.S. Geological Survey Techniques and Methods, book 6, chap. A53, 69 p., [Link]
tm6a53, [Link]
Bird, R., Stewart, W., and Lightfoot, E., 2006, Transport Phenomena: New York, New York, John Wiley &
Sons, 432 p.
Burnett, R., and Frind, E., 1987, Simulation of contaminant transport in three dimensions, 2, dimensionality
effects: Water Resources Research, v. 23, p. 695–705.
Clement, T.P., 1997, RT3D (Version 1.0) A modular computer code for simulating reactive multi-species trans-
port in 3-dimensional groundwater systems: Prepared for the U.S. Department of Energy Under Contract
DE-AC06-76RLO 1830.
Forsyth, P.A., Unger, A. J.A., and Sudicky, E.A., 1998, Nonlinear iteration methods for nonequilibrium multi-
phase subsurface fow: Advances in Water Resources, v. 21, p. 433–449.
Goode, D.J., 1990, Governing equations and model approximation errors associated with the effects of fuid-
storage transients on solute transport in aquifers: U.S. Geological Survey Water-Resources Investigations
Report 90–4156, 20 p.
Goode, D.J., 1996, Direct simulation of groundwater age: Water Resources Research, v. 32, no. 2, p. 289–296.
Goode, D.J., 1999, Age, double porosity, and simple reactions modifcations for the MOC3D ground-water
transport model: U.S. Geological Survey Water-Resources Investigations Report 99–4041, 34 p.
Guo, W., and Langevin, C.D., 2002, User’s Guide to SEAWAT: A Computer Program for Simulation of
Three-Dimensional Variable-Density Ground-Water Flow: U.S. Geological Survey Techniques of Water-
Resources Investigations book 6, chap. A7, 77 p.
Harbaugh, A.W., 2005, MODFLOW-2005, the U.S. Geological Survey modular ground-water model—the
Ground-Water Flow Process: U.S. Geological Survey Techniques and Methods, book 6, chap. A16, vari-
ously paged, accessed June 27, 2017, at [Link]
Hornberger, G., and Konikow, L.F., 2006, Use of the Multi-Node Well (MNW) package when simulating
solute transport with the MODFLOW ground-water transport process: U.S. Geological Survey Techniques
and Methods 6–A15, 34 p.
Hornberger, G., Konikow, L.F., and Harte, P., 2002, Simulating solute transport across Horizontal-Flow Bar-
riers using the MODFLOW Ground-Water Transport Process: U.S. Geological Survey Open-File Report
02–52, 28 p.
Hughes, J.D., Langevin, C.D., and Banta, E.R., 2017, Documentation for the MODFLOW 6 framework: U.S.
Geological Survey Techniques and Methods, book 6, chap. A57, 36 p., accessed August 4, 2017, at https:
//[Link]/10.3133/tm6A57.
Keating, E., and Zyvoloski, G., 2009, A stable and effcient numerical algorithm for unconfned aquifer analy-
sis: Ground Water, v. 47, no. 4, p. 569–579, accessed June 27, 2017, at [Link]
2009.00555.x.
Konikow, L.F., 2010, The secret to successful solute-transport modeling: Ground Water, v. 49, no. 2, p. 144–
159, [Link]
Konikow, L.F., and Bredehoeft, J., 1978, Computer model of two-dimensional solute transport and dispersion
in ground water: U.S. Geological Survey Techniquest of Water-Resources Investigations, book 7, chap. C2,
90 p.
R–2 Documentation for the MODFLOW 6 Groundwater Transport Model

Konikow, L.F., and Grove, D., 1977, Derivation of equations describing solute transport in ground water: U.S.
Geological Survey Techniquest of Water-Resources Investigations Report 77-19, [Revised 1984], 30 p.
Konikow, L.F., Goode, D.J., and Hornberger, G.Z., 1996, A three-dimensional method-of-characteristics
solute-transport model (MOC3D): U.S. Geological Survey Water-Resources Investigations Report 96–4267,
87 p., accessed June 27, 2017, at [Link]
Langevin, C.D., Shoemaker, W.B., and Guo, W., 2003, MODFLOW-2000 the U.S. Geological Survey Mod-
ular Ground-Water Model–Documentation of the SEAWAT-2000 Version with the Variable-Density Flow
Process (VDF) and the Integrated MT3DMS Transport Process (IMT): U.S. Geological Survey Open-File
Report 03-426, 43 p., accessed July 25, 2019, at [Link]
Langevin, C.D., Thorne Jr, D.T., Dausman, A.M., Sukop, M.C., and Guo, W., 2008, SEAWAT Version 4—
A computer program for simulation of multi-species solute and heat transport: U.S. Geological Survey
Techniques and Methods, book 6, chap. A22, 39 p., accessed June 27, 2017, at [Link]
publication/tm6A22.
Langevin, C.D., Hughes, J.D., Provost, A.M., Banta, E.R., Niswonger, R.G., and Panday, S., 2017, Documen-
tation for the MODFLOW 6 Groundwater Flow (GWF) Model: U.S. Geological Survey Techniques and
Methods, book 6, chap. A55, 197 p., accessed August 4, 2017, at [Link]
Langevin, C.D., Panday, S., and Provost, A.M., 2020, Hydraulic-head formulation for density-dependent fow
and transport: Groundwater, v. 58, no. 3, p. 349–362.
Lichtner, P.C., Kelkar, S., and Robinson, B., 2002, New form of dispersion tensor for axisymmetric porous
media with implementation in particle tracking: Water Resources Research, v. 38, no. 8, p. 21–1–21–16,
[Link]
Ma, R., and Zheng, C., 2011, Not all mass transfer rate coeffcients are created equal: Groundwater, v. 49,
no. 6, p. 772–774, [Link]
McDonald, M.G., Harbaugh, A.W., Orr, B.R., and Ackerman, D.J., 1992, A method of converting no-fow
cells to variable-head cells for the U.S. Geological Survey modular fnite-difference ground-water fow
model: U.S. Geological Survey Open-File Report 91–536, 99 p., accessed June 27, 2017, at [Link]
[Link]/publication/ofr91536.
Merritt, M.L., and Konikow, L.F., 2000, Documentation of a computer program to simulate lake-aquifer inter-
action using the MODFLOW ground-water fow model and the MOC3D solute-transport model: U.S.
Geological Survey Water-Resources Investigations Report 00–4167, 146 p., accessed June 27, 2017, at
[Link]
Morway, E.D., Langevin, C.D., and Hughes, J.D., 2021, Use of the MODFLOW 6 Water Mover Package to
represent natural and managed hydrologic connections: Groundwater, v. 59, no. 6, p. 913–924.
Narasimhan, T.N., and Witherspoon, P.A., 1976, An integrated fnite difference method for analyzing fuid
fow in porous media: Water Resources Research, v. 12, no. 1, p. 57–64, accessed June 27, 2017, at https:
//[Link]/10.1029/WR012i001p00057.
Niswonger, R.G., Panday, S., and Ibaraki, M., 2011, MODFLOW-NWT, A Newton formulation for
MODFLOW-2005: U.S. Geological Survey Techniques and Methods, book 6, chap. A37, 44 p., accessed
June 27, 2017, at [Link]
Painter, S., Başağaoğlu, H., and Liu, A., 2008, Robust representation of dry cells in single-layer MODFLOW
models: Ground Water, v. 46, no. 6, p. 873–881, accessed June 27, 2017, at [Link]
1745-6584.2008.00483.x.
Pal, M., and Edwards, M.G., 2011, Non-linear fux-splitting schemes with imposed discrete maximum princi-
ple for elliptic equations with highly anisotropic coeffcients: International Journal for Numerical Methods
in Fluids, v. 66, no. 3, p. 299–323, accessed June 27, 2017, at [Link]
 References Cited R–3

Panday, S., 2020, USG-Transport Version 1.5.0: The Block-Centered Transport Process for MODFLOW-
USG, GSI Environmental.
Panday, S., Langevin, C.D., Niswonger, R.G., Ibaraki, M., and Hughes, J.D., 2013, MODFLOW-USG ver-
sion 1—An unstructured grid version of MODFLOW for simulating groundwater fow and tightly coupled
processes using a control volume fnite-difference formulation: U.S. Geological Survey Techniques and
Methods, book 6, chap. A45, 66 p., accessed June 27, 2017, at [Link]
Panday, S., Bedekar, V., and Langevin, C.D., 2018, Impact of local groundwater fow model errors on transport
and a practical solution for the issue: Groundwater, v. 56, no. 4, p. 667–672, [Link]
12627, [Link]
Pollock, D.W., 2016, User guide for MODPATH Version 7—A particle-tracking model for MODFLOW: U.S.
Geological Survey Open-File Report 2016–1086, 35 p., accessed June 27, 2017, at [Link]
ofr20161086.
Prommer, H., Barry, D., and Zheng, C., 2003, Modfow/mt3dms-based reactive multi-component transport
modeling: Groundwater, v. 41, no. 2, p. 247–257.
Provost, A.M., Langevin, C.D., and Hughes, J.D., 2017, Documentation for the “XT3D” Option in the Node
Property Flow (NPF) Package of MODFLOW 6: U.S. Geological Survey Techniques and Methods, book 6,
chap. A56, 46 p., accessed August 4, 2017, at [Link]
Prudic, D.E., Konikow, L.F., and Banta, E.R., 2004, A New Streamfow-Routing (SFR1) Package to simulate
stream-aquifer interaction with MODFLOW-2000: U.S. Geological Survey Open-File Report 2004–1042,
104 p., accessed June 27, 2017, at [Link]
Scheidegger, A.E., 1961, General theory of dispersion in porous media: Journal of Geophysical Research
(1896-1977), v. 66, no. 10, p. 3273–3278, [Link]
Toride, N., Leij, F.J., and van Genuchten, M.T., 1993, A comprehensive set of analytical solutions for nonequi-
librium solute transport with frst-order decay and zero-order production: Water Resources Research, v. 29,
no. 7, p. 2167–2182, [Link]
Voss, C.I., and Provost, A.M., 2010, SUTRA—a model for saturated-unsaturated, variable-density ground-
water fow with solute or energy transport: U.S. Geological Survey Water-Resources Investigations Report
02–4231, 291 p.
Winston, R.B., Konikow, L.F., and Hornberger, G.Z., 2018, Volume-weighted particle-tracking method for
solute-transport modeling; Implementation in MODFLOW–GWT: U.S. Geological Survey Techniques and
Methods, book 6, chap. A58, 44 p., [Link]
Zheng, C., 1990, MT3D, A modular three-dimensional transport model for simulation of advection, dispersion
and chemical reactions of contaminants in groundwater systems: Report to the U.S. Environmental Protec-
tion Agency, 170 p.
Zheng, C., 2010, MT3DMS v5.3, Supplemental User’s Guide: Technical Report Prepared for the U.S. Army
Corps of Engineers, 51 p., [Link] .
Zheng, C., and Bennett, G.D., 2002, Applied contaminant transport modeling, 2nd edition: New York, New
York, Wiley-Interscience, 560 p.
Zheng, C., and Wang, P.P., 1999, MT3DMS—A modular three-dimensional multi-species transport model for
simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; Docu-
mentation and user’s guide: Contract report SERDP–99–1: Vicksburg, Miss., U.S. Army Engineer Research
and Development Center, 169 p.
Zheng, C., Hill, M.C., and Hsieh, P.A., 2001, MODFLOW-2000, the U.S. Geological Survey Modular
Ground-Water Model—User guide to the LMT6 package, the linkage with MT3DMS for multi-species
R–4 Documentation for the MODFLOW 6 Groundwater Transport Model

mass transport modeling: U.S. Geological Survey Open-File Report 01–82, 43 p., accessed June 27, 2017,
at [Link]
Publishing support provided by the U.S. Geological Survey
Science Publishing Network, Reston Publishing Service Center

For information concerning this publication, please contact:

Integrated Modeling and Prediction Division

U.S. Geological Survey
Mail Stop 411
12201 Sunrise Valley Drive
Reston, VA 20192
(703) 648–5001
[Link]
Langevin and others—Documentation for the MODFLOW 6 Groundwater Transport Model—Techniques and Methods 6–A61

[Link]
ISSN 2328-7055 (online)