The ECOSIM groundwater quality models: MT3D

MT3D model, is a model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater flow systems in either two or three dimensions. The model uses a mixed Eulerian-Lagrangian approach to the solution of the advective-dispersive-reactive equation, based on combination of the method of characteristics and the modified method of characteristics. This approach combines the strength of the method of characteristics for eliminating numerical dispersion and the computational efficiency of the modified method of characteristics. The model program uses a modular structure similar to that implemented in the U.S. Geologic Survey modular three-dimensional finite-difference groundwater flow model, referred to as MODFLOW, (McDonald and Harbaugh, 1988). The modular structure of the transport mode makes it possible to simulate advection, dispersion, source/sink mixing, or chemical reactions independently without reserving computer memory space for unused options; new packages involving other transport processes can be added to the model readily without having to the existing code.

The MT3D transport model was developed for use with any block-centered finite- difference flow model such as MODFLOW and is based on the assumption that changes in concentration field will not affect the flow field significantly. After a flow model is developed and calibrated, the information needed by the transport model can be saved in disk files which are then retrieved by the transport model.

The MT3D transport model can be used to simulate changes in concentration of single-species miscible contaminants in groundwater considering advection, dispersion and some simple chemical reactions, with various types of boundary conditions and external sources or sinks. The chemical reactions included in the model are equilibrium-controlled linear or non-linear sorption and first-order irreversible decay or biodegradation. Currently, MT3D accommodates the following spatial discretization capabilities and transport boundary conditions:

1. confined, unconfined or variably confined/unconfined aquifer layers;
2. inclined model layers and variable cell thickness within the same layer;
3. specified concentration or mass flux boundaries; and
4. the solute transport effects of external sources and sinks such as wells, drains, rivers, areal recharge and evapotranspiration.

Solution techniques

The advective-dispersive-reactive equation describes the transport of miscible contaminants in groundwater flow systems. Most numerical methods for solving the advective-dispersive-reactive equation can be classified as Eulerian, Lagrangian or mixed Eulerian-Lagrangian (Neuman 1984). In the Eulerian approach, the transport equation is solved with a fixed grid method such as the finite-difference or finite- element method. The Eulerian approach offers the advantage and convenience of a fixed grid, and handles dispersion/reaction dominated problems effectively. For advection-dominated problems which exist in many field conditions, however, an Eulerian method is susceptible to excessive numerical dispersion or oscillation, and limited by small grid spacing and time steps. In the Lagrangian approach, the transport equation is solved in either a deforming grid or deforming coordinate in a fixed grid. The Lagrangian approach provides an accurate and efficient solution to advection dominated problems with sharp concentration fronts. However, without a fixed grid or coordinate, a Lagrangian method can lead to numerical instability and computational difficulties in nonuniform media with multiple sinks/sources and complex boundary conditions (Yeh, 1990). The mixed Eulerian-Lagrangian approach attempts to combine the advantages of both the Eulerian and the Lagrangian approaches by solving the advection term with a Lagrangian method and the dispersion and reaction terms with an Eulerian method.

The numerical solution implemented in MT3D is a mixed Eulerian-Lagrangian method. The Lagrangian part of the method, used for solving the advection term, employs the forward tracking method of characteristics (MOC), the backward- tracking modified method of characteristics (MMOC), or a hybrid of these two methods. The Eulerian part of the method, used for solving the dispersion and chemical reaction terms, utilizes a conventional block-centered finite-difference method.

The MT3D transport model uses an explicit version of the block-centered finite- difference method to solve the dispersion and chemical reaction terms. The limitation of an explicit scheme is that there is a certain stability criterion associated with it, so that the size of time steps cannot exceed a certain value. However, the use of an explicit scheme is justified by the fact that it saves a large amount of computer memory which would be required by a matrix solver used in an implicit scheme. In addition, for many advection-dominated problems, the size of transport steps is dictated by the advection process, so that the stability criterion associated with the scheme for the dispersion and reaction processes is not a factor. It should be noted that a solution package based on implicit schemes for solving dispersion and reactions could easily be developed and added to the model as an alternative solver for mainframes, more powerful personal computers, or workstations with less restrictive memory constraints.

The implementation of the model consists of the following actions: definition of the simulation problem, specification of the initial and boundary conditions, determination of appropriate transport stepsize, preparation of global mass balance information and output of simulation results. The Flow Model Interface Package interfaces with a flow model to obtain the flow solution from the flow model.

The primary modules of MT3D are:

Basic BTN handles basic tasks that are required by the entire transport model. Among these tasks are definition of the problem, specification of the boundary and initial conditions, determination of the stepsize, preparation of mass balance information, and printout of the simulation results.

Flow Model FMI Interfaces with a flow model. Currently, the interfacing is done through an unformatted disk file containing heads and flow terms. The FMI Package reads the contents of this file and prepares heads and flow terms in the form needed by the transport model.

Advection ADV Solves the concentration change due to advection with one of the three mixed Eulerian-Lagrangian schemes included in the package: MOC, MMOC, or HMOC. Dispersion DSP Solves the concentration change due to dispersion with the explicit finite difference method.

Sink & Source SSM Solves the concentration change due to fluid Mixing sink/source mixing with the explicit finite difference method. Sink/source terms may include wells, drains, rivers, recharge and evapotranspiration. The constant-head boundary and general-head-dependent boundary are also handled as sink/source terms in the transport model.

Chemical RCT Solves the concentration change due to chemical reactions. Currently, the chemical reactions include linear or nonlinear sorption isotherms and first-order irreversible rate reactions (radioactive decay or biodegradation).

Utility UTL Contains a number of utility modules that are called upon by primary modules to perform such general-purposed tasks as input/output of data arrays.