Geology 727

MT3DMS Solution Methods and Parameter Options

Introduction

Please refer to Table 3.1 on p. 3-20 of the manual (also attached) for a summary of the solution methods available in MT3DMS. Note that the there are two main columns in the table. One column refers to methods used to solve the troublesome advection component of the advection-dispersion equation and the other column refers to solution methods for the rest of the equation (dispersion, sink/source, and chemical reaction terms). It is useful to think in terms of solving the advection component and thus to think about four basic solution approaches:

  1. MOC/MMOC/HMOC (with either explicit or implicit finite differences for the rest of the equation);
  2. Explicit finite difference method;
  3. Implicit finite difference method;
  4. TVD method (with either explicit or implicit finite differences for the rest of the equation).

Time step. The time step is automatically selected by the code to meet various stability constraints that are solution dependent. By setting the transport time step size (DT0) equal to zero, the code will automatically select a time step. You may select a different value for the time step by entering a non-zero value. You must, however, select a value less than the value selected by the code or else it will default to the code-selected value.

Courant Number. The Courant number (Cr) is used in all solution methods as a way to ensure that the time step selected by MT3DMS is appropriate for the problem. If the Courant number is set equal to one, the time step will equal the grid spacing divided by the velocity. This means that the solution will be calculated for the time it takes for the solute to move across one cell. MT3DMS selects the time step as follows:

 t = R Cr min (x/vx, y/vy, z/vz)

The Courant number is typically set to a value between 0.5 and 1.0. Sometimes a value up to two is used. But for the explicit finite difference option and the TVD method, the Courant number is also a stability constraint and must not exceed one.

Peclet Number. The Peclet number (Pe) measures the ratio of advection to dispersion. Thus, an advection dominant problem has a high Peclet number.

Pe = v x / D

or

Pe = x / 

Lower order finite difference schemes in MT3DMS (i.e., all finite difference schemes except for the TVD method) show numerical errors (known as numerical dispersion) for advection dominant problems. Hence, these methods should be used only when the Peclet number is small, usually less than four. This can be severely limiting since it limits the magnitude of the permissible grid spacing.

MOC and HMOC solution methods and the TVD method are not subject to numerical dispersion and can be used for any value of Peclet number.

MOC/MMOC/HMOC Methods

These methods can be used to solve the advection component of the advection-dispersion equation. MOC starts the simulation with a fixed number of particles per cell. MMOC uses one particle per cell. HMOC uses a combination of MOC and MMOC; MOC is used near sharp concentration fronts and MMOC is used away from such fronts. The MMOC technique should be used only in situations where sharp concentrations fronts are not present.

In the MOC/MMOC/HMOC methods either the explicit or the implicit finite difference method is used to solve the dispersion term, sink/source term, and the reaction term. If the explicit finite difference method is used, various stability constraints arise which set limits on the time step. See p. 4-20 of the manual for details.

Parameters in the MOC/MMOC/HMOC methods

WD is the concentration weighting factor, which ranges between 0 and 1. The default value is 0.5. This can be adjusted to get better mass balance. Increase toward 1.0 for advection dominant problems.

NPH is the number of particles per cell. Use a maximum of 32 for 3D problems. Can reduce to 16 for 2D problems.

Number of particles refers to the maximum particles allowed. The default value is 5000. If you keep introducing particles (solute) into the model domain without removing particles, the solution may run out of particles, i.e., all 5000 particles will be in use in the model domain and there won’t be any particles left to continue the introduction of solute. This can create numerical problems.

Particle Planes refer to the particle placement pattern at the beginning of the simulation. If set to 0, a random pattern is selected. If a value greater than zero is selected, a fixed pattern is used. See p. 6-28 of the manual for details. The parameter name is NPLANE.

ITRACK selects the particle tracking scheme. Options are a simple 1st order Euler scheme, the more accurate 4th order Runge Kutta scheme, or a combination such that Runge Kutta is used near sink and source cells and Euler is used everywhere else.

The small concentration gradient (DCEPS) is the gradient below which advective transport is considered negligible, usually set around 1E-5.

DCHMOC is the critical relative concentration gradient for controlling the selective use of either MOC or MMOC in the HMOC solution. The default value is 0.005.

Finite Difference Methods

Summary

The explicit method is influenced by numerical dispersion and has many restrictions on the time step. The implicit method is subject to numerical dispersion when the advection term is solved with upstream weighting; alternatively the method is subject to oscillations when the advection term is solved with central weighting. It has no stability limit on the time step. Numerical errors in both methods can be minimized by keeping the Peclet number small (less than four). The TVD method uses an explicit method to solve the advection term and thus has a stability constraint on the time step but it is relatively free from numerical errors and can be used for all values of the Peclet number.

Methods

The explicit method uses upstream weighting for the advection term since the central weighting scheme is unstable for this method. The explicit method has stability contraints on the time step: 1) the Courant number must be one or less; 2) see p 4-20 of the manual for additional stability constraints.

The implicit method uses either upstream weighting or central weighting for the advection term. The time step is controlled by a time step multiplier TSMULT. To use a constant value for the time step, set the multiplier equal to one. The generalized conjugate gradient (GCG) solver is used in this method. There is a head change closure criterion on the solver, for which the default value is 1E-6. There is a choice of three iterative algorithms: Jacobi, SSOR, and Modified Incomplete Cholesky (MIC). MIC usually takes fewer iterations to converge but requires more memory.

The explicit 3rd order TVD method uses an explicit finite difference approximation for the advection term. Therefore, the Courant number must be less than or equal to one. The method works well for all values of Peclet number. Either standard explicit finite differences or implicit finite differences are selected to solve the dispersion, sink/source, and chemical reaction terms. If explicit finite differences are selected, there are additional constraints on the time step. See p. 4-20 of the manual for details.

Table 3.1 from MT3DMS manual

1