MT3DMS Groundwater Contaminant Transport Simulation

Groundwater contamination from industrial solvents, agricultural nitrates, and legacy landfill leachate represents one of the most persistent environmental challenges facing hydrogeologists and remediation engineers. MT3DMS (Modular 3-D Multi-Species Transport Model) remains the industry-standard simulator for quantifying how dissolved contaminants migrate, disperse, and react within saturated aquifer systems. When coupled with MODFLOW's flow solutions, MT3DMS delivers the predictive accuracy required for regulatory compliance, pump-and-treat design, and natural attenuation monitoring.

Architecture and MODFLOW Coupling

MODFLOW and MT3DMS coupled workflow for groundwater contaminant transport simulation

MT3DMS operates as a post-processor to MODFLOW, reading the head and cell-by-cell flow budget files generated by a completed flow simulation. This decoupled architecture means the transport solver inherits a fully converged velocity field without re-solving the flow equations at each transport time step—a critical efficiency advantage for long-duration plume simulations spanning decades.

The coupling workflow follows three stages:

Flow model execution: Run MODFLOW (or MODFLOW 6 via the FloPy interface) to generate .hds (heads) and .cbc (cell-by-cell budget) files.
MT3DMS configuration: Define the Basic Transport (BTN) package specifying species, initial concentrations, and output control; configure the Advection (ADV), Dispersion (DSP), Source/Sink Mixing (SSM), and Reaction (RCT) packages.
Transport simulation: MT3DMS reads MODFLOW outputs and solves the advection-dispersion-reaction (ADR) equation for each species across the shared spatial grid.

The governing equation solved is:

∂(θC)/∂t = ∂/∂xi[θDij ∂C/∂xj] − ∂(θviC)/∂xi + qsCs + Σ Rn

where θ is porosity, C is concentration, Dij is the hydrodynamic dispersion tensor, vi is the seepage velocity, and Rn represents reaction terms.

Advection Solvers: Choosing the Right Scheme

MT3DMS advection solver comparison: TVD, MOC, standard FD, and MMOC numerical dispersion

MT3DMS provides four advection solution methods, each with distinct accuracy-stability trade-offs:

Solver	Method	Best Use Case
MOC	Method of Characteristics	Sharp fronts, low numerical dispersion
MMOC	Modified MOC	Large grids, moderate accuracy
HMOC	Hybrid MOC	Adaptive switching between MOC/MMOC
TVD	Total Variation Diminishing (ULTIMATE)	Steep concentration gradients, mass balance

For chlorinated solvent plumes (TCE, PCE) with sharp source-zone boundaries, the TVD scheme is generally preferred. It eliminates the artificial oscillations that plague standard finite-difference upwinding while maintaining strict mass conservation—a regulatory requirement for most remediation feasibility studies. The MOC scheme excels when tracking tracer fronts in heterogeneous aquifers where numerical dispersion would otherwise mask true dispersivity values.

Reaction Package: Sorption and Biodegradation

Sequential reductive dechlorination breakthrough curves and RCT package decay chain configuration

The RCT package handles the geochemical and biological processes that govern contaminant fate:

Sorption is implemented via three isotherms:

Linear (Kd model): appropriate for low-concentration organic compounds
Freundlich: captures nonlinear partitioning at higher concentrations
Langmuir: models capacity-limited sorption on mineral surfaces

For petroleum hydrocarbons and chlorinated solvents, the retardation factor R = 1 + (ρb/θ)Kd directly controls plume velocity. A TCE plume with Kd = 0.5 mL/g in a sandy aquifer (ρb = 1.8 g/cm³, θ = 0.30) travels at roughly 40% of the groundwater velocity—a critical parameter for predicting receptor exposure timelines.

Biodegradation options include:

First-order decay (single species)
Sequential first-order decay chains (parent → daughter products)
Monod kinetics (via the RT3D extension)

The sequential decay chain capability is essential for modeling reductive dechlorination: PCE → TCE → DCE → vinyl chloride → ethene. Each daughter product is tracked as a separate species with its own retardation and decay parameters, enabling regulators to assess whether natural attenuation is producing the more toxic vinyl chloride intermediate.

Practical Workflow: DNAPL Source Zone Characterization

A common MT3DMS application involves back-calculating DNAPL (dense non-aqueous phase liquid) source zone mass from downgradient monitoring well data. The procedure:

Build the MODFLOW model: Calibrate hydraulic conductivity using head observations and slug tests. Typical K values for contaminated sand aquifers range from 10⁻⁴ to 10⁻² m/s.
Define the source term: Represent the DNAPL dissolution zone as a constant-concentration cell (SSM package, Type 1 BC) or a mass-flux boundary using field-measured source concentrations.
Calibrate dispersivity: Adjust longitudinal dispersivity (αL) and the transverse/vertical ratios (αT/αL ≈ 0.1, αV/αL ≈ 0.01) to match observed plume widths at monitoring wells.
Run sensitivity analysis: Vary Kd, first-order decay rate (λ), and source concentration to bracket plume extent under uncertainty.
Forecast remediation: Simulate pump-and-treat extraction wells (specified-flux SSM cells) to evaluate mass removal rates and cleanup timelines.

FloPy's mt3d module streamlines this workflow in Python, enabling scripted parameter sweeps and automated post-processing of concentration arrays.

Integration with SEAWAT for Density-Dependent Transport

In coastal aquifers and deep saline formations, variable-density flow invalidates the standard MODFLOW/MT3DMS decoupled approach. SEAWAT (version 4) integrates MODFLOW-2000 and MT3DMS into a fully coupled variable-density solver, where salinity-driven density differences feed back into the flow equations at each time step.

SEAWAT is the tool of choice for:

Saltwater intrusion modeling in coastal municipal wellfields
CO₂ sequestration brine displacement in deep saline aquifers
Geothermal reservoir thermal plume simulation (using temperature as the transport species)

The density coupling is controlled by the VDF (Variable-Density Flow) package, which modifies the freshwater head formulation to account for fluid density variations up to 1,025 kg/m³ (seawater).

Output Analysis and Visualization

MT3DMS writes concentration arrays to the UCN (Unformatted Concentration) binary file at user-specified time steps. Post-processing options include:

FloPy: flopy.utils.UcnFile reads UCN arrays directly into NumPy for Python-based visualization
ModelMuse: USGS GUI providing plan-view and cross-section concentration maps
Groundwater Vistas: Commercial interface with built-in particle tracking and mass budget reporting
PEST/PEST++: Automated parameter estimation using MT3DMS as the forward model for concentration-based calibration

Mass balance reporting via the MAS file is essential for verifying simulation integrity. A global mass balance error exceeding 1% typically indicates time-step control issues (reduce the Courant number below 0.75 for TVD) or grid resolution problems near high-flux boundaries.

Key Resources

MT3DMS continues to underpin regulatory groundwater investigations worldwide, from Superfund site remediation feasibility studies to drinking water protection zone delineation. Its transparent, modular design and deep integration with the MODFLOW ecosystem make it an indispensable tool for any hydrogeologist working at the intersection of contaminant fate, transport, and remediation engineering.