A Three-Dimensional 

Method-of-Characteristics 

Solute-Transport Model (MOC3D)

By L.F. Konikow, D.J. Goode, and G.Z. Hornberger

U. S. GEOLOGICAL SURVEY

Water-Resources Investigations Report 96-4267

Reston, Virginia 

1996
U.S. DEPARTMENT OF THE INTERIOR 

BRUCE BABBITT, Secretary

U.S. GEOLOGICAL SURVEY 

Gordon P. Eaton, Director

The use of firm, trade, and brand names in this report is for identification purposes 

only and does not constitute endorsement by the U.S. Geological Survey.

For additional information write to:

Leonard F. Konikow 

U.S. Geological Survey 

Water Resources Division 

431 National Center 

Reston, VA20192

Copies of this report can be purchased 

from:

U.S. Geological Survey 

Information Services 

Box 25286 

Federal Center 

Denver, CO 80225
PREFACE

The computer model described in this report is designed to simulate the transport and 

dispersion of a single solute in ground water flowing through porous media. The model is 

developed as a module for the U.S. Geological Survey's (USGS) MODFLOW ground-water 

model, and it is intended to be the first of a family of alternative solution methods for the solutetransport equation that will be compatible with MODFLOW.

This model, named MOC3D, was developed through modifications of an existing twodimensional code (MOC), which was documented originally by Konikow and Bredehoeft (1978). 

Although extensive testing ofMOCSD indicates that this model will yield reliable calculations for a 

wide variety of field problems, the user is cautioned that the accuracy and efficiency of the model 

can be affected significantly for certain combinations of values for parameters and boundary 

conditions. Development of alternative codes that will optimize the accuracy and efficiency of 

solving the solute-transport equation for a broader range of conditions is planned.

The code for this model will be available for downloading over the Internet from a USGS 

software repository. The repository is accessible on the World Wide Web (WWW) from the 

USGS Water Resources Information web page at URL http://h2o.usgs.gov/. The URL for the 

public repository is: http://h2o.usgs.gov/software/. The public anonymous FTP site is on 

the Water Resources Information server (h2o.usgs.gov or 130.11.50.175) in the /pub/software 

directory. When this code is revised or updated in the future, new versions or releases will be 

made available for downloading from these same sites.

Acknowledgments. The authors appreciate the helpful model evaluation and review, 

comments provided by USGS colleagues H. I. Essaid, W. B. Fleck, and S. P. Garabedian.

in
CONTENTS

Page

ABSTRACT ................................................................................................................................ 1

INTRODUCTION.....................^ 1

THEORETICAL BACKGROUND AND GOVERNING EQUATIONS ........................................ 3

Ground-Water Flow Equation ......................................................................................... 3

Average Interstitial Velocity ............................................................................................ 4

Governing Equation for Solute Transport ....................................................................... 4

Dispersion Coefficient .................................................................................................... 6

Review of Assumptions ................................................................................................... 8

NUMERICAL METHODS .......................................................................................................... 8

Ground-Water Flow Equation ......................................................................................... 9

Average Interstitial Velocity ............................................................................................ 10

Solute-Transport Equation .............................................................................................. 11

Method of Characteristics ................................................................................... 11

Particle Tracking ................................................................................................ 14

Linear interpolation and semi-analytic integration .................................. 15

Bilinear interpolation and explicit integration ......................................... 18

Discussion Choosing an appropriate interpolation scheme .................... 20

Decay ................................................................................................................. 20

Node Concentrations .......................................................................................... 21

Finite-Difference Approximations ...................................................................... 21

Stability and Accuracy Criteria ........................................................................... 23

Mass Balance ...................................................................................................... 25

Special Problems ................................................................................................ 26

Review of MOC3D Assumptions and Integration with MODFLOW ...................... 30

COMPUTER PROGRAM............................................................................................................ 33

General Program Features ............................................................................................... 34

Program Segments .......................................................................................................... 35

MODEL TESTING AND EVALUATION .................................................................................. 41

One-Dimensional Steady Flow ........................................................................................ 41

Three-Dimensional Steady Flow ..................................................................................... 45

Two-Dimensional Radial Flow and Dispersion ................................................................ 49

Point Initial Condition in Uniform Flow .......................................................................... 50

Constant Source in Nonuniform Flow ............................................................................. 55

Relative Computational Efficiency .................................................................................. 60

CONCLUSIONS ......................................................................................................................... 62

REFERENCES ............................................................................................................................ 62

APPENDIX A

FINITE-DIFFERENCE APPROXIMATIONS ................................................................. 64

APPENDIX B

DATA INPUT INSTRUCTIONS FOR MOC3D ............................................................... 70
Page 

MODFLOW Name File .................................................................................................... 70

MODFLOW Source and Sink Packages ........................................................................... 70

MOC3D Input Data Files ................................................................................................. 70

MOC3D Transport Name File (CONC)................................................................ 71

Main MOC3D Package Input (MOC) .................................................................. 72

Source Concentration in Recharge File (CRCH) .................................................. 77

Observation Well File (OBS) ............................................................................... 77

APPENDIX C

ANNOTATED EXAMPLE INPUT DATA SET FOR SAMPLE PROBLEM .................... 78

APPENDIX D

SELECTED OUTPUT FOR SAMPLE PROBLEM .......................................................... 82

FIGURES

1-3. Diagrams showing:

1. Notation used to label rows, columns, and nodes within one layer (k) of a threedimensional, block-centered, finite-difference grid for MOC 3D ............................. 9

2. Representative three-dimensional grid for MOC3D illustrating notation for layers ... 9

3. Part of a hypothetical finite-difference grid showing relation of flow field to 

movement of points (or particles) in method-of-characteristics model for 

simulating solute transport (modified from Konikow and Bredehoeft, 1978) ......... 12

4. Graph showing representative change in breakthrough curve from time level / to M-l. 

Note that concentration changes are exaggerated to help illustrate the sequence of 

calculations. Curve for Ct+adv represents the concentration distribution at time M-l due 

to advection only. (Modified from Konikow and Bredehoeft, 1978.)........................... 13

5-11. Diagrams showing:

5. Spatial weights used in linear interpolation method to estimate Vx at location of a

particle in ce\\j,i,k ................................................................................................. 16

6. Interpolation factors used in bilinear interpolation method to estimate horizontal 

components of velocity, Vx and Vy, at the position of a particle located in the 

southeast quadrant of cell j,i,k ................................................................................ 19

7. Possible movement of a particle near an impermeable (no-flow) boundary

(modified from Konikow and Bredehoeft, 1978) ................................................... 27

8. Replacement of particles in fluid-source cells (a) for case of negligible regional

flow and (b) for case of relatively strong regional flow .......................................... 28

9. One layer of finite-difference grid illustrating use of uniformly-spaced transport 

subgrid for MOC3D within variably-spaced primary grid for MODFLOW .............. 31

VI
Page

10. Alternative MODFLOW approaches to vertical discretization of an aquifer system 

(a) consisting of two high-permeability units separated by a confining layer 

consisting of a low-permeability clay. In the fully three-dimensional 

representation (b), the clay unit would be represented by one or more model 

layers. In the quasi-three-dimensional approach (c), heads are not calculated in 

the clay unit, which is represented more simply by vertical conductance terms 

between the layers above and below it; in this case, the bottom of layer 1 coincides 

with the top of layer 2. (Modified from figs. 10-12 of McDonald and Harbaugh, 

1988.) ................................................................................................................... 32

11. Double time-line illustrating the sequence of progression in the MOC3D model for 

solving the flow and transport equations. This example is for transient flow and 

three stress periods (NPER = 3) of durations PERLENj, PERLEN2, and PERLEN3 . 

Each time step for solving the flow equation (of duration DELT) is divided into 

one or more time increments (of duration TIMV) for solving the transport equation; all particles are moved once during each transport time increment. For 

illustration purposes, the sequence of solving the two equations is labeled for the 

first five time steps of the first stress period, and the indices for counting time steps 

for flow and time increments for transport are labeled for the fourth time step ....... 34

12. Generalized flow chart for MOC3D ................................................................................ 36

13. Simplified flow chart for the transport loop, which is shown as a single element in

fig. 12 .......................................................................................................................... 37

14-16. Diagrams showing:

14. Default initial particle configurations for a one-dimensional simulation using

(a) one, (b) two, (c) three, and (d) four particles per cell......................................... 38

15. Default initial particle configurations for a two-dimensional simulation using

(a) one, (b) four, (c) nine, and (d) sixteen particles per cell .................................... 38

16. Default initial particle configurations for a three-dimensional simulation using 

(a) one, (b) eight, and (c) 27 particles per cell. For clarity, only the volume of a 

single cell is illustrated ........................................................................................... 38

17. Flow chart for MOVE and MOVEBI subroutines of MOC3D, which is shown as a single

element in fig. 13 ......................................................................................................... 39

18-24. Graphs showing:

18. Numerical (MOC3D) and analytical solutions at three different locations for solute 

transport in a one-dimensional, steady flow field. Parameter values for this base 

case are listed in Table 11 ...................................................................................... 43

19. Numerical (MOC3D) and analytical solutions at three different locations for solute 

transport in a one-dimensional, steady flow field for case of increased dispersivity 

(aL = 1.0 cm, Dxx = 0.1 cm2/s, and other parameters as defined in Table 11) ......... 44

20. Detailed view of numerical and analytical solutions for early times (t < 10 s) at the 

first node downgradient from the inflow source boundary for same problem as 

shown in fig. 19 ..................................................................................................... 44

21. Detailed view of numerical and analytical solutions for early times (t < 10 s) at the 

first node downgradient from the inflow source boundary for same conditions as 

shown in fig. 20, except that the initial number of particles per node, NPTPND, 

equals 50 ............................................................................................................... 44

Vll
Page

22. Numerical (MOC3D) and analytical solutions at three different times for same onedimensional, steady flow, solute-transport problem shown in fig. 18 ...................... 45

23. Numerical (MOC3D) and analytical solutions after 240 seconds for three different

retardation factors for same problem represented in fig. 22 ................................... 45

24. Numerical (MOC3D) and analytical solutions for four different times for solute 

transport in a one-dimensional, steady flow field for case with decay at rate of 

A = 0.01 s' 1 . All other parameters as defined in Table 11 ...................................... 45

25-32. Plots showing:

25. Concentration contours for (a) analytical and (b) numerical solutions in the horizontal plane containing the solute source (layer 1) for three-dimensional solute 

transport in a uniform steady flow field. Parameters are defined in Table 12......... 47

26. Concentration contours for (a) analytical and (b) numerical solutions in the vertical 

plane parallel to the flow direction and aligned with column 2 for threedimensional solute transport in a uniform steady flow field. Parameters are 

defined in Table 12 ............................................................................................... 48

27. Concentration contours for (a) analytical and (b) numerical solutions in the vertical 

plane transverse to the flow direction and aligned with row 12 for threedimensional solute transport in a uniform steady flow field. Parameters are 

defined in Table 12. Flow is towards the reader ..................................................... 48

28. Initial particle positions within the source cell for radial flow case (based on custom 

particle placement and NPTPND = 46). The relative coordinates on the x- and yaxes shown for the cell (1,1) are the same for any cell of the grid; this relative 

coordinate system is used for the custom definition of particle locations in the 

input file ............................................................................................................... 50

29. Contours of relative concentrations calculated using (a) analytical and (b) 

numerical models for solute transport in a steady radial flow field. Source 

concentration is 1.0 and source is located in cell (1,1). Grid spacing is 10.0 m ...... 50

30. Concentration contours for (a, c) analytical and (b, d) numerical solutions for 

transport of a point initial condition in uniform flow in the x-direction. The zcomponent of flow is zero, but there is dispersion in all three directions. Contour 

values are the log of the concentrations .................................................................. 53

31. Concentration contours for (a) analytical and (b, c, d) numerical solutions for 

transport of a point initial condition in uniform flow at 45 degrees to x and y. 

Contour values are the log of the concentrations .................................................... 54

32. Concentration contours showing effects on areas of negative concentrations related 

to decreasing CELDIS factor in MOC3D in simulation of flow at 45 degrees to 

grid having 24 rows and 24 columns of nodes: (a) CELDIS = 0.50; 

(b) CELDIS = 0.25; and (c) CELDIS = 0.10 ......................................................... 55

33. Schematic representation of transport domain and boundary conditions for

nonuniform-flow test problem of Burnett and Frind (1987); front surface represents a 

plane of symmetry (modified from Burnett and Frind, 1987, fig. 4) ............................. 56

Vlll
Page 

34-38. Plots showing:

34. Initial particle positions within a cell for the Burnett-Frind (1987) test case (based 

on custom particle placement and NPTPND = 3). The relative coordinates on the 

x- and y-axes shown for the cell (1,1) are the same for any cell of the grid; this 

relative coordinate system is used for the custom definition of particle locations in 

the input file .......................................................................................................... 57

35. Two-dimensional simulation results for nonuniform-flow test case showing plume 

positions as contours of relative concentration: (a) finite-element model (modified 

from Burnett and Frind, 1987, fig. 8a), and (b) MOC3D. Contour interval is 0.2 

relative concentration ............................................................................................. 58

36. Three-dimensional simulation results for nonuniform-flow test case in which

OCTH = 0.1 m and CX,jv = 0.01 m: (a) finite-element model (modified from Burnett 

and Frind, 1987, fig. 8c), and (b) MOC3D. Plume positions are represented by 

contours of relative concentration; contour interval is 0.2 relative concentration .... 58

37. Perspective view of MOC3D results for three-dimensional problem of constant 

source in nonuniform flow and unequal transverse dispersivity coefficients. This 

visualization of the plume was generated from particle concentrations using a 

three-dimensional visualization software package and is derived from the same 

simulation that is the basis of fig. 36b. Note that a piece of the plume near the 

source is cut away (a "chair" cut) to expose a clearer view of the degree of 

transverse spreading in the selected vertical and horizontal planes. Shading 

increments are in relative-concentration intervals of 0.20, and interval bounds 

range between 0.9 and 0.1. Concentrations less than 0.10 are transparent ............. 59

38. Three-dimensional simulation results for nonuniform-flow test case in which

(XjH = (Xjv = 0.01 m: (a) finite-element model (modified from Burnett and Frind, 

1987, fig. 9b), and (b) MOC3D, Plume positions are represented by contours of 

relative concentration; contour interval is 0.2 relative concentration ....................... 60

TABLES

\.MOC3D subroutines controlling simulation preparation ................................................. 37

2. MOC3D subroutines controlling transport time factors ................................................... 37

3.MOC3D subroutines controlling velocity calculations and output................................... 38

4. MOC3D subroutines controlling particle tracking, concentration calculations, and

output for particle and concentration data .................................................................... 38

5.MOC3D subroutines controlling dispersion calculations and output ............................... 40

6. MOC3D subroutines controlling MODFLOW source/sink package calculations .............. 40

7. MOC3D subroutines controlling cumulative calculations relating to sources and sinks ... 40

8.MOC3D subroutines controlling observation wells ......................................................... 40

9. MOC3D subroutines controlling mass balance calculations and output........................... 41

10. Miscellaneous MOC3D subroutines ................................................................................ 41

IX
Page

11. Parameters used in MOC3D simulation of transport in a one-dimensional, steady-state

flow system .................................................................................................................. 42

12. Parameters used in MOC3D simulation of transport from a continuous point source in a

three-dimensional steady-state flow system ................................................................... 47

13. Parameters used in MOC3D simulation of two-dimensional, steady-state, radial flow

case .............................................................................................................................. 49

14. Parameters used in MOC3D simulation of three-dimensional transport from a point

source with flow in the x-direction and flow at 45 degrees to x and y ............................ 52

15. Parameters used in MOC3D simulation of transport in a vertical plane from a

continuous point source in a nonuniform, steady-state, two-dimensional flow system 

described by Burnett and Frind (1987) ......................................................................... 57

16. Comparison of MOC3D simulation times for selected test cases on various computer

platforms ...................................................................................................................... 61

17. Formats associated with MOC3D print flags ................................................................... 74
ABSTRACT

This report presents a model, MOC3D, that simulates three-dimensional solute transport 

in flowing ground water. The model computes changes in concentration of a single dissolved 

chemical constituent over time that are caused by advective transport, hydrodynamic dispersion 

(including both mechanical dispersion and diffusion), mixing (or dilution) from fluid sources, 

and mathematically simple chemical reactions (including linear sorption, which is represented 

by a retardation factor, and decay). The transport model is integrated with MODFLOW, a 

three-dimensional ground-water flow model that uses implicit finite-difference methods to 

solve the transient flow equation. MOC3D uses the method of characteristics to solve the 

transport equation on the basis of the hydraulic gradients computed with MODFLOW for a 

given time step. This implementation of the method of characteristics uses particle tracking to 

represent advective transport and explicit finite-difference methods to calculate the effects of 

other processes. However, the explicit procedure has several stability criteria that may limit the 

size of time increments for solving the transport equation; these are automatically determined by 

the program. For improved efficiency, the user can apply MOC3D to a subgrid of the primary 

MODFLOW grid that is used to solve the flow equation. However, the transport subgrid must 

have uniform grid spacing along rows and columns. The report includes a description of the 

theoretical basis of the model, a detailed description of input requirements and output options, 

and the results of model testing and evaluation. The model was evaluated for several problems 

for which exact analytical solutions are available and by benchmarking against other numerical 

codes for selected complex problems for which no exact solutions are available. These test 

results indicate that the model is very accurate for a wide range of conditions and yields 

minimal numerical dispersion for advection-dominated problems. Mass-balance errors are 

generally less than 10 percent, and tend to decrease and stabilize with time.

INTRODUCTION

This report describes and documents a 

computer model (MOC3D) for calculating 

transient changes in the concentration of a 

single solute in a three-dimensional groundwater flow field. The calculations require the 

numerical solution of two simultaneous partial 

differential equations. One equation is the 

ground-water flow equation, which describes 

the head distribution in the aquifer. The second 

is the solute-transport equation, which 

describes the solute concentration within the 

flow system. By coupling the flow equation 

with the solute-transport equation, the model 

can be applied to both steady-state and transient 

ground-water flow problems.

The purpose of the simulation model is 

to compute the concentration of a dissolved 

chemical species in an aquifer at any specified 

place and time. Changes in chemical 

concentration occur within a dynamic groundwater system primarily because of four distinct 

processes: (1) advective transport, in which 

dissolved chemicals are moving with (are being 

carried by) the flowing ground water; (2) 

hydrodynamic dispersion, in which molecular 

and ionic diffusion and mechanical dispersion 

(related mostly to variations in fluid velocity 

through the porous media) cause the paths of 

dissolved molecules and ions to diverge and 

spread from the average direction of ground-
water flow; (3) fluid sources, where water of 

one composition is introduced into and mixes 

with water of a different composition; and (4) 

reactions, in which some amount of the solute 

is added to or removed from the ground water 

because of chemical, biological, and (or) 

physical reactions in the water or between the 

water and the solid aquifer materials.

The MOC3D model is integrated with 

MODFLOW, the U.S. Geological Survey's 

(USGS) modular, three-dimensional, finitedifference, ground-water flow model 

(McDonald and Harbaugh, 1988; Harbaugh 

and McDonald, 1996a and 1996b). 

MODFLOW solves the ground-water flow 

equation, and the reader is referred to the 

documentation for that model and its 

subsequent modules for complete details. In 

this report it is assumed that the reader is 

familiar with the MODFLOW family of codes. 

The numerical solution to the solute-transport 

equation is based on the method-of-characteristics, which is advantageous (relative to other 

standard numerical schemes) for transport 

problems in which advection is a dominant 

process. The MOC3D code is based largely on 

MOC, the USGS's two-dimensional, methodof-characteristics, solute-transport model 

(Konikow and Bredehoeft, 1978; Goode and 

Konikow, 1989).

This model can be applied to a wide 

variety of field problems. However, the user 

should first become aware of the assumptions 

and limitations inherent in the model, as 

described in this report. MOC3D is offered as 

a basic tool that is applicable to a wide range of 

field problems involving solute transport. 

However, there are some situations in which 

the model results could be inaccurate or model 

operation inefficient. The report includes 

guidelines for recognizing and avoiding these 

types of problems.

The computer program is written in 

FORTRAN, and has been developed in a 

modular style, similar to the MODFLOW

model. At present, the model is not compatible 

with all the modules for MODFLOW that 

describe secondary flow processes or features, 

such as streamflow routing, subsidence, and 

rewetting of dry cells. As assumed by 

MODFLOW, it is also assumed by MOC3D 

that fluid properties are homogeneous and that 

concentration changes do not significantly 

affect the fluid density or viscosity, and hence 

the fluid velocity. Within the finite-difference 

grid used to solve the flow equation in 

MODFLOW, the user is able to specify a 

window or subgrid over which MOC3D will 

solve the solute-transport equation. This 

feature can significantly enhance the overall 

efficiency of the model by avoiding calculation 

effort where it is not needed. However, 

MOC3D also requires that the horizontal (row 

and column) grid spacing be uniform within the 

subgrid, although an expanding or nonuniform 

spacing may be applied outside of the subgrid 

boundaries.

The types of reactions incorporated into 

MOC3D are restricted to those that can be 

represented by a first-order rate reaction, such 

as radioactive decay, or by a retardation factor, 

such as instantaneous, reversible, sorptiondesorption reactions governed by a linear 

isotherm and constant distribution coefficient 

(Kd). This is somewhat more restrictive than 

the two-dimensional model, MOC, which 

allowed the representation of nonlinear 

isotherms.

The report includes a detailed 

description of the numerical methods used to 

solve the solute-transport equation. The data 

requirements, input format specifications, 

program options, and output formats are all 

structured in a general manner that should be 

compatible with the types of data available for 

many field problems. We have attempted to 

maximize the use of existing MODFLOW 

output modules and styles in developing the 

MOC3D output options and features.
THEORETICAL BACKGROUND 

AND GOVERNING EQUATIONS

Mathematical equations that describe 

ground-water flow and transport processes 

may be developed from the fundamental 

principle of conservation of mass of fluid or of 

solute. A statement of conservation of mass 

(or continuity equation) may be combined with 

a mathematical description of the relevant 

process to obtain a differential equation 

describing flow or transport (see, for example, 

Bear, 1979; Freeze and Cherry, 1979; or 

Domenico and Schwartz, 1990).

Ground-Water Flow Equation

A quantitative description of groundwater flow is a prerequisite to accurately 

simulating solute transport in aquifers. A 

general form of the equation describing the 

transient flow of a compressible fluid in a 

heterogeneous anisotropic aquifer may be 

derived by combining Darcy's Law with the 

continuity equation. A general ground-water 

flow equation may be written in Cartesian 

tensor notation as:

where KIJ is the hydraulic conductivity of the 

porous media (a second-order tensor), LT' 1 ; h is 

the hydraulic (or potentiometric) head, L; Ss is 

the specific storage, L' 1 ; t is time, T; W is the 

volumetric flux per unit volume (positive for 

inflow and negative for outflow), T' 1 ; and xi are 

the Cartesian coordinates, L. The summation 

convention of Cartesian tensor analysis is 

implied in eq. 1. Equation 1 can generally be 

applied if isothermal conditions prevail, the 

porous medium deforms only vertically, the 

volume of individual grains remains constant 

during deformation, Darcy's Law applies (and 

gradients of hydraulic head are the only driving

force), and fluid properties (density and 

viscosity) are homogeneous and constant.

In general, the properties of porous 

media vary in space. Although fluid sources 

and sinks may vary in space and time, they 

have been lumped into one term (W) in the 

previous development for convenience in 

notation. Also, at any given location, more 

than one process may be adding or removing 

fluid simultaneously from the system, such as 

well withdrawals, diffuse recharge from 

precipitation, and evapotranspiration. 

However, the solution to the governing 

equation depends only on the net flux from 

sources and sinks as a function of time at each 

location.

If the principal axes of the hydraulic 

conductivity tensor are aligned with the x-y-z 

coordinate axes, the cross-product terms of the 

hydraulic conductivity tensor are eliminated; 

that is, KIJ = 0 when i *j. The ground-water 

flow equation may then be written to include 

explicitly all hydraulic conductivity terms as:

2-\K ^

dy( "dy.

\ f I TT 7 C* / O \ 

-v -**-T7 -v '' "*^n-v V '

dz dt

(1) where Kxx, Kyy, and Kzz are values of 

hydraulic conductivity along the x, y, and z 

coordinate axes, LT' 1 . Equation 2 is identical 

to eq. 1 of McDonald and Harbaugh (1988, p. 

2-1).

Under these same assumptions, 

Darcy's law may be written as:

dx (3a)

(3b)

where q is the specific discharge, LT' 1 .
Average Interstitial Velocity

The migration and mixing of chemicals 

dissolved in ground water will obviously be 

affected by the velocity of the flowing ground 

water. The specific discharge, q^ calculated 

from eq. 3 represents the volumetric flux per 

unit cross-sectional area. Thus, to calculate the 

actual average interstitial velocity of ground 

water, one must account for the actual crosssectional area through which flow is occurring. 

This is usually calculated as follows:

(4)

where V, is the average interstitial velocity (or 

seepage velocity), LT" 1 ; and e is the effective 

porosity (dimensionless) of the porous 

medium. Assuming the same grid alignment as 

stated for eq. 2, it can also be written in terms 

of Darcy's law as:

v, = - (5)

Governing Equation for Solute 

Transport

The principle of conservation of mass 

requires that the net mass of solute entering and 

leaving a specified volume of aquifer during a 

given time interval must equal the accumulation 

or loss of mass stored in that volume during the 

interval. This relation may be expressed 

mathematically in a general governing equation 

for solute transport in three dimensions in an 

incompressible fluid flowing through a porous 

medium as (see Bear, 1979, p. 239-243; and 

Goode and Konikow, 1989):

where e is porosity, C is volumetric concentration (mass of solute per unit volume of fluid, 

ML"3 ), pt> is the bulk density of the aquifer 

material (mass c>f solids per unit volume of 

aquifer, ML"3), C is the mass concentration of 

solute sorbed on or contained within the solid 

aquifer material (mass of solute per unit mass 

of aquifer material, MM" 1 ), V is a vector of 

interstitial fluid velocity components (LT" 1 ), D 

is a second-rank tensor of dispersion 

coefficients (L2^ 1 ), W is a volumetric fluid 

sink (W<0) or fluid source (Wi>0) rate per unit 

volume of aquifer (T" 1 ), C' is the volumetric 

concentration in the sink/source fluid (ML"3 ), 

and A is the decay rate (T" 1 ).

The decay term in eq. 6 typically 

represents radioactive decay of both the free 

and sorbed solute. A radioactive decay rate is 

usually expressed as a half-life (^2)- The 

half-life is the time required for the 

concentration to decrease to one-half of the 

original value, and is related to the decay rate 

as:

_ (In 2) (7)

(6)

In limited cases, the decay term can also 

adequately represent chemical decomposition or 

biodegradation. However, if in these latter 

cases there is also a sorbed phase present, it 

must be assured that the decay process occurs 

at the same rate for both the dissolved and 

sorbed phases, an assumption that is true for 

radioactive decay.

The concentration in the fluid leaving 

the aquifer at fluid sinks is commonly assumed 

to have the same concentration as the fluid in 

the aquifer (that is, C = C for W<0). The 

summation for the sink/source term allows for 

multiple fluid sinks and sources having 

different associated source concentrations. The 

assumption of fluid incompressibility means 

that all changes in fluid storage are represented 

by changes in porosity in the three-dimensional 

transport equation.
Equation 6 contains velocity divergence 

terms that can be eliminated (Konikow and 

Grove, 1977). This removes several 

derivatives from the transport equation, which 

may reduce errors during the numerical 

solution, as well as removing the direct effect 

of fluid sinks from the governing equation. 

The accumulation and divergence terms in eq. 6 

can be expanded:

dt

= { 

dt dt

(8)

Substituting these expressions into eq. 6, 

adding (WC-WC) = 0, and rearranging terms 

yields

.ac '& 3(pfcC) 

3*

dC

c = 0. (10)

The last term on the left side of eq. 10 contains 

a bracketed term, which is an expression of 

fluid continuity. If fluid continuity is satisfied, 

then the bracketed term is zero, leaving:

+

3(pfo C) __ dC .....r° ' _L v.

d f fD- n *0 v -\

OX;

\ J J

= 0. (11)

This form of the governing equation can be 

designated as "flow-equation-removed" and is 

advantageous, though not required, for the 

method-of-characteristics numerical solution 

because the divergence of velocity does not 

appear. Under these assumptions, including

incompressible fluid, the porosity is not 

constant unless the flow system is in a steadystate equilibrium (Goode, 1990b). Therefore, 

we update the porosity to account for elastic 

changes in aquifer volume caused by transient 

changes in hydraulic head using the method 

proposed by Goode (1990b). Testing and 

evaluation of this approach for incorporation 

into MOC3D is documented by Goode and 

Konikow (1991).

The governing equation can be further 

simplified for the case of reversible, 

instantaneous, equilibrium sorption of the 

solute governed by a linear isotherm. For this 

case, the sorbed concentration, C, is given by:

"r=KjC (\2\ ^ ^d^ \ 1 ^)

where K^ is the sorption coefficient, or 

distribution coefficient, which is assumed to be 

constant in time. The accumulation in the 

sorbed phase can be expressed as:

= K,

dt

(13)

if it is assumed that the aquifer bulk density is 

constant in time. Substituting eqs. 12 and 13 

into eq. 11 gives:

£ I 3f

/

1 +Mi 1 = 0. (14)

The terms controlling sorption can be combined 

into a single parameter the retardation factor 

(/?/), which is defined as:

o =

7

(15)

Rf may be slightly variable in time if the 

porosity changes due to transient flow effects. 

We ignore this possible minor variability and
assume that the retardation factor is constant in 

time. Substituting eq. 15 into eq. 14 yields:

ac y_ dC

dt R dxt dxt J J

eR,

(16)

This is the form of the governing equation 

solved in MOC3D.

We convert eq. 16 from an Eulerian 

framework to a Lagrangian one through the 

material derivative, yielding a simpler form of 

the governing equation (for example, see 

Konikow and Bredehoeft, 1978, p. 6) for the 

concentration of a reference point moving with 

the retarded velocity (V/R/):

(17)

Although this concentration is now that of a 

moving point in space, we retain the same 

symbol, C, as a matter of convenience.

The mathematical solution of the 

governing equations requires the specification 

of certain initial and boundary conditions. 

Because the transport equation is always solved 

for transient conditions, the initial concentration 

must be specified throughout the domain within 

which solute transport occurs (which may be 

equal to or smaller than the domain in which 

the flow equation is applied and solved).

The specification of a constantconcentration boundary condition at one or 

more nodes for the transport equation would be 

analogous to the use of a constant-head 

boundary condition for the flow equation. 

Although this is mathematically and 

numerically feasible, it is rare that a field 

environment would be consistent with such a

constant-concentration condition. Therefore, 

we have not implemented the use of this type of 

boundary condition in this model. Instead, 

input concentrations must always be associated 

with a fluid flux.

For the transport equation, two 

specified mass-flux boundary conditions are 

used in this model. At no-flow boundaries for 

the flow equation, the solute mass flux is also 

required to be zero. The second type of 

specified mass boundary condition is applied 

when the transport subdomain is within a flow 

domain. That is, the boundaries of the 

transport subdomain do not coincide with the 

flow domain boundary. In this case, solute 

mass movement into and out of the transport 

subdomain is assumed to be by advection only; 

no dispersive solute flux can occur across a 

subdomain boundary, which is mathematically 

equivalent to a zero gradient in concentration 

across the boundary.

The effects of all other external fluid 

sinks and sources on transport are incorporated 

through the fluid source/sink terms (W) in eq. 

17. For a fluid source (W>0), denoted W+, 

C' in eq. 17 is the specified source 

concentration of the incoming fluid. For a fluid 

sink (W<0), denoted W_, C in eq. 17 is 

assumed to equal the concentration in the 

aquifer, C, at the location of the sink. In this 

case W(C -C) = (W_)(C- C) = 0. However, 

if the fluid sink is associated with evaporation 

or transpiration, it is assumed that the fluid 

discharge mechanism will exclude dissolved 

chemicals; for this special case C" = 0 and 

W(C -C) = (WL)(0 - C) , which results in an 

increase in concentration at the location of the 

fluid sink.

Dispersion Coefficient

The third term in eq. 16 represents the 

change in concentration due to hydrodynamic 

dispersion. This expression is analogous to
Pick's Law describing diffusive flux. This 

Fickian model assumes that the driving force is 

the concentration gradient and that the 

dispersive flux occurs in a direction from 

higher concentrations towards lower 

concentrations. The dispersion coefficient may 

be related to the velocity of ground-water flow 

and to the nature of the aquifer using 

Scheidegger's (1961) equation:

V m Vn

M

/j,m,n=l,2,3 (18)

where (Xijmn is a component of the dispersivity 

tensor (L), a fourth rank tensor, Vm and Vn are 

components of the velocity vector in the m and 

n directions, respectively, and I VI is the

magnitude of velocity:

\v\ = 'I + Vl + K2 . (19)

Scheidegger (1961) further shows that 

for an isotropic aquifer the dispersivity tensor 

can be defined in terms of the longitudinal and 

transverse dispersivities, (XL and (XT- This 

yields two dispersion coefficients oriented with 

the direction of flow, the longitudinal 

dispersion coefficient, DL, and the transverse 

dispersion coefficient, D? :

DL =aL \V\; = aT \V\. (20)

The directional dispersion coefficients 

in the x, y, and z coordinates are derived by a 

transformation from the L-T coordinates. 

Additionally, an isotropic diffusion coefficient, 

Dm, can be included to account for molecular 

diffusion. The diffusion coefficient includes 

the effects of tortuosity. These manipulations 

yield the components of the dispersion tensor, 

D, which includes diffusion:

(21)

where 8ij = 1 if i = j and Sy = 0 if i ^ j. This 

short-hand notation can be written explicitly as:

YL \v\

Yl

M

YL 

M

v v

\v\ \v\

v v

\v\ T M

(22a)

(22b)

rr + ar-rr + Ai, (22c) 

\V\ T \V\ m

(22d)

V V

(22f)

A number of field studies have 

indicated that transverse spreading in the 

vertical direction is much smaller than 

transverse spreading in the horizontal direction 

(for example, see Robson, 1974; Robson, 

1978; Garabedian and others, 1991; and Gelhar 

and others, 1992). To allow modeling of this 

observed spreading pattern, Burnett and Frind 

(1987) made an ad hoc modification to the 

transverse terms in the preceding expressions, 

which we incorporate into MOC3D:

IL YLM

vl

v2

f Z- \v\

~H YL M

vl

YL 

M

v?

(23a) v /

(23b) v /

v

7

\v\

v-

\v\ -T7 + aT -+- + Dm (23c)

(23d) M
( \ n = n = (a T - ar } yz zy \ L Tv )

VV

(23e)

(23f)

where aT is the horizontal transverse 

dispersivity and aT is the vertical transverse 

dispersivity. If aTfj = aTy , then eq. 23 

reduces to eq. 22.

Review of Assumptions

A number of assumptions have been 

made in the development of the previous forms 

of the governing equations. Following is a list 

of the main assumptions that must be carefully 

evaluated before applying the model to a field 

problem.

1. Darcy's law is valid and hydraulic-head 

gradients are the only significant driving 

mechanism for fluid flow.

2. The hydraulic conductivity of the aquifer 

system is constant with time. Also, if the 

system is anisotropic, it is assumed that the 

principal axes of the hydraulic conductivity 

tensor are aligned with orthogonal frame of 

reference, so that the cross-product terms of the 

hydraulic conductivity tensor are eliminated.

3. Gradients of fluid density, viscosity, and 

temperature do not affect the velocity 

distribution.

4. No chemical reactions occur that affect the 

fluid or aquifer properties.

5. The dispersivity coefficients are constant 

with time, and the aquifer is isotropic with 

respect to longitudinal dispersivity.

As noted by Konikow and Bredehoeft 

(1978), the nature of a specific field problem 

may be such that not all of these underlying 

assumptions are valid. The degree to which 

field conditions deviate from these assumptions 

will affect the applicability and reliability of the

model for that problem. If the deviation from a 

particular assumption is significant, the 

governing equations and the numerical code 

will have to be modified to account for the 

appropriate processes or factors.

NUMERICAL METHODS

Because aquifers are heterogeneous and 

have complex boundary conditions, exact 

analytical solutions to the governing equations 

can not be obtained directly. Instead, 

numerical methods are used, in which the 

continuous variables of the governing 

equations are replaced with discrete variables 

that are defined at grid blocks (or cells or 

nodes). Thus, the continuous differential 

equation, which defines hydraulic head or 

solute concentration everywhere in the system, 

is replaced by a finite number of algebraic 

equations that defines the hydraulic head or 

concentration at specific points. This system of 

algebraic equations generally is solved using 

matrix techniques.

However, numerical methods yield 

only approximate (rather than exact) solutions 

to the governing equation (or equations); they 

require discretization of space and time. The 

variable internal properties, boundaries, and 

stresses of the system are approximated within 

the discretized format. In general, the finer the 

discretization, the closer the numerical solution 

will be to the true solution.

The notation and conventions used in 

this report and in the MOC3D code to describe 

the grid and to reference (or to number) nodes 

are described in figs. 1-2. The indexing 

notation used here is consistent with that used 

in the FORTRAN code for MOD FLOW by 

McDonald and Harbaugh (1988), although not 

necessarily the notation used in the text of their 

report. Our indexing notation maintains 

conformity between the text of this report and 

the FORTRAN code in MOC3D, and the index
order corresponds to an x,y,z sequence, which 

is standard in numerical models. However, 

our notation differs from that used in some 

other ground-water models in that the xdirection is indexed by "j" and increases from 

left to right along a row to indicate the column 

number. Our use of Ax and Ay is synonymous 

with the use of Ar and Ac, respectively, by 

McDonald and Harbaugh (1988). The ydirection is indexed by "i" and increases from 

the top of the grid to the bottom within a 

column to indicate the row number. Thus, in a 

map view of any one horizontal layer, as 

illustrated in fig. 1, the node representing a cell 

in the first row and first column of the grid 

would lie in the upper left corner of the grid. 

The z-direction represents layers and is indexed 

by "k." As indicated in fig. 2, the first layer 

(k = 1) in a multilayer grid would be the top 

(or highest elevation) layer. The saturated 

thickness of a cell (£/,/,&) is equivalent to Az.

Ax:-

U) I

QC

"t

COLUMNS

EXPLANATION

© Node of finite-difference cell

Ax: = Ar: = cell dimension in row direction

Ay i = Aci = cell dimension in column direction

Figure 1. Notation used to label rows, columns, 

and nodes within one layer (/c) of a three-dimensional, block-centered, finite-difference grid for 

MOC3D.

LAYERS

c

r

^

S S S S S S

s s s s

f ~ / / f S

S f S S /

/ / / /" / /

'

'

'

/ 

/ 

/ 

X

s

.=b

Figure 2. Representative three-dimensional grid 

for MOC3D illustrating notation for layers.

Ground-Water Flow Equation

A numerical solution of threedimensional ground-water flow equation is 

obtained by the MODFLOW code using 

implicit (backward-in-time) finite-difference 

methods. The model was coded in FORTRAN 

in a modular style to allow and encourage the 

development of additional packages or modules 

that can be added on or linked to the original 

code. Many such packages or modules, which 

typically allow additional ground-water 

processes, hydrogeological features, solution 

algorithms, or input/output options to be 

represented or used, have been documented 

since MODFLOW was first released. Most of 

these are summarized by Appel and Reilly 

(1994).

MODFLOW is based on use of a blockcentered finite-difference grid that allows 

variable spacing of the grid in three 

dimensions. Flow can be steady or transient. 

Layers can be simulated as confined, 

unconfined, or a combination of both. Aquifer 

properties can vary spatially and hydraulic 

conductivity (or transmissivity) can be 

anisotropic. Flow associated with external 

stresses, such as wells, areally distributed 

recharge, evapotranspiration, drains, and 

streams, can also be simulated through the use 

of specified head, specified flux, or headdependent flux boundary conditions.
MODFLOW offers several options to solve the 

implicit finite-difference equations, including 

the Strongly Implicit Procedure (SIP), SliceSuccessive Overrelaxation (SSOR) methods, or 

Preconditioned Conjugate-Gradient matrix 

solvers (for example, Hill, 1990). Successful 

use of MOC3D, which is programmed as a 

module to MODFLOW, requires a thorough 

familiarity with the use of MODFLOW. 

Comprehensive documentation of MODFLOW 

is presented by McDonald and Harbaugh 

(1988), Harbaugh and McDonald (1996a and 

1996b), and the various reports for additional 

implemented modules.

Average Interstitial Velocity

Because advective transport and 

hydrodynamic dispersion both depend on the 

velocity of ground-water flow, the solution of 

the transport equation requires knowledge of 

the velocity (or specific discharge) field. 

Therefore, after the head distribution has been 

calculated for a given time step or steady-state 

flow condition, the specific discharge across 

every face of each finite-difference cell within 

the transport subgrid is calculated next.

The specific discharge components are 

calculated for each face using a finite-difference 

form of eq. 3. For example, the specific 

discharge in the horizontal plane in the xdirection at the block interface between cells 

j,i,k and 7 + 1, i, fc is (after McDonald and 

Harbaugh, 1988):

= -K. (24)

Ax

where Kxx^+l/2iifk ^ is the interblock hydraulic 

conductivity in the jc-direction on the forward 

face of the cell and the hydraulic gradient is 

based on implicitly calculated heads at the 

adjacent nodes. (Note that the interblock

hydraulic conductivity is commonly defined by 

the harmonic mean, but in MODFLOW the user 

can specify alternative methods for calculating 

the interblock hydraulic conductivity in the 

Block-Centered Flow [BCF] package.)

The seepage velocity at any point within 

a cell must be defined to represent advective 

transport. It is calculated at a point of interest 

within a finite-difference cell based on the 

interpolated estimate of specific discharge at 

that point divided by the effective porosity of 

the cell in which the point is located (see eqs. 

4-5). The interpolation methods are discussed 

later in the section on "Particle Tracking."

Because the velocity is one of the most 

important factors controlling solute transport, it 

is necessary to examine closely the calculated 

velocity field to understand the patterns and 

rates of solute spreading. Therefore, the 

MOC3D model offers the user several options 

to output the velocity data. These options are 

described in Appendix B. The ;c-, y-, and zcomponents of the velocity vector at the nodes 

of the finite-difference cells can be computed 

with equivalent orders of accuracy for all three 

components at a single location. For example, 

the velocity in the jc-direction at node (j,i,k) 

would be computed as

'==^-. (25)

Analogous expressions are used to compute the 

velocity components in the y- and z-directions.

Note that the specific discharge 

components themselves are computed on cell 

faces with a higher order of accuracy than are 

the velocity components at the nodes. For 

example, qx(j+\/2,i,k) * s based on the head 

difference over Ax', whereas Vx at the node is 

based on the head difference over 2Ax. 

However, an estimate of Vy or Vz at the same 

location (j + l/2,i,k) would have to be at a 

lower order of accuracy (for example, 2AxAy 

for Vy). Therefore, the more accurate estimates

10
on cell faces are used in the model as the basis 

for interpolation of particle velocities, but the 

velocities at nodes are written in the output files 

either to enable direct inspection by the user or 

to facilitate postprocessing with visualization 

software.

Solute-Transport Equation

The solute-transport equation is, in 

general, more difficult to solve accurately using 

numerical methods than is the ground-water 

flow equation, largely because the mathematical 

properties of the transport equation vary 

depending upon which terms in the equation 

are dominant in a particular system. Where 

solute transport is dominated by advective 

transport, as is common in many field 

problems, then the transport equation 

approximates a hyperbolic type of equation 

(similar to equations describing the propagation 

of a wave or of a shock front). In contrast, 

where a system is dominated by dispersive 

fluxes, such as might occur where fluid 

velocities are relatively low and aquifer 

dispersivities are relatively high, then the 

transport equation becomes more parabolic in 

nature (similar to the transient ground-water 

flow equation). To further complicate matters, 

because system properties and fluid velocity 

may vary significantly, the dominant process 

(and the mathematical properties of the 

governing equation) may vary from point to 

point and over time within the same domain.

Dispersion-dominated solute-transport 

problems (and parabolic equations in general) 

are quite amenable to accurate and efficient 

solution using standard finite-difference and 

finite-element numerical methods. However, 

in solving advection-dominated transport 

problems, in which a relatively sharp front (or 

steep concentration gradient) is moving through 

a system, it is numerically difficult to preserve 

the sharpness of the front. In such cases, 

numerical solutions often will include either

erroneous oscillations (overshoot and 

undershoot) or calculate a greater dispersive 

flux than would occur by physical dispersion 

alone or than would be indicated by an exact 

solution of the governing equation (for 

example, see Finder and Gray, 1977). That 

part of the calculated dispersion introduced 

solely by the numerical solution algorithm is 

called numerical dispersion.

Method of Characteristics

The method of characteristics was 

developed to solve hyperbolic differential 

equations. A major advantage is that the 

method minimizes numerical dispersion (or 

even eliminates it in limited cases) (Garder and 

others, 1964; Finder and Cooper, 1970; 

Reddell and Sunada, 1970; Bredehoeft and 

Finder, 1973; Konikow and Bredehoeft, 1978; 

Zheng, 1990). The approach taken by the 

method of characteristics is not to solve eq. 16 

directly, but rather to solve an equivalent 

system of ordinary differential equations. 

Equation 16 can be rearranged to obtain:

d£

dt

1 d

eR<

dc

-AC.

Vi dC 

Rf dxt

(26)

Equation 26 describes the change in 

concentration over time at fixed reference 

points within a stationary coordinate system, 

which is referred to as an Eulerian framework. 

An alternative perspective is to consider 

changes in concentration over time in 

representative fluid parcels as they move with 

the flow of the fluid past fixed points in space. 

This, in effect, is a moving coordinate system, 

which is referred to as a Lagrangian 

framework. We convert eq. 16 from an 

Eulerian framework to a Lagrangian one 

(essentially, a framework of a moving grid)

11
through the material derivative. The material 

derivative of concentration with respect to time, 

dC/dt, describes the change in concentration in 

a parcel of water moving at the average 

interstitial velocity of water; it may be defined 

for a three-dimensional system as:

dC_dC dCdx dCdy dC dz ~^~~~^~ + ~^ "i""1""^ T""1"^! 7~' 

dt dt dx dt dy dt dz dt

The last three terms on the right side include the 

material derivatives of position, which are 

defined by the velocity in the x-, y-, and zdirections. We then have:

dt Rf

dt Rf

dt R

(28)

(29)

(30)

f

and a simpler form of the governing equation 

(for example, see Konikow and Bredehoeft, 

1978, p. 6) for the concentration of a reference 

point moving with the retarded velocity (V/R/) 

is obtained by substituting the right sides of 

eqs. 26 and 28-30 for the corresponding terms 

in eq. 27:

dC 1 - = - d £Dn ij-^ dC

dt £Rf dX; \ dX;

eR<

(31)

Although this concentration is now that of a 

moving point in space, we retain the same 

symbol, C, as a matter of convenience.

The solutions of the system of 

equations comprising eqs. 28-31 may be given 

as x = x(t), y = y(t), z = z(t), and C = C(t), 

and are called the characteristic curves of eq. 

26. Given solutions to eqs. 28-31, a solution 

to the original partial differential equation may

be obtained by following the characteristic 

curves. This may be accomplished by 

introducing a set of moving points (or reference 

particles) that can be traced within the 

stationary coordinates of a finite-difference 

grid. Each particle corresponds to one 

characteristic curve, and values of x, y, z, and 

C are obtained as functions of t for each 

characteristic (Carder and others, 1964). Each 

point has a concentration and position 

associated with it and is moved through the 

flow field in proportion to the flow velocity at 

its location (see fig. 3). Equation 31 can be 

solved using any one of several approaches, 

including random-walk methods (for example, 

see Prickett and others, 1981; and Engineering 

Technologies Associates, Inc., 1989). We

EXPLANATION

Flow line and direction of flow 

 Computed path of particle 

Node of finite-difference cell 

New location of particle 

Initial location of particle

Figure 3. Part of a hypothetical finite-difference grid showing relation of flow field to 

movement of points (or particles) in methodof-characteristics model for simulating solute 

transport (modified from Konikow and Bredehoeft, 1978).

12
choose to solve eq. 31 using explicit finitedifference approaches because of its efficiency 

and relative simplicity. However, as noted by 

Reddell and Sunada (1970) and Konikow and 

Bredehoeft (1978), this approach also requires 

the consideration of stability criteria, which will 

be discussed in the section "Stability and 

Accuracy Criteria." In some cases, the 

efficiency may be adversely affected because of 

restrictions imposed by the stability criteria.

As noted by Konikow and Bredehoeft 

(1978), the processes of advection, dispersion, 

mixing, and reactions are occurring 

continuously and simultaneously. Therefore, 

eqs. 28-31 should be solved simultaneously, 

but for practical reasons, they are solved 

sequentially. However, the results will be 

sensitive to the order in which they are solved 

because the change in concentration due to 

dispersion depends on the concentration 

gradient, and the concentration gradient at any

location may change significantly because of 

advection during a time increment. This is 

illustrated in fig. 4. As the position of a 

concentration front or breakthrough curve 

advances with time, say from time t at the 

beginning of a time increment, to time t+l at 

the end of a time increment, the concentration 

gradient at any particular location or reference 

point is continuously changing. For example, 

in fig. 4 the concentration gradient at the 

location where the relative distance equals 0.4 

is relatively steep at time f, as indicated by the 

slope of the curve labeled O. However, if the 

solute migrated in the direction of flow to a 

new location due to advection only without 

being subjected to any dispersive flux, as 

indicated by the curve labeled Ct+adv, then the 

concentration gradient at the same point in 

space would be much smaller. If the change in 

concentration caused by dispersion were 

computed by solving eq. 31 first in the

0.0 0.4 0.6 

RELATIVE DISTANCE

0.8 1.0

Figure 4. Representative change in breakthrough curve from time level t to f+1. Note that 

concentration changes are exaggerated to help illustrate the sequence of calculations. Curve 

for Ct+adv represents the concentration distribution at time f+1 due to advection only. 

(Modified from Konikow and Bredehoeft, 1978.)

13
sequence and using the spatial concentration 

gradients at the beginning of a time step, the net 

results would be different from those that 

would be computed from solving eq. 31 last in 

sequence and using concentration gradients 

after advection. Sensitivity to the sequence of 

solving the characteristic equations is 

eliminated by solving eq. 31 using 

concentration gradients based on the average of 

the concentrations at each node before and after 

advection. This effectively gives equal weight 

to the concentration gradients before and after 

advection while computing the solute flux due 

to dispersion. The averaged concentrations, 

designated as Cy- 1- ^, are calculated as:

-; ; f, J,l,K (32)

where Cj£jfv is the concentration at the new 

time level after advection alone.

Particle Tracking

Advection in flowing ground water is 

simulated by particle tracking. The other terms 

in the governing equation dispersion, 

sources, and decay are accounted for by 

adjusting the concentrations associated with 

each particle. The concentration changes 

caused by dispersion and fluid sources are 

computed in the fixed-in-space finite-difference 

grid, whereas concentration changes caused by 

decay are calculated directly on the particles, as 

described below. Initial particle locations are 

defined through model input, and subsequent 

particle positions are integrated in time using 

the spatially and temporally variable velocity 

field.

The fluid velocity is derived from a 

block-centered finite-difference solution of the 

three-dimensional flow equation (McDonald 

and Harbaugh, 1988). The components of the 

velocity vector are approximated by difference 

expressions at interfaces between adjacent

blocks occupying three-dimensional space (see 

eq. 24). In the quasi-three-dimensional 

approach taken by McDonald and Harbaugh 

(1988), horizontal fluxes are computed using 

the transmissivity of each layer and the vertical 

fluxes are computed using the vertical 

conductance, which is the vertical hydraulic 

conductivity divided by the vertical distance 

between layer centers. Thus, the components 

of the specific discharge vector can be 

represented by:

-T

Ax

 T

"

yy(j,i+l/2,k) Ay

(33a)

(33b)

(33c)

Note that in eqs. 33a-b, a value of b is 

not explicitly calculated on the face of a cell. 

Instead, we calculate a value for the product bq 

on the face using a mean transmissivity. 

Having computed specific discharge 

components at block interfaces, spatial 

interpolation is used to estimate the specific 

discharge and velocity at the locations of all 

particles (using values of b and £ for the 

particular cell in which the particle of interest is 

located).

MOC3D uses either linear or bilinear 

spatial interpolation, as specified by the user. 

These two alternatives are described below in 

the sections on interpolation. Depending on the 

velocity interpolation method chosen, changes 

in particle position are computed either 

explicitly for bilinear interpolation, or semianalytically for linear interpolation. The 

integration of particle position in time is also 

described in the sections on interpolation.

14
A moving point in a ground-water flow 

system will change velocity as it moves due to 

both spatial variation in velocity and temporal 

variations during transient flow. During a flow 

time step, advection is determined from 

velocities computed at the end of the flow time 

step. Temporal changes in velocity are 

accounted for by a step change in velocity at the 

start of each new flow time step. After the 

flow equation is solved for a new time step, 

specific discharges are recomputed on the basis 

of the new head distribution, and the movement 

of particles during this flow time step is based 

only on these specific discharges.

Linear interpolation and semi-analytic 

integration

Linear interpolation of specific 

discharge is formally consistent with the blockcentered finite-difference flow solution in that 

the governing flow equation is satisfied locally 

within each cell (Goode, 1990a). In this case, 

each velocity component is linearly interpolated 

in the direction of the component of interest. 

The solute velocity at any particular point or 

location within a cell, which may be retarded 

with respect to the average interstitial velocity 

of water, is:

-1

where &c, Sy, and 8z are the spatial weights for 

interpolation. The spatial weights are given by:

Ax

oy s =

Ay

Az

(35a)

(35b)

(35c)

where the subscript p indicates the particle 

number and Sx, Sy, and 8z can range from 0.0 

to 1.0.

Vx at p is calculated from eqs. 34a and 

3 5a. If a layer has a uniform porosity and 

thickness, and the solute were nonreactive, 

then V^ at p can be calculated directly from the 

following simplification of eq. 34a as:

(34c)

_}/2M} . (36)

The linear interpolation scheme is 

illustrated for Vx in fig. 5, which shows a 

particle p located in cell y,/,&. Vx at p depends 

on the relative position within the cell in the redirection only, and does not vary with changes 

in particle position in the y- or z-directions 

within the cell.

For simplicity, a dimensionless local 

coordinate is used to facilitate particle tracking. 

Internally to the program, the particle 

coordinates are represented, for example, in the 

^-direction by:

xp =(j + xp )Ax (37)

where j is the column index of the block where 

a particle is located and -0.5<Jcp<0.5 is the 

particle position in the ^-direction relative to the

15
J

^

I

',i,k)

EXPLANATION

<~t»* Component of velocity vector 

© Node of finite-difference cell 

o Initial location of particle

Figure 5. Spatial weights used in linear interpolation method to estimate Vx at location of a 

particle in cell j,i,k.

node (that is, xp =(xp /Ax}- y). The value 

j+ xp is stored in the model and, for example, 

ranges from 4.5 to 5.5 for particles within 

column 5. Thus, the spatial weights can be

represented in terms of the particle positions 

relative to the node, and the velocity at the 

location of a particle is estimated using linear 

interpolation as:

[(Rf } k ( £b^j,i, 

) = [(Rf)k ( £)j,i,

-1

Assuming that the velocity is constant 

in time during a time increment for solving the 

transport equation, the particle position can be 

integrated analytically within each cell (Goode, 

1990a). The integration of particle position is 

accomplished in a multi-step procedure. First, 

the velocity components in each direction at the 

starting position of the particle are computed 

using eqs. 38a-c. Then the gradient of velocity 

within the cell for each direction is computed 

from the last term of eqs. 38a-c and

(38a> 

(38b)

denominator terms appropriate for the particular 

direction. For example, the velocity gradient in 

the ^-direction, gx, is calculated by dividing the 

last term of eq. 38a by the product of the 

retardation factor, porosity, saturated 

thickness, and grid spacing in the ^-direction, 

as:

16
It is implicit in the linear interpolation approach 

that

dy dz

= n0. (40)

Given the solutions to eqs. 38a-c and 

39, the time of travel to each cell boundary, in 

the x-, y-, and z-directions, can be computed 

analytically. For example, the time for a 

particle to reach or cross the ^-boundary, Atx , 

is given by:

8

In

V.x,0

(41)

where V^o is the *-velocity at the starting 

location of the particle, *o is the initial position 

of the particle, and xe is the ^-coordinate of the 

exit face of the cell. By definition, if the 

starting velocity is positive, then the only 

possible ^-boundary for the particle to exit is 

that in the positive ^-direction. If gx is zero, or 

practically zero (that is, less than 10'20), then 

Vx is assumed to be uniform in space and 

constant during the particle motion, and the 

travel time to the exit boundary can be 

computed directly from:

"V .__ "V

(42)

Vx,0

The travel time for a particle to exit the 

cell, Ate , is the minimum of Atx , Aty , and 

Atz . If the cell exit travel time ( Ate ) is greater 

than the remaining time for this time increment, 

then the ending particle position is within the 

same cell that it started during that time 

increment, and the particle coordinates at the 

end of the time increment are determined from 

the velocity at the initial location, the rate of 

change in velocity in the direction of 

movement, and the length of that time 

increment. For example, the new ^-coordinate 

at the end of the move would be:

where the subscript "1" denotes the end of the 

particle movement for the time increment. If 

gx , the gradient in ^-velocity, is zero or 

practically zero (again, less than 10~20), then 

the ending coordinate is simply:

On the other hand, if the time to reach any 

cell boundary is less than the remaining time 

for this time increment (At), then the particle 

exits the cell during the time increment. In this 

case, the particle exits the cell on the boundary 

for the direction having the minimum travel 

time (Atx , Atyt or Atz ). For example, if Atx 

is less than Aty , Atz , and At, then the particle 

leaves the cell by crossing the ;c-cell boundary 

after the "cell exit travel time" ( Ate = Atx ) has 

elapsed.

The particle position must be evaluated 

when it leaves a cell during a time increment. 

By definition, the position in the direction of 

the exit face is known a priori. The other two 

coordinates are computed by inserting the cell 

exit travel time, Ate , into eq. 43 instead of the 

entire time increment, At. For example, for a 

particle that leaves a cell on the Xj+2/2 face, the 

y- and z-coordinates of that particle where it 

crosses the cell face are given by:

y\ = (45a)

(45b)

If the gradient terms in eq. 45a or eq. 45b are 

zero, or essentially zero, then the y and z 

equivalents of eq. 44 are used with At = Ate .

Also, if Ate < At, the time remaining in 

the time increment, At' , is the original time 

increment minus the cell exit travel time:

A? = At-At.. (46)

V,

Xl =XQ

ft

\\ (43) Subsequent particle movement is computed as

17
outlined above, but with a reduced time interval 

( At'} for its path in the new cell, until the entire 

original time increment (At) is exhausted. Note 

that in some cases, such as when a particle is 

located near a corner of a cell at the start of a 

move, it is possible that during that particular 

time increment this sequence of calculations 

may have to be repeated a third time before the 

move is completed for that particle.

Note that linear interpolation yields a 

discontinuous velocity field (Goode, 1990a). 

Equations 39 and 40 show that velocity 

components vary as a function of distance only 

in the direction of that particular component of 

the velocity vector. Thus, when linear 

interpolation is used, specific discharge and 

velocity components change abruptly at block 

or cell interfaces parallel to the direction 

indicated by that component, even for 

homogeneous aquifers. For example, the 

component of flux in the jc-direction, qx, is 

continuous at cell interfaces in the jc-direction 

(at j-l/2,i,k andy+l/2,/,fc), but it can change 

abruptly at the other four cell faces that are 

perpendicular to the y and z axes, even when K 

is uniform.

Bilinear interpolation and explicit 

integration

In a homogeneous aquifer the velocity field 

would be continuous and smoothly varying, 

except at the locations of strong sources and 

sinks, unlike the discontinuous velocity field 

calculated using linear interpolation. Bilinear 

interpolation generates such a continuous 

specific discharge field, as shown by Goode 

(1990a) for a case of two-dimensional flow. 

MOC3D includes an option to use bilinear 

interpolation in the x-y (horizontal) plane to 

determine the jc- and y-velocity components, 

which then will be continuous and smoothly 

varying in the x-y plane. Because of 

stratigraphic layering, many aquifer systems 

are more heterogeneous in the vertical direction 

than within a horizontal plane. Hence, specific 

discharge and velocity might be more often 

expected to change abruptly (rather than 

smoothly) in the vertical direction. Therefore, 

MOC3D always uses linear interpolation for 

calculating the z-component of particle velocity.

Following Konikow and Bredehoeft (1978) 

and Goode (1990a), the jc- and ^-velocity 

components are bilinearly interpolated in the 

horizontal plane by:

;,* ,*] {I1 ~ 

i/2,f*,*) + (x

and

Vy(P) =

-1

where the local particle coordinates (xp and 

yp } are as defined previously in eq. 37, and the 

interpolation factors (F,-) and adjacent node 

locators 0'* and /*) are defined as:

(47)

(48)

Fy = 0.0

Fy=yP /* = /+!

for yp <0 (49a) 

for yp =0 (49b) 

for yp >0 (49c)

18
and

FX=-XP j* = 7-1

Fx = 0.0

= xr j* = 7+1

for xp < 0 (49d) 

for jcp = 0 (49e) 

for jL>0.(49f)

The bilinear interpolation scheme to 

estimate V* and V-y is illustrated in fig. 6, 

which shows a particle /? located in cell j,i,k. 

Vx and Vy at /? depend on the relative position 

within the cell in both the jc- and y-directions. 

Vx and Vy at/? are calculated from eqs. 47 and 

48 using the weighting factors defined in eqs. 

49a-f. Note that ( xp + 0.5) of eq. 47 equals Sx 

in fig. 6, and (y + 0.5) of eq. 48 equals 8y in 

fig. 6. Note that the interpolated values of Vx 

and Vy at p are determined using specific 

discharge values at differing adjacent cell faces.

Also, if the cell face adjacent to the

quadrant containing a particle of interest 

represents a no-flow boundary, then the 

component of the flow parallel to that boundary 

should not change in the area of the half-cell 

between the node and that face. Therefore, in 

such a case, linear interpolation is used within 

that half-cell area for that component only. For 

example, in fig. 6, if cell (/,/+!,£) were a noflow cell, then Vx at p would be estimated by 

linear interpolation between (j-l/2,i,k) and 

(j+l/2,i,k). Vy at/? would still be estimated 

using bilinear interpolation, as shown in fig. 6, 

properly recognizing that Vy = 0.0 at 

(M+l/2,*).

Bilinear interpolation yields velocity 

components that depend on both the x- and ypositions, hence the semi-analytic integration of 

position used above is not applicable. For 

bilinear interpolation, particle movement is 

computed explicitly using the velocity of the

EXPLANATION

Area of influence for 

interpolating Vx at particle p

Area of influence for 

interpolating Vy at particle p

Component of velocity vector

o Node of finite-difference cell 

o Initial location of particle

Figure 6. Interpolation factors used in bilinear interpolation method to estimate horizontal 

components of velocity, Vx ar\6 Vy, at the position of a particle located in the southeast quadrant of 

cell j,i,k.

19
particle at its starting position and the length of 

the time increment:

Vx>0At (50a)

(50b)

The times required for a given particle to reach 

the jt- and v-faces of the cell in which it is 

located are computed using:

4f,=

V.;c,0

and

V,

(51)

(52)

For the z-direction, all factors are computed as 

described previously in eqs. 34-45.

Also, following the same procedure as 

implemented in the previously described linear 

interpolation scheme, a particle that would 

cross a cell boundary during a move interval 

(or transport time increment) is temporarily 

stopped at whichever cell boundary (or face) it 

reaches first. At that location the particle 

velocity is updated using the porosity and 

thickness of the cell that the particle is entering, 

and the particle is then moved an appropriate 

distance in the new cell for the remainder of the 

time increment (At', as calculated from eq. 

46), or until it reaches another cell face. 

Computation of particle movement continues 

until the entire movement time interval is 

exhausted.

Discussion Choosing an appropriate 

interpolation scheme

As noted by Goode (1990a), selecting the 

best interpolation scheme to determine particle 

velocity in a ground-water flow model depends 

in part on the conceptualization of aquifer 

heterogeneity. Linear interpolation is directly 

consistent with the block-centered finitedifference solution of the flow equation. 

However, linear interpolation produces a 

discontinuous velocity field, even for 

homogeneous media (Goode, 1990a). In the 

presence of strong heterogeneities between 

adjacent cells within a layer, it would usually 

be preferable to select the linear interpolation 

scheme.

If transmissivity within a layer is 

homogeneous or smoothly varying, bilinear 

interpolation of velocity yields more realistic 

pathlines for a given discretization than does 

linear interpolation. In such cases, the bilinear 

interpolation scheme may be preferable because 

it will yield pathlines more consistent with the 

conceptualization. For example, Goode 

(1990a) shows that bilinear interpolation yields 

more accurate travel times and pathlines for a 

case in which the interblock hydraulic 

conductivity is computed using the logarithmic 

mean (Goode and Appel, 1992) corresponding 

to a linear variation of hydraulic conductivity 

between nodes. Because vertical heterogeneity 

is significant in typical aquifer systems, the 

MOC3D code always uses linear interpolation 

for velocity in the vertical direction.

In general, the choice of interpolation 

scheme will only make a small difference in the 

final solution. Using a finer grid will 

invariably have a larger effect and yield a more 

accurate definition of the flow field, particularly 

in areas where the hydraulic properties are 

changing in space.

Decay

Decay is simulated by reducing the 

particle concentrations after advection. At this 

point the particle concentration has not yet been 

adjusted for dispersion and sources. However, 

the change in particle position accounts for 

advection up to time increment t+l. The loss 

of solute mass during a given transport time 

increment (At) because of decay processes is

20
accounted for by computing the decayed 

particle concentration, Cp :

f-,d _ (~<t -kAt x*-o\ 

^P ~ ^P ? V3)

where Cp is the particle concentration at the 

start of the move interval (and before 

advection).

A major advantage of calculating the effect 

of decay directly on the particles, rather than on 

the nodal concentrations, is that this procedure 

eliminates any possible reduction in precision 

(or numerical dispersion) caused by the 

framework transfer between the moving grid 

and the fixed grid (that is, from particle 

concentrations to average nodal concentrations 

and back again to particle concentrations), as is 

done in computing changes in concentration 

due to dispersion. The dispersion calculation is 

then based on average concentration gradients 

during a time increment, an average that is 

based on the advected and decayed 

concentrations at the end of the time increment 

and the particle concentrations at the start of the 

time increment.

As noted by Goode and Konikow (1989), 

this exponential formulation has no numerical 

stability restrictions associated with it. 

However, if the half-life is on the order of the 

transport time increment or smaller, then some 

accuracy will be lost because of the explicit 

decoupling of decay and other transport 

processes.

When a solute subject to decay enters the 

aquifer through a fluid source, it is assumed 

that the fluid source contains the solute in the 

concentration specified by C". The governing 

equation and the MOC3D model assumes that 

decay only occurs within the ground-water 

system, and not within the source reservoir. In 

other words, for a given stress period, the 

model assumes that C remains constant in

time and does not decay. If the problem being 

simulated requires that the source fluid itself 

undergo decay, then the code will have to be 

modified to allow this.

Node Concentrations

After all particles have been moved, the 

concentration at each node is temporarily 

assigned the average concentration of all 

particles then located within the volume of that 

cell; this average concentration is denoted as

(54) iadv _ p \

I 

p=l

=J,l ,' ; = =k\7 \ 

where the S function is 1 if the particle is within 

the cel\j,i,k, and is zero otherwise. The time 

index is labeled "adv" because this temporarily 

assigned average concentration represents the 

new time level only with respect to advective 

transport and decay. With respect to the finitedifference grid, the effect of advective transport 

is to move particles with differing 

concentrations into and out of each cell.

Finite-Difference Approximations

The divergence of dispersive flux is 

normalized by the retardation factor and 

porosity to yield the rate of change in 

concentration. In addition, in a quasi-3D 

approach, changes in saturated thickness are 

incorporated for horizontal flux terms. The rate 

of change in concentration due to dispersion 

and mixing in cells having a fluid source can be 

written:

21
(55)

where subscript m is a summation index for the 

dispersion term. The//, k subscripts in eq. 55 

denote the spatial finite-difference grid 

indexing, as discussed previously in the section 

"Numerical Methods." The superscript "*" 

indicates that the terms depend on the average 

of the concentration at the old time level, and 

the concentration at the new time level after 

advection (see eq. 32). These averaged 

concentrations are used to calculate the solute 

flux terms indicated by the superscript "*".

The components of the dispersive flux in 

each direction across cell faces are calculated 

using finite-difference approximations that are 

centered-in-space and explicit (forward-intime). A detailed description of these finitedifference approximations is given in Appendix 

A.

The explicit finite-difference approach is 

conceptually straightforward, but only 

applicable for certain conditions. Reddell and 

Sunada (1970, p. 62) show that for an explicit 

finite-difference solution of eq. 31 to be stable, 

the following constraint must be met:

i +

D»At . 1

<f. (56)

(AxY (AyY (AzY 2

Inspection of eq. 56 shows that the constraint 

will be most readily met for relatively small 

values of the dispersion coefficient (that is, for 

advection-dominated problems). Note that 

stability does not necessarily assure accuracy. 

If the constraint is not met for a given set of 

physical parameters, then either the grid

spacing must be increased (with a consequent 

loss of accuracy) or the time increments must 

be decreased until the condition is satisfied.

Note that the dispersive fluxes calculated by 

solving the eqs. A1-A3 include contributions 

related to the cross-product terms of the 

dispersion tensor. However, the constraint 

expressed in eq. 56 only includes the principal 

components of the dispersion tensor. In 

circumstances where the cross-product terms 

are relatively large, this can lead to a calculation 

of a negative concentration at a node if the 

calculated solute flux out of a cell during a time 

increment is greater than the solute mass in that 

cell at the end of the previous time increment. 

Considering eqs. 22a and 22d, for example, 

we can see that the cross-product terms (such 

as D and D ) would be large relative to a 

principal diagonal coefficient (such as D^) 

when both Vx and Vy are significant and aL is 

significantly greater than aT . For a given 

velocity and fixed dispersivities, the crossproduct terms are maximized when Vx \ = Vy 

and flow is at 45 degrees to the grid. Negative 

concentrations are more likely to occur using 

this formulation in the presence of steep 

concentration gradients.

Consider the solute mass flux due to 

dispersion across one face of a cell, say in the 

jc-direction across the cell face at (j+l/2,i,k). 

The change in solute mass, Mffj+m^^k), is 

equal to the rate of dispersive flux, expressed 

in eq. Al, multiplied by the length of the time 

increment and by the width of the cell:

22
i v +1

(DXX ) \ xxjj+i/2tit Ax

 te /--(* /--I* /--I* yo,* ~/.i'+U + S'+U+U ~ ~/,i-U ~ Cj+l,i-l

4Ay

,* A

j+l/2,i,k 2

This mass can be compared with the solute 

mass in the cell at the start of the time 

increment, Mj^, which is given by:

M^^AxAyebC'^. (58)

The criterion in eq. 56 is equivalent to 

requiring that the sum of the solute mass fluxes 

across all faces of a cell must be less than or 

equal to Mjti>k . However, because the 

criterion does not include the dispersive flux 

related to the cross-product terms, it can lead to 

oscillations in the solution that yield negative 

concentrations, although it would rarely lead to 

a strict stability problem. Experience shows 

that when oscillations do occur because of the 

cross-product terms, they are usually small and 

tend to damp out over time.

An ad hoc procedure was formulated and 

implemented in MOC3D to help minimize and 

limit the occurrence of negative concentrations 

due to the cross-product dispersive flux. The 

approach is to limit the mass flux during one 

time increment across each cell face, as 

calculated by eq. 57, to the mass available in 

the cell, as calculated by eq. 58. The flux 

across each cell face is checked independently, 

so this constraint will not completely eliminate 

negative concentrations, but our experience 

indicates that it will often reduce their 

occurrence significantly. Because the same 

check will be applied from both sides of a 

given face, if the constraint is applied, it will be 

applied equally from both adjacent nodes, so 

the procedure will not affect the global mass

2BM* ;

(57)

balance. That is, the solute mass flux into one 

cell always corresponds to the mass flux out of 

the adjacent cell, whether or not this limiting 

procedure is implemented for that particular cell 

face.

Stability and Accuracy Criteria

As noted by Konikow and Bredehoeft 

(1978), the explicit numerical solution of the 

solute-transport equation has a number of 

stability criteria associated with it. These may 

require that the time step used to solve the flow 

equation be subdivided into a number of 

smaller time increments to accurately solve the 

solute-transport equation.

First, consider the explicit finite-difference 

solution to calculate changes in concentration 

due to dispersion. Solving eq. 56 for At, and 

accounting for the effects of retardation, we see 

that

At < Mm (59)

(over 

grid)

0.5

Dxx

Rf(Ax)2 ' 1

Dyy , D^

lf(Ay)2 Rf (Az)2

Because the solution to eq. 31 is actually 

written as a set of N equations for N nodes, the 

maximum permissible time increment is the 

smallest At computed for any individual node 

in the entire transport grid. The smallest At 

will then occur at the node having the largest 

value of

23
DXX

Dyy D.zz

Rf(AxY

Next consider the effects of mixing ground 

water of one concentration with injected or 

recharged water of a different concentration, as 

represented by the source terms in eq. 31. The 

change in concentration in a source node cannot 

exceed the difference between the source 

concentration (C'j>i>k ) and the concentration in 

the aquifer ( Cj>i>k ), and the maximum possible 

change occurs when a source completely 

flushes out the volume of water in an aquifer 

cell at the start of a time increment. Konikow 

and Bredehoeft (1978) show that this condition 

translates to

(60)

Solving eq. 60 for At at all nodes yields the 

following criterion:

At < Min

(over grid)

eR,

(61)

A third type of criterion involves the 

movement of points to simulate advective 

transport. The distance a particle moves during 

a time increment is equal to (or approximately 

so in cases where particles cross a cell face and 

the adjacent cells have different properties) the 

velocity at the location of the particle times the 

length of the time increment. In effect, this 

constitutes a linear spatial extrapolation of the 

position of a particle from one time increment 

to the next. Konikow and Bredehoeft (1978) 

note that where streamlines are curvilinear, the 

extrapolated position of a particle will deviate 

from the streamline on which it was previously 

located. This deviation introduces an error into 

the numerical solution that is proportional to 

At. An accurate computation of concentration 

changes caused by advective transport requires 

the maintenance of a relatively uniformly 

spaced field of marker particles that are moving 

along relatively smooth and continuous

pathlines. The degree of curvilinearity of 

streamlines in the calculated head field is 

constrained by the grid spacing, as the finitedifference solution to the flow equation 

inherently assumes linear variations in head 

between adjacent nodes. Also, if the distance a 

particle moves in any direction during one time 

increment is greater than the grid spacing in that 

direction, it might be possible for a particle to 

cross a no-flow boundary (or even leave the 

model domain) during one time increment. 

Thus, for a given velocity field and grid, some 

restriction must be placed on the size of the 

time increment to assure that the distance a 

particle moves in the x-, y-, or z-directions 

during one time increment does not exceed 

some critical distances, which can be related to 

the grid spacing in each direction.

These critical distances can be related to the 

grid dimensions by

& Vx(p) < jAx 

At V < jAy

and

At V.z(p)

(62a)

(62b)

(62c)

where y is the fraction of the grid dimensions 

that particles will be allowed to move (normally, 0</<1). Note that these accuracy 

criteria are equivalent to requiring that the 

Courant number be less than or equal to 1. 

However, the model is designed to allow the 

user to specify the value of /(named CELDIS 

in the code and input instructions).

Because these criteria are governed by the 

maximum velocities in the system, and since 

the computed velocity of a tracer particle will 

always be less than or equal to the maximum 

velocity components computed at cell 

boundaries, we have to check only the latter. 

Substituting grid velocity components and 

solving eq. 62 for At results in

24
At< (63a)

max

At< , jAy (63b)

max

and

At< (63c)

max

If the time step used to solve the flow 

equation exceeds the smallest of the time limits 

determined by eqs. 59, 61, or 63, then the time 

step will be subdivided into an appropriate 

number of equal-sized smaller time increments 

to solve the solute-transport equation so that 

none of these limits are exceeded. To help the 

user analyze the results and the grid design, the 

model output will include a statement clarifying 

which of the several criteria were limiting and 

at which node the limiting condition occurred.

Mass Balance

Mass balance calculations are performed to 

help check the numerical accuracy and 

precision of the solution. The principle of 

conservation of mass requires that the net mass 

flux (cumulative sum of mass inflows and 

outflows plus any mass lost or removed by 

reactions) must equal the mass accumulation 

(or change in mass stored). The difference 

between the net flux and the mass accumulation 

is the mass residual (Rm) and is one measure of 

the numerical accuracy of the solution. 

Although a small residual does not prove that 

the numerical solution is accurate, a large error 

in the mass balance is undesirable and may 

indicate the presence of a significant error in the 

numerical solution (Konikow and Bredehoeft, 

1978).

As part of the mass balance calculations, 

the solute fluxes contributed by each distinct 

hydrologic component of the flow and

transport model are accumulated and itemized 

separately to produce a solute budget for the 

system being modeled. The budget is a 

valuable assessment tool because it provides a 

measure of the relative importance of each 

component to the total solute budget. The 

budgets should always be reviewed for 

consistency with the conceptual model and as a 

"reality check" on the model calculations.

In the method of characteristics, the 

accuracy of the solution is associated with the 

concentrations being tracked on the particles. 

However, it is computationally difficult to 

compute a mass stored in the system directly 

from the particle concentrations because their 

relative positions are constantly changing and 

they do not explicitly track a solute or fluid 

mass. Therefore, the mass in storage at any 

time is calculated from the concentrations at the 

nodes of the transport subgrid of the finitedifference mesh. In that sense, the calculated 

and printed mass balance values are themselves 

only an approximation.

The mass residual is computed as

Rm =AMs -Mf (64)

where AMS is the change in mass stored in the 

aquifer, and My is the net mass flux.

The above two mass balance terms are 

evaluated using the following equations:

,k) (65)

where Cjik is the initial concentration at node 

(/.*,&)» M/L3 , and C" /it is the concentration at 

that node at the end of the time increment; and

Mf = j

k i j n

(66)

For cases where W in eq. 66 represents a fluid 

source, C is the specified source concentration 

for that node. Where W represents a fluid 

sink, C" is assumed to equal the average 

concentration in cell (j,i,K) at the beginning of 

that time increment (or move interval) for 

solving the transport equation.

25
The percent error (E) in the mass balance is 

calculated by relating the residual to an 

appropriate measure of the solute flux or mass 

accumulation in the system. However, the 

appropriate basis (or denominator) is problem 

dependent. Thus, the model compares the 

residual with the cumulative mass inflow to the 

system (Mi) or mass outflow from the system 

(M0), whichever is dominant in a particular 

problem. However, if the solute mass entering 

or discharging the system with fluid sources 

and sinks is zero or very small, as it would be, 

for example, in a problem simulating the 

movement of an initial slug within the system, 

then comparing the residual to the mass flux in 

or out of the system could indicate a very large 

error when the numerical solution is actually 

quite accurate. Therefore, in these cases, the 

error will be computed by comparing the 

residual with the mass of solute stored in the 

aquifer, as described by Konikow and 

Bredehoeft (1978). The model will calculate, 

print, and label whichever of the following 

three measures of error are appropriate for the 

problem being simulated.

100- OR- (67a) E,=

M,

_ 100.0Ra 

100.0Rm

(67b) 

(67c)

ES is calculated only if the mass flux in or out 

of the system is less than 50 percent of the 

initial mass stored.

Errors in the mass balance for the flow 

model should generally be less than 0.1 

percent. However, because the solutetransport equation is more difficult to solve 

numerically than is the flow equation, the 

mass-balance error for a solute is often greater 

than for the fluid. Also, because the particles 

that represent advection in the method of

characteristics are discrete in nature and 

because the concentrations tracked on particles 

are translated to the finite-difference nodes for 

the purpose of computing the mass balance, the 

mass balance error will typically exhibit an 

oscillatory behavior over time. However, this 

is not a cumulative type of error; it is usually 

largest for the first few time increments and 

then tends to balance out over time. As long as 

the oscillations remain within a steady range, 

not exceeding about ±10 percent as a guide, 

then the error probably does not represent a 

bias and is not a serious problem. Rather, the 

oscillations only reflect the fact the mass 

balance calculation is itself just an 

approximation. Thus, it is recommended that 

users examine the solute budget and residuals 

for their particular problem. The significance 

of the residual and rates of change in the 

residual should be assessed qualitatively 

relative to the nature of a particular problem, 

and not merely on the basis of the magnitude of 

the error at any one time increment.

Special Problems

As noted by Konikow and Bredehoeft 

(1978), a number of special problems are 

associated with the use of the method of 

characteristics to solve the solute-transport 

equation. Some of these problems are 

associated with the movement and tracking of 

particles, whereas other problems are related to 

the computational transition between the 

concentrations of particles within a cell and the 

average concentration at that node. This 

section describes the more significant problems 

and the procedures used to minimize errors that 

might result from them.

One possible problem is related to no-flow 

boundaries. Neither water nor solutes can be 

allowed to cross a no-flow boundary. 

However, under certain conditions it might be 

possible for the interpolated velocity at the

26
location of a particle near a no-flow boundary 

to be such that the particle will be advected 

across the boundary during one time increment. 

Figure 7 illustrates such a situation, which 

arises from the deviation between the curvilinear flow line and the linearly projected 

particle path. Figure 7 also shows the 

correction scheme built into the MOC3D 

model. If a particle is advected across a noflow boundary, then it is relocated within the 

aquifer by reflection across the boundary. This 

correction thus will tend to relocate the particle 

closer to the true flow line. However, 

extensive testing indicates that it is unlikely that 

a particle will ever cross a no-flow boundary 

unless CELDIS > 1, which is not recommended.

The maintenance of a reasonably uniform 

and continuous spacing of particles requires 

special treatment in areas where strong fluid 

sources and sinks dominate the flow field. 

Strong fluid sources and sinks cause significant 

convergence and divergence in the flow field, 

which will degrade the desired uniform spacing 

of particles. Without special provisions, 

particles will continually move out of a cell that 

represents a strong fluid source, but few or 

none will move in to replace them and thereby 

maintain a continuous stream of particles. 

Thus, whenever a particle that originated in a 

strong fluid-source cell moves out of that 

source cell, a new particle is introduced into the 

source cell to replace it. Placement of new 

particles in a source cell is compatible with and 

analogous to the generation of fluid and solute 

mass at the source. On the other hand, if a 

fluid source or sink is very weak, it will not 

induce significant divergence or convergence in 

the flow field and have any noticeable effect on 

particle spacing. For cells representing weak 

fluid sources or sinks, particles need not be 

added or removed. The model user must 

specify explicitly whether fluid sources and 

sink are to be flagged as either weak or strong, 

so that particle tracking is implemented 

appropriately. Source/sink cells are identified

EXPLANATION

Zero-transmissivity cell (or 

no-flow boundary)

Flow line and direction of flow 

 Computed path of particle 

© Node of finite-difference cell 

x Previous location of particle 

0 Computed new location of particle 

A Corrected new location of particle

Figure 7. Possible movement of a particle near 

an impermeable (no-flow) boundary (modified 

from Konikow and Bredehoeft, 1978).

as either strong or weak for purposes of 

particle control by the user-specified value of 

the IGENPT array (see Appendix B).

The procedure used to replace particles in 

source cells is illustrated in figure 8. A steady, 

uniformly spaced stream of particles is 

maintained by generating a new particle in the 

source cell at the original location of the particle 

that left the source cell. When a relatively 

strong fluid source is imposed on a relatively 

weak regional flow field, as illustrated in figure 

8a, then radially divergent flow will be 

maintained in the area of the source, and all 

initial and replacement points will move 

symmetrically away from node j,i,k. For 

example, after particle 7 moves from cell (j,i,k) 

at the start of a time increment to cell 

(j+l,i-l,k) at the end of that time increment, the 

replacement particle (particle 18 in fig. 8a) is 

positioned at time n in cell (j,i,k) at the same 

location as the initial position of particle 7.

27
time n-7 time n

13

12 Og

16

10'

°4

12

11

'13

(a)

'14

15

time n-1 time n

(b)

EXPLANATION

Strong fluid source

Particle pathline for 

one time increment

O

A

Node of finite-difference grid 

Initial location of particle p 

New location of particle p 

Location of replacement particle p

Figure 8. Replacement of particles in fluid-source cells (a) for case of negligible regional flow 

and (b) for case of relatively strong regional flow.

Although we normally expect particles to be 

advected out of fluid-source cells, figure 8b 

demonstrates the possibility that particles may 

sometimes enter a source cell. When a 

relatively weak fluid source is imposed on a 

relatively strong regional flow field, the 

velocity distribution within the source cell does

not possess radial symmetry, and the velocity 

within the upgradient part of the source cell is 

lower than the velocity within the downgradient 

part of the source cell. It is possible then that 

particles originating in upgradient cells (such as 

particle 2) will migrate through the designated 

source cell. This can also occur when two or

28
more source cells of different strengths are 

adjacent to each other. Particles that leave a 

source cell, but did not originate in it, are not 

replaced because that would ultimately lead to 

mass balance errors in downgradient areas as 

proportionally too many particles will exist in 

downgradient areas relative to the volume of 

water being added at the source cell. Some 

particles that originated in a strong source cell 

may take more than one time increment to leave 

that cell (for example, particles 6 and 7 in fig. 

8b). Although not illustrated in fig. 8, when 

these lower-velocity particles do leave the 

source cell at time n+l, the replacement 

particles will be placed at the original positions 

of the particles (which was at time n-\) rather 

than at their positions immediately prior to 

leaving (which was at time ri).

Hydraulic sinks also require some special 

treatment. Particles will continually move into 

a cell representing a strong sink, but few or 

none will move out. To avoid the resultant 

accumulation and stagnation of tracer particles, 

any particle moving into a strong-sink cell is 

removed from the flow field after the 

calculations for that time increment have been 

completed. The numerical removal of particles 

that enter sink cells is analogous to the 

withdrawal of fluid and solute mass through 

the hydraulic sink. If the relative magnitude of 

discharge from a sink cell is not strong enough 

to maintain radially convergent flow, and a 

particle exits from the sink cell, then a 

replacement particle will be placed at the center 

of the cell. The combination of creating new 

particles at sources and destroying old particles 

at sinks will tend to maintain the total number 

of points in the flow field at a nearly constant 

value.

Both the flow model and the transport 

model assume that sources and sinks act 

uniformly over the entire area or volume of the 

cell surrounding a source or sink node. Thus, 

in effect, heads and concentrations computed at 

source or sink nodes represent average values

over the area or volume of the cell. Part of the 

total concentration change computed at a source 

node represents mixing between the source 

water at one concentration and the ground 

water at a different concentration (eq. 17). It 

can be shown from the relation between the 

source concentration (C//^) and the aquifer 

concentration at the start of a time increment 

( Cjjtk\ that the following constraints generally 

must be met in a source cell:

and

for

for

:"ji (68a)

If it is assumed that there is complete 

mixing between the source fluid and the 

resident fluid within the volume of a strong 

source cell, then these same constraints should 

also apply to all points within the cell. Because 

of the possible deviation of the concentrations 

of individual particles within a source cell from 

the average concentration, the change in 

concentration computed at a source node 

should not be applied directly to each of the 

particles in the cell. Rather, at the end of each 

time increment for solving the transport 

equation the concentration of each particle in a 

strong source cell is updated by setting it equal 

to the final nodal concentration. Although this 

may introduce a small amount of numerical 

dispersion by eliminating possible concentration variations among particles within the 

volume of the cell, it prevents the adjustment of 

the concentration for any individual particle in 

the source cell to a value that would violate the 

constraints indicated by equation 66.

In areas of divergent flow, a problem may 

arise because some cells can become void of 

particles where pathlines become spaced widely 

apart. This can occur, even in the absence of a 

strong fluid source, because of heterogeneities 

or boundary conditions. This would result in a 

calculation of no change in concentration at a 

node due to advective transport, although the

29
nodal concentration would still be adjusted for 

changes caused by hydrodynamic dispersion. 

Also, some numerical dispersion is generated at 

nodes in and adjacent to the cells in which the 

advective transport of solute was 

underestimated because of the resulting error in 

the concentration gradient. This might not 

cause a serious problem if only a few cells in a 

large grid became void or if the voiding were 

transitory (that is, if upgradient points were 

advected into void cells during later or 

subsequent time increments). Figure 8a 

illustrates radial flow, which represents the 

most severe case of divergent flow. Here it can 

be seen that when four points per cell are used 

to simulate advective transport, then in the 

numerical procedure four of the eight 

surrounding cells would erroneously not 

receive any solute by advection from the 

adjacent source. If eight uniformly spaced 

particles per cell were used initially, then at a 

distance of two rows or columns from the 

source only 8 of 16 cells would be on pathlines 

originating in the source cell. So while 

increasing the initial number of points per cell 

would help, it is obvious that for purely radial 

flow, an impractically large initial number of 

points per cell would be required to be certain 

that at least one particle pathline passes from 

the source through every cell in the grid. 

Because the MOC3D model is based on a 

rectangular Cartesian coordinate system, it is 

not recommended for applications to a purely 

radial flow problem. However, if radial flow 

is localized within a predominantly nonradial 

regional flow field, then satisfactory results 

should be achievable.

The problem of cells becoming void of 

particles can be minimized by limiting the 

number of void cells to a small fraction of the 

total number of active cells that represent the 

aquifer. The user specifies this fraction 

(FZERO) in the MOC3D input data file (see 

Appendix B). If the limit is exceeded, the 

numerical solution to the solute-transport 

equation is halted temporarily at the end of that

time increment and the "final" concentrations at 

that time are saved. Next the problem is 

reinitialized at the time of termination by 

regenerating the initial particle distribution 

throughout the grid and assigning the "final" 

concentration at the time of termination as new 

"initial" concentrations for nodes and particles. 

The solution to the solute-transport equation is 

then simply continued in time from this new set 

of "initial" conditions until the total simulation 

period has elapsed. This procedure preserves 

the mass balance within each cell but also 

introduces a small amount of numerical 

dispersion by eliminating variations in 

concentration within individual cells.

Review of MOC3D Assumptions and 

Integration with MODFLOW

Following is a brief summary of model 

application assumptions that have been 

incorporated into the MOC3D model. These 

are relevant to both grid design and model 

implementation. Efficient and accurate use of 

MOC3D requires the user to be aware of all of 

these assumptions and options.

 MOC3D is integrated with MODFLOW96 (Harbaugh and McDonald, 1996a) and will 

not work with earlier versions. The main 

MODFLOW subroutine is replaced with the 

MOC3D main subroutine. In addition, several 

MOCJD-specific source code files must be 

compiled and linked to the MODFLOW code.

 Particle velocities are interpolated 

spatially, but not over time. That is, we 

assume that the head distribution calculated for 

the end of a given time step applies during that 

entire time interval.

 The transport model is applied to a 

"window" of the grid used to solve the flow 

equation. This subgrid can be equal in size or 

smaller than the primary MODFLOW grid.

 Within the area of the transport subgrid, 

row and column discretization must be

30
uniformly spaced (that is, Ax and Ay must be 

constant, although they need not be equal to 

each other). The spatial discretization or rows 

and columns beyond the boundaries of the 

subgrid can be nonuniform, as allowed by 

MODFLOW, to permit calculations of head 

over a much larger area than the area of interest 

for transport simulation (see fig. 9). Vertical 

discretization, defined by the cell thickness, can 

be variable in all three dimensions. However, 

large variability may adversely affect numerical 

accuracy (as discussed in second item below).

COLUMNS

EXPLANATION

Area of transport subgrid 

within primary MODFLOW grid

Figure 9. One layer of finite-difference grid 

illustrating use of uniformly-spaced transport 

subgrid for MOC3D within variably-spaced 

primary grid for MODFLOW.

 Retardation factor values and all 

dispersivity values are constant in each layer. 

Values for porosity may vary within a layer and 

are defined for each node (also see discussion 

in next item).

 The particle-tracking algorithm inherently 

assumes that all particles represent the 

concentration in an equal volume of water in a 

cell, where the volume equals (ebAxAy)^. 

Thus, although Ax and Ay are uniform, it is 

also very important that the variations in the 

product of porosity and thickness within the 

transport subgrid remain relatively small.

Otherwise, when a particle moves into a cell 

having a very different volume from the cell in 

which it originated, the estimate of the average 

concentration in the new cell may become 

biased, which will also have an adverse effect 

on the overall mass balance for the solute.

 MODFLOW offers flexibility to the user 

in the conceptualization of vertical discretization 

(see McDonald and Harbaugh, 1988, Ch. 2). 

As illustrated in fig. lOa, it is common in 

applications of MODFLOW to represent the 

resistance to flow in a low hydraulic 

conductivity unit by lumping the vertical 

hydraulic conductivity and thickness of the 

confining unit into the vertical conductance 

term between the adjacent layers (fig. lOc). 

However, because transport simulation 

requires that travel distances be known 

explicitly in all directions, three-dimensional 

transport simulation requires fully threedimensional flow simulation (rather than a 

quasi-three-dimensional analysis) within the 

area of the transport subgrid. That is, even if 

the solution to the flow equation is insensitive 

to heads and storage releases in the clay layer, 

it must still be represented by one or more 

model layers for the solution to the transport 

equation (fig. lOb). If not, it can be seen in the 

quasi-three-dimensional analysis (fig. lOc) that 

any solute crossing the lower boundary of layer 

1 would immediately be located in and 

influenced by the properties of layer 2, and 

would never have been subject to the relatively 

long travel time through the clay.

 In MODFLOW, layers may be defined as 

representing confined or unconfined 

conditions, or allowed to switch between these 

two types, depending on the specification of 

the value of the parameter LAYCON. If 

LAYCON>0, then the user also must specify 

information about the elevations of the top and 

bottom of the layer. For layers having 

LAYCON=0 or LAYCON=2, MODFLOW 

assumes that transmissivity and saturated 

thickness remains constant; hence, the

31
SYSTEM

CLAY

THREE-DIMENSIONAL 

MODEL

QUASI-THREE-DIMENSIONAL 

MODEL

Layer 1 V 

Layer 2

Layer 1 ^X'vVv'/v'/v

The clay layer is represented by vertical 

conductance terms between layers 1 and 2

Layer 3 / Layer 2 :y/£$"ify.

(b) (c)

Figure 10. Alternative MODFLOW approaches to vertical discretization of an aquifer system (a) 

consisting of two high-permeability units separated by a confining layer consisting of a lowpermeability clay. In the fully three-dimensional representation (b), the clay unit would be 

represented by one or more model layers. In the quasi-three-dimensional approach (c), heads are 

not calculated in the clay unit, which is represented more simply by vertical conductance terms 

between the layers above and below it; in this case, the bottom of layer 1 coincides with the top of 

layer 2. (Modified from figs. 10-12 of McDonald and Harbaugh, 1988.)

thickness values read from the MOC3D input 

files are used. However, for layers having 

LAYCON=1 or LAYCON=3, MODFLOW 

allows the transmissivity to change as a 

function of changes in saturated thickness in 

each cell on the basis of elevation data that are 

input for the MODFLOW BCF package. In 

these cases, thickness for a cell is defined as 

"TOP" minus "EOT" if the layer is confined 

and "HEAD" minus "EOT" if it is unconfmed; 

the thickness values specified in the MOC3D 

input file are ignored in the calculations. 

However, note that thickness values must still 

be specified in the MOC3D input file in all 

cases.

 Concentrations associated with fluid 

sources are read directly from MODFLOW

source/sink package files. In each of the 

package files used, the options 

"CBCALLOCATE" and "AUXILIARY 

CONC" must be included on the first line of 

input data. See MODFLOW and MOC3D input 

instructions.

 If the evapotranspiration package (EVT) is 

implemented, MODFLOW will calculate a fluid 

discharge (or sink) rate that is typically 

associated with an evapotranspirative process 

that removes water but excludes dissolved 

solids, which are retained in the remaining fluid 

at a consequently higher concentration. 

MOC3D assumes that for any such calculated 

flux, the associated source concentration ( C' in 

eq. 17) will equal 0 rather than equaling the 

concentration at the node, as is assumed

32
normally for a fluid sink. This will induce an 

appropriate increase in concentration at a cell 

representing a fluid sink due to 

evapotranspiration. Note that this MOC3D 

assumption should be viewed as a first-order 

approximation because in actuality (1) the 

evapotranspirative process may not be 100 

percent effective in excluding solutes, 

depending on the particular chemical species, 

and (2) calculated solute concentrations may 

exceed the upper limits of solubility for a 

particular chemical constituent (and MOC3D 

does not simulate mineral precipitation).

 When the solute of interest is subject to 

decay, it is assumed that the solute in liquid and 

solid phases will decay at the same rate. If a 

fluid source contains the decaying solute, it is 

subject to decay after it enters the ground-water 

system, but is not decayed within its "source 

reservoir."

 All unit numbers specified in the name 

files for a particular simulation must be unique. 

Unit numbers 99, 98, and 97 are reserved for 

the MODFLOW name file, the batch mode 

input file, and the batch mode output file, so 

cannot be specified for any other use. 

However, unit numbers may be reused in 

separate simulation runs in batch mode.

 The model includes output options to 

create separate binary data files (Ftypes CNCB, 

VELB, and PRTB); when implemented, the 

model will write calculated values from the 

simulation for the selected variables as 

unformatted data. The concentration and 

velocity files (CNCB and VELB) use the 

MODFLOW module ULASAV to write the data 

(see MODFLOW documentation). When the 

velocity option VELB is implemented, the code 

will first write the velocities in the column 

direction at all nodes, then all velocities in the 

row direction, and finally all of the velocities in 

the layer direction. The velocity and concentration arrays are dimensioned to the size of the 

transport subgrid only. When particle data are 

written to a separate binary file (PRTB), the file

begins with a header line that includes the move 

number, number of particles, and length of the 

transport time increment. A record for each 

particle in sequence follows the header line and 

contains the location and concentration of each 

particle in the following order: column 

coordinate, row coordinate, layer coordinate, 

and particle concentration.

COMPUTER PROGRAM

MOC3D is implemented as a package for 

MODFLOW. MOC3D uses the flow components calculated by MODFLOW to compute 

velocities across each cell face in the transport 

domain. The computed velocities are used in 

an interpolation scheme to move each particle 

an appropriate distance and direction to 

represent advection. The effects of fluid 

sources, dispersion, and decay on concentration are then applied to the particles.

A separate executable version of 

MODFLOW, which is adapted to link with and 

use the MOC3D module, must first be 

generated and then used to run MOC3D 

simulations. The MOC3D code is written in 

standard FORTRAN, and it has been 

successfully compiled and executed on multiple 

platforms, including 486- and Pentium-based 

personal computers, Macintosh personal 

computers, and Data General and Silicon 

Graphics Unix workstations. FORTRAN 

compilers for each of these platforms vary in 

their characteristics and may require the use of 

certain options to successfully compile 

MOC3D. For instance, the compiler should 

initialize all variables to zero. Depending on 

the size of the X-array (defined by LENX in 

the MODFLOW source code), options to 

enable the compiler to handle a large array may 

be needed.

Implementing MOC3D requires the use of a 

separate "name" file similar to the one used in 

MODFLOW. The principal MOC3D input data

33
(such as subgrid dimensions, hydraulic 

properties, and particle information) are read 

from the main MOC3D data file. Other files are 

used for observation wells, concentrations in 

recharge, and several input and output options. 

Detailed input data requirements and 

instructions are presented in Appendix B. 

Also, a sample input data set for a test problem 

is included in Appendix C.

MOC3D output is routed to a main listing 

file, separate from the MODFLOW listing file. 

There are also several options for writing 

specific data to separate output files, which will 

facilitate graphical postprocessing. Appendix 

D contains output from the sample data set 

described in Appendix C.

General Program Features

Because the model assumes that changes in 

concentration do not affect the fluid properties 

(such as density and viscosity), the head

distribution and flow field are independent of 

the solution to the solute-transport equation. 

Therefore, the flow and transport equations can 

be solved sequentially, rather than simultaneously. But because transport depends on fluid 

velocity, which is calculated from the solution 

to the flow equation, the sequence order must 

be to solve the flow equation first. This 

sequence is illustrated in figure 11 for a 

hypothetical problem involving transient flow 

and three stress periods. The numbered 

sequence from 1 through 16, which starts at the 

left edge of the double time line, illustrates the 

order of solving equations as the simulation 

progresses through the first five time steps of 

the first stress period in this hypothetical 

example. This figure also helps to illustrate the 

nomenclature used for time parameters in 

MODFLOW and MOC3D, as well as the 

relation between them.

The implicit solution to the flow equation in 

MODFLOW generally allows the use of time 

steps of increasing length during a given stress

At

(Elapsed time, 

flow) 

TOT1M 

(Elapsed time, 

transport) 

SUMTCH -

-NMOV=3 

 IMOV=1

-TIME INCREMENT 

FOR TRANSPORT (TIMV)

tfotal

Figure 11. Double time-line illustrating the sequence of progression in the MOC3D model for solving 

the flow and transport equations. This example is for transient flow and three stress periods (NPER = 

3) of durations PERLENf, PERLEN2, and PERLENs. Each time step for solving the flow equation (of 

duration DEL7) is divided into one or more time increments (of duration TIMV) for solving the transport 

equation; all particles are moved once during each transport time increment. For illustration purposes, 

the sequence of solving the two equations is labeled for the first five time steps of the first stress 

period, and the indices for counting time steps for flow and time increments for transport are labeled 

for the fourth time step.

34
period. The length of the first time step for 

solving the flow equation is calculated by 

MODFLOW on the basis of user-defined 

values for the number of time steps (NSTP), a 

time-step multiplier (TSMULT), and the length 

of the stress period (PERLEN). After the flow 

equation is solved for the first time step (At]), 

the model compares the length of the time step 

for the flow equation with the limitations 

imposed by the stability and accuracy criteria 

for solving the transport equation. If any 

criteria are exceeded, MOC3D will subdivide 

the time step into the fewest number of equalsized time increments that meet all of the 

criteria. In the example shown in figure 11, the 

first two time steps are small enough so that the 

transport equation can be solved for a single 

time increment of the same duration as the flow 

time step (that is, TIMV = DELT). As this 

equation-solving sequence progresses and is 

repeated for increasingly long time steps, the 

stability criteria are eventually exceeded. 

Figure 11 shows that for the third and 

subsequent time steps, the transport equation 

had to be solved over shorter time increments. 

Note that because time increments for transport 

are the same length (TIMV) during any given 

time step for flow, the length of the transport 

time increments will generally be slightly 

different between any two different flow time 

steps. For example, the length of the three 

transport time increments during the fourth 

flow time step (m = 4) are slightly different 

than the lengths of the four time increments 

during flow time step 5. At any point during 

the progress of the simulation, the elapsed time 

for transport is always less than or equal to the 

elapsed time for flow.

Transport may be simulated within a 

subgrid, which is a "window" within the 

primary MODFLOW grid used to simulate flow 

(see fig. 9). The grid dimensions are limited 

only by the size of the "X" array (see "Space 

Allocation" in the MODFLOW documentation). 

Within the subgrid, the row and column

spacing must be uniform, but thickness can 

vary within a small range from cell to cell and 

layer to layer.

Many MOC3D subroutines are linked 

closely with MODFLOW counterparts. When 

possible, MOC3D follows MODFLOW 

subroutine structure. In general, data are 

defined, space is allocated in the "X" array, and 

simulation parameters are read just as in 

MODFLOW. The overall structure of the 

MOC3D code and its integration with 

MODFLOW art illustrated in figure 12, which 

shows a flow chart for the main program 

(excluding details of the transport calculations).

A more detailed flow chart of the program 

segments controlling transport calculations is 

shown in figure 13. The fluxes that 

MODFLOW calculates within the transport 

subgrid are processed by MOC3D subroutines 

to generate a transient solution to the solutetransport equation.

Program Segments

MOC3D input and output utilizes the 

standard MODFLOW array reading and writing 

utilities as much as possible. MOC3D also 

takes advantage of new features in 

MODFLOW, such as the option for auxiliary 

parameters in the source and sink packages and 

storing budget flows from each of those 

packages, as documented by Harbaugh and 

McDonald (1996a and 1996b). However, 

many subroutines in MOC3D do not fit into a 

MODFLOW module class. For those model 

users who are interested in more details about 

the internal structure and organization of the 

code, Tables 1-10 list and describe briefly each 

of the subroutines in MOC3D that are used for 

ten different categories of functions.

Tables 1 and 2 include subroutines that 

initialize and set up the transport simulation. 

MOC3D data are read and checked for 

consistency with each other and with several 

MODFLOW parameters. Seepage velocities

35
are calculated on the basis of hydraulic 

gradients determined by the MODFLOW 

solution to the flow equation. The VELO 

subroutine (Table 3) is called up to three times 

after each solution to the flow equation is

obtained once for each dimension of the 

simulation.

Table 4 includes subroutines controlling 

particle tracking, concentration calculations, 

and related output. The SMOC5GP subroutine

Define 

problem

 » Allocate 

space

Read and 

prepare data

Q. 

O 

_O

TJ 

.0

CD 

Q.

CO

Calculate time step 

length and set heads

Set initial saturated thickness 

for water-table cells

Calculate finitedifference equation 

coefficients

Calculate initial 

solute mass

Generate initial 

particles

Q. 

o 

.g

o

Calculate flow and transport 

terms from storage

Calculate flow and transport 

terms from fixed heads

Calculate flow across 

boundary of subgrid

Compute flow between 

adjacent cells

Compute and 

print velocities

Formulate and print dispersion

equation coefficients; check

stability criteria

Calculate flow and transport 

terms from source/sink cells

Yes

Yes

Transport loop 

(See Fig. 13)

Figure 12. Generalized flow chart for MOC3D.

36
Compute transport 

time increment

8.

(0 

05

Move particles and

account for decay

(see Fig. 17)

I

Compute intermediate 

node concentrations

Approximate change in concentration 

due to sources/sinks

Approximate change in concentration 

due to dispersion

Compute new node and 

particle concentrations

Compute solute 

mass balance

Print concentrations and 

mass balance (optional)

ore moves to 

complete time 

step?

Figure 13. Simplified flow chart for the transport 

loop, which is shown as a single element in fig. 

12.

generates the initial coordinates for each tracer 

particle in the transport subgrid. Default 

particle placement patterns include several 

configurations for simulations in one, two, or 

three dimensions. Each of the configurations 

distributes the particles uniformly in space with 

respect to dimension. Figures 14 through 16 

illustrate the distribution of particles with each 

of the default number of particles per node. 

Note that in some cases, it may be necessary to 

customize the initial positions of the particles

(see the Two Dimensional Radial Flow and 

Dispersion test case in the Model Testing and 

Evaluation chapter). For a given transport time 

increment, particles are moved a distance and 

direction on the basis of the estimated velocity 

at the location of each particle and the length of 

the time increment. The particle velocity is 

estimated by interpolation from the velocities 

on adjacent cell faces to the location of a 

particular particle. Either linear or bilinear 

interpolation is used (in subroutines MOVE or 

MOVEBI, respectively) based on the userselected value for the INTRPL flag (see 

MOC3D Input Instructions). However, when 

the bilinear interpolation option is used, particle 

velocity in the vertical (or "layer") direction will 

still be interpolated linearly. A flow chart 

describing details of the MOVE and MOVEBI 

subroutines is presented in figure 17.

Table 1. MOC3D subroutines controlling simulation preparation

Subroutine Description

SMOC5O Open MOC data files

MOC5DF Define subgrid and other key 

parameters

MOC5AL Allocate space in "X" array to 

store MOC data

MOC5RP Read MOC parameters

MOC5CK Check MOC data for 

consistency

SMOC5Z Set an array to a specified 

constant value

Table 2. MOC3D subroutines controlling transport 

time factors

Subroutine Description

MOC5ST Check limiting stability 

criteria; compute time 

increment and number of 

moves for solute transport

MOC5AD Update elapsed transport time

37
Table 3. MOC3D subroutines controlling velocity 

calculations and output

Subroutine Description

VELO Calculate velocities from 

flows across cell faces

SMOC5V Output velocity data

Table 4. MOC3D subroutines controlling particle 

tracking, concentration calculations, and output for 

particle and concentration data

Subroutine Description

SMOC5GP Generate initial particle 

distribution

SMOC5P Output particle locations and 

concentrations

MOVE and Advect particles; compute 

MOVEBI concentrations at the end of 

each move; decay particle 

concentrations (MOVE uses 

linear interpolation of 

velocity and MOVEBI uses 

bilinear interpolation)

MOVTM Compute time for particle to 

reach boundary of cell

MOC5AP Compute new node and

particle concentrations at 

end of move

SMOC5C Output node concentrations

(a) (b)

(c) (d)

Figure 14. Default initial particle configurations 

for a one-dimensional simulation using (a) one, 

(b) two, (c) three, and (d) four particles per cell.

(a) (b)

(c) (d)

Figure 15. Default initial particle configurations 

for a two-dimensional simulation using (a) one, 

(b) four, (c) nine, and (d) sixteen particles per 

cell.

(a)

r

- __

(b) (c)

Figure 16. Default initial particle configurations for a three-dimensional simulation using (a) one, 

(b) eight, and (c) 27 particles per cell. For clarity, only the volume of a single cell is illustrated.

38
Subroutines related to dispersion 

calculations are included in Table 5. 

Dispersion coefficients are determined on cell 

faces. However, to facilitate improved 

efficiency in the code the dispersion

coefficients are lumped with the porosity, 

thickness, and an appropriate grid dimension 

factor of the cell into combined parameters 

called "dispersion equation coefficients." For 

example, the dispersion equation coefficient for

Move particle to 

cell boundary

Will particle 

reach boundary?

Compute minimum time 

to cell boundaries Select next particle 

and get coordinates

Calculate velocity of

particle (linear or

bilinear interpolation)

Is particle in 

an active cell? Move particle to 

new location

Sum decayed concentrations of

all particles in cell; increment

number of particles in current cell

Is particle 

outside of 

subgrid?

No ^ Did particle cross 

cell boundary?

Remove particle from

active array; place

in inactive buffer Is old cell

a strong fluid

source or sink?

Did particle 

leave an inflow 

boundary cell?

Is new

location a discharge 

undary or strong 

sink?

Place new particle 

at same relative position 

as old particle in new cell

Place new particle in 

old cell at node location

Sum decayed concentrations of

all particles in cell; increment 

number of particles in current cell

Did particle

originate in

that cell?

Place new particle 

in old cell at same 

originating position

Figure 17. Flow chart for MOVE and MOVEBI subroutines of MOC3D, which is shown as a single 

element in fig. 13.

39
the j+ 1/2, i, k face in the column direction is 

(ebD\

V **//+l/2,i.fc

Ax

These combined coefficients are the ones that 

are written to the output files by SMOC5D. 

Appendix A includes a more detailed 

description of the dispersion equation 

coefficients.

Table 5. MOC3D subroutines controlling dispersion calculations and output

Subroutine Description

DSP5FM Calculate dispersion 

coefficients

SMOC5D Output dispersion equation 

coefficients

DSP5AP Use explicit finite-difference 

formulation to compute 

changes in concentration 

due to dispersion

Subroutines that link the MOD FLOW 

source/sink package calculations of fluid flux to 

the MOC3D calculations of solute concentration 

and solute flux are listed in Table 6. 

MODFLOW source and sink packages contain 

an option called CBCALLOCATE. When 

used, the package will save the cell-by-cell 

flow terms across all faces of every source or 

sink cell. MOC3D uses these fluid fluxes to 

calculate solute flux to or from the source/sink 

nodes. Because these individual solute fluxes 

are required to compute the solute mass 

balance, the CBCALLOCATE option must 

always be selected when using MOC3D. 

Calculations of concentration changes at nodes 

caused by mixing with fluid sources are 

controlled by the "SRC" subroutines listed in 

Table 7.

Subroutines controlling observation well 

features are listed in Table 8. Table 9 lists 

subroutines related to the solute mass balance

calculations. Table 10 lists subroutines related 

to calculating fluid storage terms in the solutetransport equation.

Table 6. MOC3D subroutines controlling 

MODFLOW source/sink package calculations

Subroutine Description

CDRN5FM Calculate solute flux to drains

CEVT5FM Calculate solute flux to 

evapotranspiration

CGHB5FM Calculate solute flux to/from 

general head boundary cells

CRCH5AL Allocate space in "X" array 

for concentrations 

associated with recharge

CRCH5FM Calculate solute flux from 

recharge

CRTV5FM Calculate solute flux to/from 

river cells

CWEL5FM Calculate solute flux to/from 

well cells

Table 7. MOC3D subroutines controlling cumulative calculations relating to sources and sinks

Subroutine Description

SRC5FM Initialize source/sink array; 

compute terms at fixed 

heads

SRC5AP Calculate changes in

concentration due to flux at 

sources and sinks

Table 8. MOC3D subroutines controlling observation wells

Subroutine Description

OBS5DF Read number of observation 

wells

OBS5AL Allocate space in "X" array 

for observation well data

OBS5RP Read observation well 

locations

SOBS5O Output observation well data

40
Table 9. MOC3D subroutines controlling mass 

balance calculations and output

Subroutine Description

SMOC5IM Calculate initial solute mass 

stored

SMOC5BY Calculate flows across

boundaries of subgrid (for 

mass balance)

SMOC5BD Compute cumulative solute 

mass balance

SMOC5M Output mass balance 

information

Table 10. Miscellaneous MOC3D subroutines

Subroutine Description

SMOC5TK Set initial saturated thickness 

for water-table cells

SMOC5UP Update fluid storage terms for 

transport equation

MODEL TESTING AND 

EVALUATION

In developing and documenting a new 

numerical model, it must be demonstrated that 

the generic model can accurately solve the 

governing equations for various boundary 

value problems. This is accomplished by 

demonstrating that the numerical code gives 

good results for problems having known 

solutions, such as those for which an analytical 

solution is available.

The accuracy of numerical solutions is 

sometimes sensitive to spatial and temporal 

discretization. Therefore, even a perfect 

agreement for selected test cases proves only 

that the numerical code can accurately solve the 

governing equations, not that it will under any 

and all circumstances.

Analytical solutions generally require that 

an aquifer can be assumed to have simple

geometry, uniform properties, and idealized 

boundary and initial conditions. A major 

advantage of numerical methods is that they 

relax the simplifications required by analytical 

methods and allow the representation of more 

realistic field conditions, such as heterogeneous 

and anisotropic properties, irregular geometry, 

mixed boundary conditions, and multiple 

stresses that vary in time and space. However, 

analytical solutions approximating these 

complexities are unavailable for comparison. 

Therefore, it is difficult to prove that the 

numerical models can accurately solve the 

governing equations for the very situations for 

which they are most needed. For such cases, 

we are limited to relatively simple tests, such as 

benchmarking and evaluating the global massbalance error. In the benchmarking approach, 

we compare the results of the MOC3D model 

for selected complex problems to results of 

other well accepted models. Although 

benchmarking is useful to improve confidence 

in the model, it is largely a measure of 

consistency and does not guarantee or measure 

accuracy. Overall, we have attempted to test 

and evaluate the MOC3D model for a range of 

conditions and problem types so that the user 

will gain an appreciation for both the strengths 

and weaknesses of this particular code. Additional testing and benchmarking of MOC3D 

is documented in Goode and Konikow (1991).

One-Dimensional Steady Flow

Wexler (1992) presents an analytical 

solution for one-dimensional solute transport in 

a finite-length aquifer system having a thirdtype source boundary condition. The 

governing equation is subject to the following 

boundary conditions:

VC' = VC-D dC_ 

dx

(69)

41
dx

x = L (70)

and the following initial condition:

C = 0, 0<x<L. (71)

For this test problem we assumed that the 

length of the system, L, is equal to 12 cm, 

(7 = 1.0, and V = 0.10 cm/s. The analytical 

solution is given by equations 52 and 53 of 

Wexler (1992, p. 17). In generating an 

equivalent solution using MOC3D, we set up a 

one-dimensional grid having 122 cells (nodes) 

in the ^-direction within which the flow 

equation was solved. The solute-transport 

equation was solved in a 120-cell subgrid to 

assure a constant velocity within the transport 

domain and to allow an accurate match to the 

boundary conditions of the analytical solution. 

The grid spacing was Ax = 0.1 cm. The 

numerical solution was implemented using 

three initial particles per cell (NPTPND = 3) 

and a CELDIS factor of 0.5. The input 

parameters for the model simulation are 

summarized in Table 11.

Two different values of dispersion 

coefficients were evaluated in the first set of 

tests. The values were Dxx = 0.1 and 0.01 

cm2/s, which are equivalent to aL = 1.0 and 

0.1 cm, respectively. Breakthrough curves 

showing concentration changes over time at 

three different locations as calculated with both 

the analytical and numerical solutions for the 

lower dispersion case are compared in figure 

18. To improve clarity, this plot only shows 

every fourth data point for the numerical model 

results, except for the curve for x = 0.05 cm, 

where every data point is shown for times less 

than 10 seconds. Note that this distance (x = 

0.05) is the first node downgradient from the 

source location. With the possible exception of 

very early time at locations very close to the 

source, there is essentially an exact fit between 

the numerical and analytical solutions. At early 

times and short distances the numerical solution 

exhibits some nonsmoothness and oscillation

Table 11. Parameters used in MOC3D simulation 

of transport in a one-dimensional, steady-state 

flow system

Parameter Value

= T

yy

&TH ~ aTv

PERLEN (length of stress 

period)

Vx

Vy = Vz

Initial concentration (Co)

Source concentration ( C" )

Number of rows

Number of columns

Number of layers

DELR(Ax)

DELC(Ay)

Thickness (b)

NPTPND (Initial number of 

particles per cell)

CELDIS

INTRPL (Interpolation 

scheme)

0.01 cm2/s 

0.1

0.1 cm 

0. 1 cm 

1 20 s

0.1 cm/s

0.0 cm/s

0.0

1.0

1 22

1

1

O.lcm

O.lcm

1 .0 cm

3

0.5

2

about the mean, which is related to the discrete 

nature of the particles used to represent the 

advection process. However, this small loss of 

precision is not a cumulative error, as it 

vanishes after moderate travel times or 

distances.

The results for the higher dispersion case 

are presented in figure 19. Because dispersion 

is a limiting stability criterion and the 

dispersion coefficient is ten times higher in the 

same grid, the transport simulation takes many 

more time increments (or particle moves). 

Thus, in figure 19 only every 100th point is 

shown (as small circles) for the numerical 

solution at the two larger distances. For these 

two curves, the match between the analytical

42
1.0 - OOP O Oi

0

Analytical 

MOC3D

x=11.05

20 40 60 80 

TIME (SECONDS)

100 120

Figure 18. Numerical (MOC3D) and analytical solutions at three different locations for 

solute transport in a one-dimensional, steady flow field. Parameter values for this base 

case are listed in Table 11.

solution and the numerical solution is almost 

perfect. For the location that is close to the 

source, every point is plotted as a small dot to 

illustrate again the small loss of precision for 

short travel distances and times.

To help assess the significance of the 

oscillations and loss of precision at nodes very 

close to the source, the early (times less than 10 

seconds) part of that breakthrough curve is 

replotted at a larger scale in figure 20. The 

oscillation is caused by the fact that the stability 

requirements related to the explicit solution for 

the dispersive flux causes the time increment 

for solving the transport equation to be so small 

that particles used to track the advective flux 

can only move a small fraction of the width of a 

cell during a given time increment. Because the 

distance that the particles move during one time

increment (about 0.005 cm in this case) is 

smaller than the spacing between particles (onethird of Ax for this case in which three particles 

per cell are used, or 0.033 cm), particles only 

cross cell boundaries after every seven moves 

in this case. Therefore, the change in 

concentration caused by advection is 

underestimated during six moves when no 

particles cross a cell boundary and 

overestimated during the seventh move when 

one particle does cross the cell boundary. 

However, there is essentially no cumulative 

error and the numerical solution oscillates 

regularly in a small and decreasing range about 

the true solution. Also, the magnitude of the 

oscillations diminishes over time as dispersion 

reduces the local concentration gradients. To 

check this explanation and to demonstrate that

43
20 40 60 80 

TIME (SECONDS)

100 120

Figure 19. Numerical (MOC3D) and analytical 

solutions at three different locations for solute 

transport in a one-dimensional, steady flow field 

for case of increased dispersivity (ai = 1.0 cm, 

DXX= 0.1 cm2/s, and other parameters as defined 

in Table 11).

the error shown in figure 20 is primarily an 

artifact of having too few initial particles per 

cell for this particular combination of 

parameters, and that it is not a generic 

deficiency in the algorithm, MOC3D was run 

for the same problem using an initial particle 

density of 50 particles per cell. The results for 

the same location for the first 10 seconds are 

shown in figure 21 for comparison. When 50 

particles are used, the distance that each particle 

moves during one time increment (again about 

0.005 cm) is greater than the spacing between 

adjacent particles (0.002 cm). For this case, 

the agreement between the analytical solution 

and the MOC3D results are much closer than in 

figure 20 and the oscillations are almost entirely 

eliminated.

The results of these tests can also be 

presented in the form of breakthrough curves 

that plot concentration against distance for 

various times. Figure 22 shows the results for 

the same set of parameters as shown in figure 

18 (that is, the low dispersion case). For

clarity in figure 22, only every fourth data 

point is plotted for the numerical results, 

except every data point is shown for distances 

less than 1.5 cm on the curve for t = 6 

seconds. The results show an almost perfect 

agreement between the analytical and numerical 

solutions.

1.0

0.8

ir o.e iz

LU 

^0.4

O

o

0.2

0.0

 Analytical

 MOC3D

4 6 

TIME (SECONDS)

10

Figure 20. Detailed view of numerical and 

analytical solutions for early times (f < 10 s) at the 

first node downgradient from the inflow source 

boundary for same problem as shown in fig. 19.

1.0

0.8

<0.6

Iz

LU 

O0.4

Z

o 

o

0.2

0.0

 Analytical

 -MOC3D

x=0.05

4 6 

TIME (SECONDS)

10

Figure 21. Detailed view of numerical and 

analytical solutions for early times (t < 10 s) at the 

first node downgradient from the inflow source 

boundary for same conditions as shown in fig. 

20, except that the initial number of particles per 

node, NPTPND, equals 50.

44
The effect of incorporating a retardation 

factor to represent a linear, equilibrium, 

reversible, sorption process is illustrated in 

figure 23. This shows a comparison between 

the analytical solution and the MOC3D results 

for the same low-dispersion problem 

represented in figures 18 and 22, except that 

the elapsed time is 240 seconds and the three 

different curves are for cases in which Rf= 2, 

4, and 40. Only every fourth data point is 

plotted in figure 23 for the numerical results, 

except for the case for Rf = 40, where every 

data point is shown for distances less than 1.5 

cm. The agreement is excellent. Note that

1.0

z 0.8 

o

<0.6

H0.4

O 

0 0 .2

0.0

t = 60 sec

t = 120 sec

0

Analytical 

MOC3D

468 

DISTANCE (cm)

10 12

Figure 22. Numerical (MOC3D) and analytical 

solutions at three different times for same onedimensional, steady flow, solute-transport 

problem shown in fig. 18.

468 

DISTANCE (cm)

10 12

Figure 23. Numerical (MOC3D) and analytical 

solutions after 240 seconds for three different 

retardation factors for same problem 

represented in fig. 22.

because the net effect of the retardation factor is 

to transform the time scale, the three sets of 

curves and data points shown in figure 23 for t 

= 240 s are identical to the three sets of curves 

and data points shown in figure 22 for shorter 

times andRf = 1.0.

The effect of decay is evaluated by 

specifying the decay rate as A = 0.01 S' 1 for the 

same low-dispersion, no sorption, problem as 

defined for figures 18 and 22. These results 

are presented in figure 24, which shows 

excellent agreement between the analytical and 

numerical solutions. Only every fourth data 

point is plotted in figure 24 for the numerical 

results.

1.0r

468

DISTANCE (cm)

10 12

Figure 24. Numerical (MOC3D) and analytical 

solutions for four different times for solute 

transport in a one-dimensional, steady flow field 

for case with decay at rate of A = 0.01 s~ 1 . All 

other parameters as defined in Table 11.

Three-Dimensional Steady Flow

To further evaluate and test MOC3D for 

three dimensional cases, we compare the 

numerical results with those of the analytical 

solution developed by Wexler (1992) for the 

case of three-dimensional solute transport from 

a continuous point source in a steady, uniform 

flow field in a homogeneous aquifer of infinite 

extent. Relative to the previous tests, and 

because the flow field is aligned with the grid

45
in this test case also, this evaluation primarily is 

a test of the accuracy of the calculated 

dispersive flux in three directions. The 

analytical solution of Wexler (1992, p. 42-43) 

is subject to the following boundary conditions:

= 0, 

= 0,

# = ±00

y = +00

(72a)

(72b)

C = 0, z = ±°° (72c) 

and the following initial condition (at t = 0):

C = 0,

 oo < x < °°, oo < y < oo? oo < 2 < oo. (73)

For the analytical solution to this test problem, 

we assumed that Vx = 0.1 m/d, Vy = Vz = 0.0 

m/d, Dx = 0.06 m2/d, Dy = 0.003 m2/d, Dz = 

0.0006 m2/d, e = 0.25, and that the source 

strength (or solute mass flux) is QCf - 10.0 

g/d. The analytical solution is given by 

equation 105 of Wexler (1992, p. 47), and 

assumes that the fluid source does not affect the 

flow field.

Whereas the analytical solution assumes an 

infinite aquifer, the numerical solution can only 

be applied to a finite system. In generating an 

equivalent solution using MOC3D, we aimed to 

use a grid that was sufficiently large so as to 

minimize any effects of the boundaries on the 

solution. Because of the symmetry of the 

problem, we only simulated one quadrant of 

the cross-sectional area of the aquifer 

downgradient from the point source. The 

three-dimensional transport subgrid had 30 

rows of cells (nodes) at a grid spacing of 3 m 

that are parallel to the ^-direction of the 

analytical solution, 12 columns at a grid 

spacing of 0.5 m in the j-direction, and 40 

layers at a grid spacing of 0.05 m in the zdirection, within which the transport equation 

was solved. Boundary conditions, values of 

heads on boundaries, hydraulic conductivity, 

and porosity were specified to assure that the 

velocity would equal 0.1 m/d in the ^-direction. 

Identical values for the dispersion coefficients

were generated by specifying on = 0.6 m, OTH 

= 0.03 m, and aTv = 0.006 m. The point 

source was represented in the numerical model 

by a combination of g = 1.0xlO~6 m3/d and 

C" = 2.5xl06 g/m3 (note that g/m3 i s 

equivalent to mg/L), which together yield onefourth of the source flux assumed in the 

analytical solution. The very small fluid 

injection rate at the point source assures that the 

fluid has only a negligible effect on the flow 

field, as required for consistency with the 

analytical solution. However, another small 

but unavoidable difference between the 

analytical and numerical solutions is that the 

former has the solute source at a true point but 

the numerical model inherently assumes the 

solute source is within the volume of one cell. 

Both models were run for a total elapsed time 

of 400 days. The numerical solution was 

implemented using three initial particles per cell 

(NPTPND = 3) and a CELDIS factor of 0.1. 

The input parameters for the model simulation 

are summarized in Table 12.

The results of the analytical solution are 

compared graphically with those of MOC3D 

for three different planes in figs. 25-27. Figure 

25a shows the concentrations in the x-y plane 

of the point source as calculated using the 

analytical solution and figure 25b shows the 

same for the MOC3D results. In this view, the 

left edge of each map is a line of symmetry, so 

that the map represents only half of the true 

problem domain. Both sets of contours were 

generated for an identical number of points and 

locations to eliminate any differences 

attributable solely to the contouring procedure. 

The results agree very closely, although a 

slightly greater distance of migration or 

spreading is evident in the MOC3D results, 

both upstream as well as downstream of the 

source. However, a large part of this small 

difference can be explained simply by the fact 

that the source is applied over a larger area in 

the horizontal plane of the MOC3D model, in 

which the length of the source cell is 3 m in the 

direction parallel to flow.

46
Figure 26 shows a comparison of the 

results in a vertical plane parallel to the flow 

direction and aligned with column 2 of the 

numerical grid (and y = 0.75 m in the analytical 

solution). In this view, the top edge is a line of 

symmetry, so that the cross section represents 

only the lower half of the true problem domain. 

Overall the agreement is very close, although 

the numerical results show slightly more 

upstream dispersion, particularly for the lower 

concentrations (less than 10). Figure 27 shows 

a comparison of the results in a vertical plane 

transverse to the flow direction and aligned

Table 12. Parameters used in MOC3D simulation 

of transport from a continuous point source in a 

three-dimensional steady-state flow system

Parameter Value

= T

yy

PERLEN (length of stress 

period)

Source concentration (C") 

Q (at well) 

Source location

Number of rows 

Number of columns 

Number of layers 

DELR (Ax) 

DELC (Ay) 

Thickness (b)

NPTPND (Initial number 

of particles per cell)

CELDIS

INTRPL (Interpolation 

scheme)

0.0125 m2/day 

0.25 

0.6m 

0.03m 

0.006 m 

400 days

0.1 m/day 

0.0 m/day 

0.0 m/day 

2.5 x 106 g/m3 

1.0 x 10-6 m3/d

row 8, column 

1,layer 1

30

12

40

3m

0.5 m

0.05m

3

0.1 

1

with row 12 of the numerical grid (and x = 

34.5 m in the analytical solution). In this view, 

both the left and top edges are lines of 

symmetry, so that this cross section represents 

only one quarter of the true problem domain. 

Once again the overall agreement is excellent. 

The only noticeable difference is a minor one. 

In the numerical results there is slightly more 

lateral spreading, as indicated by a 

displacement of the contours for the lower 

values of concentration by less than half a cell 

distance away from the axis of symmetry, but 

only in the upper half of the grid.

10

100

15

20

25

30

(a) 

Analytical

30

10

Flow Direction

i i i i i i i i i i i 30

5 10 

Column

5 10 

Column

Vertical Exaggeration = 0.33

Figure 25. Concentration contours for (a) analytical and (b) numerical solutions in the horizontal 

plane containing the solute source (layer 1) for 

three-dimensional solute transport in a uniform 

steady flow field. Parameters are defined in 

Table 12.

47
10

20

30

40

(a) 

Analytical

Flow Direction

I i i i i i i i i

10 15 

Row

20 25 30

10

20

t

30

40

(b) 

MOC3D

Flow Direction

10 15 

Row

20 25 30 

Vertical Exaggeration « 33.5

Figure 26. Concentration contours for (a) analytical and (b) numerical solutions in the vertical 

plane parallel to the flow direction and aligned with column 2 for three-dimensional solute 

transport in a uniform steady flow field. Parameters are defined in Table 12.

10

20 i_

Q>

J3 30 

40

-100 (a) 

Analytical

6 7 

Column

10 11 12

10

_ 20

I J3 30

(b) 

MOC3D

40 L_L

6 7 

Column

8 9 10 11 12 

Vertical Exaggeration « 2.2

Figure 27. Concentration contours for (a) analytical and (b) numerical solutions in the vertical 

plane transverse to the flow direction and aligned with row 12 for three-dimensional solute 

transport in a uniform steady flow field. Parameters are defined in Table 12. Flow is towards the 

reader.

48
Two-Dimensional Radial Flow and 

Dispersion

A radial dispersion problem was used to 

compare the MOC3D solution to the analytical 

solution given by Hsieh (1986) for a finiteradius injection well in an infinite aquifer of 

two dimensions. The problem is equivalent to 

flow from a single injection well; the velocities 

vary in space and are inversely related to the 

distance from the injection well. The 

governing equation for the analytical solution is

Table 13. Parameters used in MOC3D simulation 

of two-dimensional, steady-state, radial flow case

dc_

dt

A^dC_ 

r dr

Ad2C 

r dr2

r>rw t>0 (74)

where A = Q/27tbe, r is the radial distance 

from the center of the well, rw is the radius of 

the injection well, a is the longitudinal 

dispersivity in radial flow, Q is the volumetric 

rate of the well injection, b is the thickness of 

the aquifer, and e is porosity. The initial and 

boundary conditions are

(75a)

C(rw,t) =

r>r

r>0

f>0.

(75b) 

(75c)

The radius of the well (rw) was set to 1.0 

(dimensionless) and the concentration of the 

injected tracer (C') was 1.0 (dimensionless) at 

the well.

The problem was modeled using a grid 

having 30 cells in the ^-direction and 30 cells in 

the y-direction, representing one quadrant of 

the radial flow field (90 of 360 degrees). The 

initial concentration was set to 1.0 at the well 

node (1,1), defined by a specified flux of 

56.25 m3 /h. The input parameters for the 

model simulation are summarized in Table 13. 

Initial particle positions were defined using the 

custom particle placement option in the input 

data set and were aligned in a quarter circle in 

two-degree increments, equidistant from the 

center point of radial symmetry in the upper left 

corner of the grid (see fig. 28). Using a large

Parameter Value

TXX 

£

UTV

PERLEN (length of stress 

period)

Q (at well)

Source concentration (C") 

Number of rows 

Number of columns 

Number of layers 

DELR (Ax) 

DELC (Ay) 

Thickness (b)

NPTPND (Initial number of 

particles per cell)

CELDIS

INTRPL (Interpolation 

scheme)

3.6 m2/hour 

0.2

10.0 m 

10.0m 

10.0m 

1000 hours

56.25 

m3/hour

1.0 

30 

30 

1

10.0 m 

10.0 m 

10.0 m 

46

0.5

2

number of particles per cell in this 

configuration ensures that all cells in the grid 

will include at least one particle pathline 

emanating from the source cell. If one of the 

MOC3D default options was selected to place 

fewer particles in a regular geometric pattern 

where strong flow divergence exists, then 

some cells located just a few rows or columns 

away from the source cell would never receive 

particles originating in the source cell. This 

would make it impossible to calculate 

accurately concentration changes caused by 

advection. This custom configuration also 

parallels the expected solute distribution pattern 

in a radially divergent flow field; standard 

particle positioning (linear, quadratic, or cubic 

configurations) results in extreme spreading of

49
particles at small distances from the source, 

distorting contours of concentration in that 

area. The same pattern of particles shown for 

the source cell in figure 28 is repeated in all 

cells of the grid at the start of the simulation.

Figure 29 shows a comparison of the 

concentrations calculated in one quadrant after 

1000 hours using both the analytical (a) and 

numerical (b) solutions. The two solutions are 

almost identical.

Point Initial Condition in Uniform 

Flow

A test problem for three-dimensional solute 

transport from an instantaneous point source, 

or Dirac initial condition, in a uniform flow 

field was used to test the MOC3D model. An 

analytical solution for an instantaneous point 

source in a homogeneous infinite aquifer is 

given by Wexler (1992), and he presents a 

code (named POINTS) for a related case for a 

continuous point source. The POINTS code 

was modified to solve for the desired case of an 

instantaneous point source. Test problems

-0.5

X

0.0

-0.5 0.5

0.0

0.5

(1,1)

EXPLANATION

 Node location

 Initial particle location

Figure 28. Initial particle positions within the 

source cell for radial flow case (based on 

custom particle placement and NPTPND = 46). 

The relative coordinates on the x- and y-axes 

shown for the cell (1,1) are the same for any 

cell of the grid; this relative coordinate system 

is used for the custom definition of particle 

locations in the input file.

10 -

20 -

30

10 20 

X (nodes)

to

CD 

TJ 

O

c

10

20

30

(b) 

MOC3D

-ri i i i i i i i I i i i i i i i i i I i i i i i i i i i

10 20 

X (nodes)

30

Figure 29. Contours of relative concentrations calculated using (a) analytical and (b) numerical 

models for solute transport in a steady radial flow field. Source concentration is 1.0 and source is 

located in cell (1,1). Grid spacing is 10.0 m.

50
were designed to evaluate the numerical 

solution for two cases one in which flow is 

parallel to the grid (in the ^-direction) and one 

in which flow occurs at 45 degrees to the xand y-axes. This allows us to evaluate the 

accuracy of the numerical model for this basic 

type of problem, and also to evaluate the

sensitivity of the numerical solution to the 

orientation of the flow relative to that of the 

grid.

The governing equation and boundary 

conditions for an instantaneous point source are 

(seeWexler, 1992, p. 42):

dc_ 

dt

WLC

where V is the velocity in the direction of flow 

(assumed to be the ^-direction in eq. 76), K is 

the decay rate (K = 0 for this problem), Q is the 

injection rate for the well, 8 is a dirac delta 

function, Xc, Yc, and Zc are coordinates of the 

point source, and t' is the time at which the 

instantaneous point source activates. 

The initial condition (at t' = 0) is:

C = 0,

 oo < X < oo, oo < _y < +00, oo < i < +00. (78)

For the test case of flow in the ^-direction, we 

assumed Vx = 1.0 m/d, and Vy = Vz = 0.0 

m/d. For flow at 45 degrees to x and y, we 

assumed Vx = Vy = 1.0 m/d, and Vz = 0.0 

m/d. For both cases, the distance the plume 

travels in the ^-direction is the same for equal 

simulation times. Note, however, that the 

magnitude of velocity is higher in the latter 

case; therefore, there will be more dispersion in 

that problem during an equivalent time interval.

Several different grid spacings were used 

for the model simulations to help show the 

relation between discretization and the accuracy 

of the numerical results. The coarsest transport

d2C ,.dC

-Yc )S(z-Zc )S(t-t')

c=o,

c=o,

c=o,

X = ±oo

y = ±00

Z = ±oo

(76)

(77a) 

(77b) 

(77c)

subgrid used 24 rows, 24 columns, and 24 

layers with a grid spacing of 10.0 m in each 

direction. Subsequent runs doubled and tripled 

the number of rows and columns; the grid 

spacing was reduced accordingly so that the 

model domain (the aquifer volume being 

simulated) remained the same for all grids, and 

changes in accuracy would be attributable only 

to changes in the spatial discretization. 

Because advective transport occurs only in the 

x-y plane, the number of layers was held 

constant at 24. The input parameters for the 

simulations are presented in Table 14, which 

includes the several different values used for 

grid dimensions and spacing.

The results for both the analytical and 

numerical solutions for the case in which flow 

occurs only in the ^-direction are shown in 

figure 30. The MOC3D results for the coarsest 

grid (fig. 30b) clearly show too much 

spreading transverse to flow (that is, in the ydirection) relative to the analytical solution (fig. 

30a). Note that the analytical solution was 

contoured for values only at the same exact 

locations as the nodes in the grid used for the 

numerical solution to which it is being

51
compared (fig. 30b); this eliminates differences 

for a visual comparison due only to artifacts of 

the contouring procedure, and is the reason that 

the analytical solution appears less smooth than 

it should. To improve numerical and 

contouring accuracy, a finer grid was used for 

both the analytical and MOC3D solutions (figs. 

30c and 30d). Using three times as many 

nodes in each direction in the horizontal plane, 

the numerical dispersion is reduced 

significantly and the numerical solution very 

closely matches the analytical solution. (Note 

that the analytical solution mapped in fig. 30c is 

based on values calculated at nine times as 

many points as that in fig. 30a.)

The results of the test problem for flow at 

45 degrees to the grid are shown in figure 31. 

The analytical solution (fig. 3la), which 

provides the basis for the evaluation, was 

solved on a 72x72 grid, and the MOC3D 

solutions are shown for a 24x24x24 grid (fig. 

31b), a 48x48x24 grid (fig. 31c), and a 

72x72x24 grid (fig. 31d). Unlike the 

previous case (where flow is aligned with the 

grid), the numerical results in figure 31 show a 

noticeable difference in the shape of the plume 

relative to the analytical solution. The 

numerically calculated "hourglass" shape is 

characteristic of a grid-orientation effect and is 

related to the cross-product terms of the 

dispersion tensor. When flow is oriented 

parallel to the grid, or when longitudinal and 

transverse dispersivities are equal, the crossproduct terms of the dispersion equation are 

zero. Because flow is at 45 degrees to the grid, 

the cross-product terms of the dispersion 

equation are nonzero. The model estimates the 

concentration gradients associated with the 

cross-product terms less accurately than those 

associated with the diagonal terms, and 

therefore the overall solution is less accurate.

The magnitude of this effect is minimized 

by using a finer grid. Overall, the coarsest grid 

exhibits too much spreading, but the next finer 

grid results in minimal numerical dispersion

(although the grid-orientation effect is not 

eliminated). Further reducing the grid spacing 

(fig. 3 Id) does not significantly further reduce 

numerical dispersion, thereby indicating the 

desired characteristic of grid convergence.

Table 14. Parameters used in MOC3D simulation 

of three-dimensional transport from a point source 

with flow in the x-direction and flow at 45 degrees 

to x and y

Parameter Value

e 

OCL

PERLEN (length of 

stress period)

Vx

Vy

Source concentration (C')

Source location

Number of rows 

Number of columns 

Number of layers 

DELR (Ax)

DELC (Ay}

Thickness (b)

NPTPND (Initial number 

of particles per cell)

CELDIS

INTRPL (Interpolation 

scheme)

10.0 m2/day 

0.1 

1.0m 

0.1 m 

O.lm 

90 days

1.0 m/day 

0.0 m/day* 

0.0 m/day 

1 x 106

x = 30 m, 

y = 120 m, 

z = 40 m**

24, 48, and 72 

24, 48, and 72 

24

10.0 m, 5.0 m, 

and 3.33 m

10.0 m, 5.0 m, 

and 3.33 m

10.0m 

8

0.5

2

* For flow at 45 degrees to x and y, Vy = 1.0 m/day 

** For flow at 45 degrees to x and y, the source 

location is x = 30 m, y = 30 m, z = 120 m.

52
Each of the MOC3D results also shows a 

slight asymmetry in the shape of the plume in 

the direction of flow (that is, there is slightly 

less forward spreading compared to backward 

spreading), which is inconsistent with 

symmetrical spreading indicated by the

analytical solution. This is caused by the 

sequence in which the dispersive and advective 

terms of the transport equation are solved.

The numerical errors are exaggerated in 

figures 30 and 31 because the concentrations 

are contoured on a logarithmic scale. Some of

10

20

Initial location of point source

(a) Analytical (calculated 

at 24x24 points)

10

20

(b) MOC3D (24x24 grid)

10

X (nodes)

20 10 20

X (nodes)

CO 

CD -a 

o

20

40

60 (c) Analytical (calculated 

at 72x72 points)

20

40

60

(d) MOC3D (72x72 grid)

!_________I_________|_________I

20 40

X (nodes)

60 20 40

X (nodes)

60

Figure 30. Concentration contours for (a, c) analytical and (b, d) numerical solutions for transport of 

a point initial condition in uniform flow in the x-direction. The z-component of flow is zero, but there 

is dispersion in all three directions. Contour values are the log of the concentrations.

53
20

40

60

Initial location of point source

(a) Analytical (calculated 

at 72x72 points)

10

20

(b)MOC3D (24x24 grid)

20 40

X (nodes)

60 10

X (nodes)

20

eg 20h

 o

O

c

40 -

(c) MOC3D (48x48 grid)

20

40

60

(d) MOC3D (72x72 grid)

20

X (nodes)

40 20 40

X (nodes)

60

Figure 31. Concentration contours for (a) analytical and (b, c, d) numerical solutions for transport of 

a point initial condition in uniform flow at 45 degrees to xand y. Contour values are the log of the 

concentrations.

the discrepancies may also be due to the 

contouring program used to visually represent 

the solutions.

When the flow is at an angle to the grid, as 

for the case illustrated in fig. 31, then negative 

concentrations are most likely to occur. In this 

case, some small areas of slightly negative 

concentrations were calculated, but are not

evident in fig. 31 because they were filtered out 

during the contouring process to allow a clear 

depiction of the position of the plume. 

However, to indicate the extent of the area of 

negative concentrations, we have replotted the 

central part of the domain illustrated in fig. 3 Ib 

(for the 24x24 grid) in fig. 32, in which all 

areas where the relative concentration is less

54
5.0

7.5

10.0

12.5

15.0

17.5

5.0 7.5 10.0 12.5 15.0 17.5

5.0 

7.5 

10.0 

1Z5 

15.0 

17.5

5.0 7.5 10.0 12.5 15.0 17.5

5.0 

7.5 

10.0 

1Z5 

15.0 

17.5

5.0 7.5 10.0 12.5 15.0 17.5

Area of negative concentration

Figure 32. Concentration contours showing effects on areas of negative concentrations related 

to decreasing CELDIS factor in MOC3D in simulation of flow at 45 degrees to grid having 24 rows 

and 24 columns of nodes: (a) CELDIS = 0.50; (b) CELDIS = 0.25; and (c) CELDIS = 0.10.

than -0.05 are shaded. Figure 32a represents 

the same solution as shown in fig. 31b. We 

tested the sensitivity of the solution and of the 

extent of negative concentrations to the size of 

the transport time increment by adjusting the 

value of CELDIS. The area in which negative 

concentrations occurred at the same elapsed 

simulation time was slightly smaller for the 

smallest value of CELDIS. In all three cases 

the mass balance errors were about the same 

and always less than 0.2 percent relative to the 

initial solute mass stored.

Constant Source in Nonuniform 

Flow

Burnett and Frind (1987) used a numerical 

model to analyze a hypothetical problem having 

a constant source of solute in a finite area at the 

surface of an aquifer having homogeneous 

properties, but nonuniform boundary 

conditions, which result in nonuniform flow. 

Because an analytical solution is not available 

for such a complex system, we use their results 

for this test case as a benchmark for comparison with the results of applying MOC3D to 

the same problem. Burnett and Frind (1987)

used an alternating-direction Galerkin finiteelement technique to solve the solute-transport 

equation in both two and three dimensions. 

Their model also includes the capability to vary 

0(T as a function of direction, thereby allowing 

that feature of MOC3D to be evaluated in the 

same problem set.

A simplified diagram of the problem is 

illustrated in figure 33. The left (x = 0) and 

bottom surfaces are no-flow boundaries, 

representing a ground-water divide and an 

impermeable base, respectively. The top and 

right surfaces are constant-head boundaries, 

representing the water table and a discharge 

boundary, respectively. The front and back 

vertical faces are no-flow boundaries, 

representing streamlines or flow paths parallel 

to those surfaces. The length of the domain is 

200 m and the saturated thickness varies from 

21 m on the left no-flow boundary to 20 m on 

the right constant-head boundary. The heads 

on the upper surface are specified as a onequarter cosine from 1 m on the left to 0 m on 

the right. Heads are fixed at a slightly negative 

value (-0.00736 m) at all nodes in the right 

column of cells, so that the head on the right 

side of the transport domain (at x = 200 m) will 

almost exactly equal 0 m, and the heads and

55
CO 

CO

LLI

o

I

20m

200m

Figure 33. Transport domain and boundary conditions for nonuniform-flow test problem of 

Burnett and Frind (1987); front surface represents a plane of symmetry (modified from Burnett 

and Frind, 1987, fig. 4).

hydraulic gradients throughout the domain 

calculated using MOC3D would be consistent 

with those in the analysis of Burnett and Frind 

(1987). The aquifer properties are assumed to 

be homogeneous and isotropic. These 

boundary conditions yield a two-dimensional 

flow field, which has components of flow in 

the jc- and z-directions only. Therefore, the 

specified width of the domain is varied for each 

particular simulation to accommodate the 

plume, and depends on the dimensionality of 

the simulation (whether two- or threedimensions) and on the values of the 

dispersivity coefficients. The solute source is 

located between 18.25 and 32.50 m from the 

left side of the aquifer and has a width of 10 m, 

extending 5 m on either side of the plane of 

symmetry at y = 0.

Cases of both two- and three-dimensional 

transport were examined for this basic 

problem. The grids used in the MOC3D 

simulations were designed to maximize 

compatibility with the results of the finiteelement models used by Burnett and Frind 

(1987), so that comparisons of results would 

represent a reasonable benchmarking exercise. 

However, some differences in discretization 

could not be avoided because the finite-element 

methods allow specifications of values at 

nodes, which can be placed directly on

boundaries. Nodes in MOC3D are located at 

the centers of cells, and values specified at 

nodes are always one-half of the grid spacing 

away from the edge of the model domain. 

Among the small differences arising from the 

alternative discretization schemes is that (1) the 

modeled location of the 14.25 m long source 

area is offset by 0.225 m towards the right in 

the MOC3D grid, and (2) the total length of the 

domain is 199.5 m in the MOC3D grid.

The first analysis of this test case focused 

on the simplest one a two-dimensional 

analysis. The input data values for this 

analysis are listed in Table 15. The MOC3D 

grid consisted of 141 columns and 91 layers. 

However, the top layer and right column of 

cells are devoted to assuring consistent 

boundary conditions between the two models, 

and they were considered to lie outside of the 

domain of the transport problem for 

benchmarking purposes. That is, 140 columns 

at a spacing of 1.425 m yields a transport 

domain length of 199.5 m, which is 

approximately equal to the desired total 

horizontal distance of 200 m. Similarly, the 

appropriate height of the domain is assured by 

setting the thickness of each cell within a given 

column equal to l/90th of the height of the 

column, where the height (or total saturated 

thickness) varies from 21 m on the left to 20 m

56
on the right. (In comparison, Burnett and 

Frind used a variable spacing in which Ax 

ranged from 2.5 to 6.0 m and Az ranged from 

0.75 to 1.25 m.) The top flow layer consisted 

of constant-head nodes and the solute source. 

Because of the symmetry in the flow field, we 

were able to increase the efficiency of the 

simulation by using a custom initial particle 

placement of only three particles in each cell, as 

shown in figure 34, and still achieve reasonably accurate results. Burnett and Frind (1987) 

report that their solution yielded an average

Table 15. Parameters used in MOC3D simulation of transport in a vertical plane from a 

continuous point source in a nonuniform, 

steady-state, two-dimensional flow system 

described by Burnett and Frind (1987)

-4

Parameter Value

K 

e

Dm

PERLEN (length of 

stress period)

Source concentration (C')

Number of rows 1 

Number of columns 

Number of layers 1 

DELR (Ax) 

DELC (Ay) 

Thickness (b)

NPTPND (Initial number 

of particles per cell)

CELDIS

INTRPL (Interpolation 

scheme)

1.0 m/day 

0.35 

3m 

0.01 m 

10-4 m2/day 

12,000 days

1.0

141

1

91

1.425 m

1.0m

0.2222-0.2333 m

3

1.0 

1

1 One row and layer were allocated to defining 

boundary conditions, so concentrations calculated in 

only 140 rows and 90 layers were used for 

benchmarking.

areal recharge rate of about 8x10 m/day; the 

MOC3D solution yielded an average areal

-4 recharge rate of about 7.8 x 10 m/day. This 

agreement is evidence that the parameters and 

boundary conditions for the two model 

analyses are similar enough to permit a 

benchmarking comparison.

-0.5

o.o

(1,1)

0.5 

-0.5 0.0 0.5 

X

EXPLANATION 

Node location 

Initial particle location

Figure 34. Initial particle positions within a cell 

for the Burnett-Frind (1987) test case (based 

on custom particle placement and NPTPND = 

3). The relative coordinates on the x- and zaxes shown for the cell (1,1) are the same for 

any cell of the grid; this relative coordinate 

system is used for the custom definition of 

particle locations in the input file.

Results for the two-dimensional case from 

the MOC3D model closely match those of 

Burnett and Frind (see fig. 35). Both models 

represent the solution on the plane of symmetry 

(that is, on the front face of the block shown in 

fig. 33). The concentration contours are 

located in almost exactly the same positions for 

both models. However, in the MOC3D 

results, the contours lag slightly behind those 

of Burnett and Frind (1987). This may be

57
80

2D Finite-Element Model -,

120 160 200

(b) MOC3D Model

120 160 200

Figure 35. Two-dimensional simulation results for nonuniform-flow test case 

showing plume positions as contours of relative concentration: (a) finiteelement model (modified from Burnett and Frind, 1987, fig. 8a), and (b) 

MOC3D. Contour interval is 0.2 relative concentration.

attributable to small differences in the numerical 

treatment of the source between the two models 

and (or) to the slightly lower flux (and velocity) 

in the MOC3D solution. Note that the solution 

to the two-dimensional case in the x-z plane 

does not depend on the value of &THFor the three-dimensional analyses, the 

MOC3D grid is expanded to 15 rows having 

Ay of 1.0 m. The source is applied over the 

first 5 rows, taking advantage of the symmetry 

along the j-axis to account for the 10 m width

of the source. Note that because of symmetry, 

the flow fields are identical in the two- and 

three-dimensional cases. Figure 36 shows the 

transport results for both models for the case in 

which (%TV =0.01 m and O,TH = 0.1 m. In the 

MOC3D results (fig. 36b) the vertical plane in 

the first row is contoured. Note that there will 

be a slight discrepancy in the basis of 

comparison because concentrations from 

MOC3D are evaluated at the center of the block 

(1/2 of a cell width from the plane of

20

10

(a) 3D Finite-Element Model

40 80 120 160 

(b) MOC3D Model

200

40 80 120 160 200

Figure 36. Three-dimensional simulation results for nonuniform-flow test case 

in which ajH = 0.1 m and ajv = 0.01 m: (a) finite-element model (modified 

from Burnett and Frind, 1987, fig. 8c), and (b) MOC3D. Plume positions are 

represented by contours of relative concentration; contour interval is 0.2 

relative concentration.

58
symmetry), whereas those from Burnett and 

Frind (1987) are evaluated on the cell faces 

(directly on the plane of symmetry). A threedimensional visualization of the MOC3D 

results are presented in fig. 37, which was 

generated from the concentrations on particles 

(as opposed to nodal concentrations in fig. 

36b) using a three-dimensional visualization 

software package. (The plume definition from 

particle concentrations and locations is 

inherently more precise than can be defined 

from averaged nodal concentrations.) This 

perspective shows more clearly the entire 

plume and the magnitude of lateral spreading of 

the plume. It also illustrates the impact of 

specifying unequal transverse dispersivity 

coefficients, in that the transverse spreading in 

the horizontal direction increases more rapidly 

than does the transverse spreading in the

vertical direction, as evident in the "chair" cut 

near the source location and on the face of the 

plume at the downgradient boundary. Figure 

38 shows the results for the case in which the 

vertical transverse dispersivity is increased by a 

factor of ten, so that CCTH = arv = 0.1 inOverall, the MOC3D results (figs. 35b, 36b 

and 38b) agree closely with those of Burnett 

and Frind (1987) (figs. 35a, 36a and 38a).

Comparison of the three-dimensional 

results with the two-dimensional analysis 

shows that all contours are closer to the source 

of solute in the three-dimensional cases. This 

is expected because the contaminant source has 

a finite length and consideration of the 

additional dimension for the dispersion process 

allows spreading of solute in the ^-direction, 

which means that less solute mass will remain 

in the vertical plane being contoured. This

Figure 37. Perspective view of MOC3D results for three-dimensional problem of constant source 

in nonuniform flow and unequal transverse dispersivity coefficients. This visualization of the 

plume was generated from particle concentrations using a three-dimensional visualization 

software package and is derived from the same simulation that is the basis of fig. 36b. Note that a 

piece of the plume near the source is cut away (a "chair" cut) to expose a clearer view of the 

degree of transverse spreading in the selected vertical and horizontal planes. Shading 

increments are in relative-concentration intervals of 0.20, and interval bounds range between 0.9 

and 0.1. Concentrations less than 0.10 are transparent.

59
40

(a) 3D Finite-Element Model ,

160

(b) MOC3D Model

200

80 120 160 200

Figure 38. Three-dimensional simulation results for nonuniform-flow test case in 

which CCJH = aTv = 0.01 m: (a) finite-element model (modified from Burnett and 

Frind, 1987, fig. 9b), and (b) MOC3D. Plume positions are represented by 

contours of relative concentration; contour interval is 0.2 relative concentration.

makes it appear that the plume has not spread 

as far in the ^-direction during the same elapsed 

time. In the three-dimensional simulations, the 

lower-value and higher-value contours are 

slightly closer to the 0.5 contour in the MOC3D 

results, and the contours defining the lateral 

edges of the plume are spaced closer together. 

These characteristics indicate that the MOC3D 

results may include less numerical dispersion 

than the finite-element results. However, this 

minor difference may simply be an outcome of 

having used a finer grid spacing in MOC3D.

Relative Computational Efficiency

The computational effort required by the 

MOC3D code is strongly dependent on the size 

of the problem being solved, as reflected 

primarily by the total number of nodes, total 

number of particles, and total number of time 

increments. Analyses indicate that the greatest 

computational effort, as measured by CPU 

time, is typically expended in the particle 

tracking

routines. For a given problem, the efficiency 

of the code may vary significantly as a function 

of the characteristics of the particular computer 

on which the simulation is performed, and to 

some extent on which FORTRAN compiler 

(and which compiler options) were used to 

generate the executable code.

To provide a qualitative indication of these 

relations, we have run all of the sample 

problems described in this report on a variety 

of computers. The relative running times for 

each problem on a variety of different 

computers are presented in table 16. The run 

times are measured as CPU time in seconds. 

As indicated in table 16, the efficiency for a 

given problem may vary by more than a factor 

of ten, depending on which of the tested 

computers were used. However, for the given 

test problems, the efficiency was much more 

sensitive to the overall size of the problem, and 

the CPU time on a given computer varied by 

about four orders of magnitude between the 

simple one-dimensional problem and the more 

complex three-dimensional problem. The 

model user should be aware that in some cases, 

model efficiency may be a serious constraint.

60
Table 16. Comparison of MOC3D simulation times for selected test cases on various co

o\

Run Time in CPU-secon

Problem 

Description

One-Dimensional

Steady Flow

Number of 

Nodes

120

Number of 

Moves

201

Maximum 

Number of 

Particles

360

Silicon 

Graphics1

0.71

Data 

General 

Server2

2.6

Data 

General 

Workstation 3

5.2

PC 

(90 MH

Pentium

3.5

Three-Dimensional 14,400 

Steady Flow

Two-Dimensional 900 

Radial Flow and 

Dispersion

Point Initial Condition 13,824 

in Uniform Flow (flow 

at 45 degrees to grid)

Constant Source in 12,831 

Nonuniform Flow 

(Two-Dimensional)

Constant Source in 192,465 

Nonuniform Flow 

(Three-Dimensional)

207 43,200 82 350 650

596 52,666 240 750 1,450

19 110,729 35 240 450

320

620

140

4,218 68,027 2,290 13,450 18,310 9,060

4,218 881,260 51,720 139,300 NT NT

NT = Not Tested

1 Silicon Graphics server with an R6000 chip running Irix 6.0.1 with 576MB RAM and a 90 MHz processor. MIPSpro F77 wa

2 Data General server with a Motorola 88110 chip running DG Unix 5.4R3.10 with 256MB RAM and a 45 MHz processor. Gr

FORTRAN-88000 was used to compile MOC3D.

3 Data General AVIION 530 running DG Unix 5.4R3.10 with 32MB RAM and a 33MHz processor. Green Hills Software FOR

compile MOC3D.

4 IBM-compatible Pentium PC running MS-DOS 6.2 with 8MB RAM and a 90 MHz processor. Powerstation 1.0 FORTRAN w

5 IBM-compatible Pentium PC running Windows95 with 32MB RAM and a 133 MHz processor. Powerstation 1.0 FORTRAN 

6 Macintosh 9500 PC with 120 MHz PowerPC 604 processor and 32MB of RAM. MOC3D was compiled using Fortner Rese

for Power Mac.
CONCLUSIONS

The MOC3D model described in this report 

can simulate the transient, three-dimensional, 

transport and dispersion of a solute subject to 

decay and retardation. The solute-transport 

model is integrated with the MODFLOW 

model, which is used to solve the ground-water 

flow equation for either steady-state or transient 

flow. The numerical methods used to solve the 

governing equations allow their application to 

systems having heterogeneous properties and 

complex boundary conditions. The package 

thus has broad general application and 

flexibility for application to a wide range of 

hydrogeological problems.

The accuracy and precision of the numerical 

results were tested and evaluated by 

comparison of the MOC3D results with 

analytical solutions for several relatively simple 

and idealized problems and by benchmarking 

comparisons against the results of other 

numerical codes for more complex problems 

for which no analytical solutions are available. 

These tests indicate that the model can 

successfully and accurately simulate the threedimensional transport and dispersion of a 

solute in flowing ground water. This 

implementation of the method of characteristics 

is not strictly mass conservative and the method 

of calculating a solute mass balance is 

inherently an approximation. Therefore, 

calculated mass-balance errors may be nonzero, 

but are generally less than 10 percent and often 

decrease and stabilize with time. For some 

problems, the accuracy and precision of the 

numerical results may be sensitive to the initial 

number of particles placed in each cell. The 

efficiency of the solution is sensitive to the total 

number of particles used and to the size of the 

transport time increment, as determined by the 

stability criteria for the solute-transport 

equation. An advantage of this method is that, 

in general, its accuracy and efficiency are 

greatest for advection-dominated problems,

which is a characteristic of many ground-water 

contamination problems that pose serious 

environmental risks.

REFERENCES

Appel, C.A., Reilly, T.E., 1994, Summary of 

computer programs produced by the U.S. 

Geological Survey for simulation of groundwater flow and quality, 1994: U.S. Geological 

Survey Circular 1104, 98 p.

Bear, Jacob, 1979, Hydraulics of Groundwater: 

McGraw-Hill, New York, 567 p.

Bredehoeft, J.D., and Finder, G.F., 1973, Mass 

transport in flowing groundwater: Water 

Resources Research, v. 9, no. 1, p. 194-210.

Burnett, R.D., and Frind, E.G., 1987, Simulation of 

contaminant transport in three dimensions, 2. 

Dimensionality effects: Water Resources 

Research, v. 23, no. 4, p. 695-705.

Domenico, P.A., and Schwartz, F.W., 1990, Physical 

and Chemical Hydrogeology: John Wiley & 

Sons, New York, 824 p.

Engineering Technologies Associates, Inc., 1989, Three 

dimensional solute transport using the randomwalk algorithm and MODFLOW, User's manual 

for RAND3D: Engineering Technologies 

Associates, Inc., Ellicott City, MD, 102 p.

Freeze, R.A., and Cherry, J.A., 1979, Groundwater: 

Prentice-Hall, Englewood Cliffs, N.J., 588 p.

Garabedian, S.P., LeBlanc, D.R., Gelhar, L.W., and 

Celia, M.A., 1991, Large-scale natural gradient 

tracer test in sand and gravel, Cape Cod, 

Massachusetts, 2, Analysis of spatial moments 

for a nonreactive tracer: Water Resources 

Research, v. 27, no. 5, p. 911-924.

Garder, A.O., Peaceman, D.W., Pozzi, A.L., 1964, 

Numerical calculation of multidimensional 

miscible displacement by the method of 

characteristics: Soc. Petroleum Eng. Jour., v. 4, 

no.l, p. 26-36.

Gelhar, L.W., Welty, C., and Rehfeldt, K.R., 1992, A 

critical review of data on field-scale dispersion in 

aquifers: Water Resources Research, v. 28, no. 

7, p. 1955-1974.

Goode, D.J., 1990a, Particle velocity interpolation in 

block-centered finite difference groundwater flow 

models: Water Resources Research, v. 26, no. 5, 

p. 925-940.

62
Goode, D.J., 1990b, Governing equations and model 

approximation errors associated with the effects 

of fluid-storage transients on solute transport in 

aquifers: U.S. Geological Survey WaterResources Investigations Report 90-4156, 20 p.

Goode, D.G., and Appel, C.A., 1992, Finite-difference 

interblock transmissivity for unconfined aquifers 

and for aquifers having smoothly varying 

transmissivity: U.S. Geological Survey WaterResources Investigations Report 92-4124, 79 p.

Goode, D.J., and Konikow, L.F., 1989, Modification of 

a method-of-characteristics solute-transport model 

to incorporate decay and equilibrium-controlled 

sorption or ion exchange: U.S. Geological 

Survey Water-Resources Investigations Report

89-4030, 65 p.

Goode, D.J., and Konikow, L.F., 1991, Testing a 

method-of-characteristics model of threedimensional solute transport in ground water: in 

Lennon, G. P., ed., Symposium on Ground 

Water, Proceedings of the International 

Symposium, Nashville, Tenn., American 

Society of Civil Engineers, New York, p. 21-27.

Harbaugh, A.W., and McDonald, M.G., 1996a, User's 

documentation for MODFLOW-96, an update to 

the U.S. Geological Survey modular finitedifference ground-water flow model: U.S. 

Geological Survey Open-File Report 96-485, 56

PHarbaugh, A.W., and McDonald, M.G., 1996b, 

Programmer's documentation for MODFLOW96, an update to the U.S. Geological Survey 

modular finite-difference ground-water flow 

model: U.S. Geological Survey Open-File 

Report 96-486, 220 p.

Hill, M.C., 1990, Preconditioned Conjugate-Gradient 2 

(PCG2) A computer program for solving 

ground-water flow equations: U.S. Geological 

Survey Water-Resources Investigations Report

90-4048, 43 p.

Hsieh, P. A., 1986, A new formula for the analytical 

solution of the radial dispersion problem: Water 

Resources Research, v. 22, no. 11, p. 1597- 

1605.

Konikow, L.F., and Bredehoeft, J.D., 1978, Computer 

model of two-dimensional solute transport and 

dispersion in ground water: U.S. Geological 

Survey Techniques of Water-Resources 

Investigations, Book 7, Chapter C2, 90 p.

Konikow, L.F., and Grove, D.B., 1977, Derivation of 

equations describing solute transport in ground 

water: U.S. Geological Survey Water-Resources

Investigations Report 77-19, [Revised 1984], 30 

PMcDonald, M.G., and Harbaugh, A.W., 1988, A 

modular three-dimensional finite-difference 

ground-water flow model: U.S. Geological 

Survey Techniques of Water-Resources 

Investigations, Book 6, Chapter Al, 586 p.

Pinder, G.F., and Cooper, H.H., Jr., 1970, A numerical 

technique for calculating the transient position of 

the saltwater front: Water Resources Research, 

v. 6, no. 3 , p. 875-882.

Pinder, G.F., and Gray, W.G., 1977, Finite Element 

Simulation in Surface and Subsurface Hydrology: 

Academic Press, New York, N.Y, 288 p.

Prickett, T.A., Naymik, T.G., and Lonnquist, C.G., 

1981, A "random-walk" solute transport model 

for selected groundwater quality evaluations: 

Illinois State Water Survey Bulletin 65, 103 p.

Reddell, D.L., and Sunada, O.K., 1970, Numerical 

simulation of dispersion in groundwater aquifers: 

Colorado State University, Ft. Collins, CO, 

Hydrology Paper 41, 79 p.

Robson, S.G., 1974, Feasibility of digital water-quality 

modeling illustrated by application at Barstow, 

California: U.S. Geological Survey WaterResources Investigations Report 46-73, 66 p.

Robson, S.G., 1978, Application of digital profile 

modeling techniques to ground-water solutetransport at Barstow, California: U.S. 

Geological Survey Water-Supply Paper 2050, 28 

PScheidegger, A.E., 1961, General theory of dispersion 

in porous media: Jour. Geophys. Research, v. 

66, no. 10, p. 3273-3278.

Wexler, E.J., 1992, Analytical solutions for one-, two-, 

and three-dimensional solute transport in groundwater systems with uniform flow: U.S. 

Geological Survey Techniques of WaterResources Investigations, Book 3, Chapter B7, 

190 p.

Zheng, C., 1990, MT3D: A modular three-dimensional 

transport model: S.S. Papadopulos and 

Associates, Inc., Bethesda, MD.

63
APPENDIX A: FINITE-DIFFERENCE APPROXIMATIONS

Applying finite-difference approximations that are centered-in-space and explicit (forwardin-time), the component of the dispersive flux in the ^-direction across the cell face at (j+l/2,i,k) 

(equivalent to My(7-+1/2,t,£)/^ Ay from eq. 57) may be written:

= -(£bD\

\ xx /7+l/2,f, Ax

[CU^'xy)j+l/2,i,

S^i* S^i*

l,i+l,k ~ Cj,i-\,k ~ Cj+l,i-l,

4Ay

-(fifrl>_y+1 -

V xz 'j+\/2,i,k 2

t+\

(Al)

where 2Bjj i> jc = bj^k + 1/2 (bj^k-1 + bj^+i) is the vertical distance between nodes (j,i,k+l) and 

(j,i,k-l). The superscript "*" indicates the use of an average concentration, as defined in eq. 55.

Similarly, the y-component of the dispersive flux vector atj,i+l/2,k is approximated by:

dC i \t+i 

= -(£bDyy )_ M/2

I \f+l 

-K*l>l/2,

Ay

4Ax

t+\ ]_

2

. (A2)

The z-component does not include the saturated thickness and is approximated by:

N *

^ dC}

3m -»

oxt

[^zxjj^k+m

S~t*

Cj+\,i,k

~*"

4Ay

(

c* + c* c* c* 

J ' '______J' '________J* '______J* '

4Ay

(A3)

Applying centered finite-difference approximations, the change in concentration due to dispersion, 

neglecting the sink/source term, can be written:

64
dC 

dr

A*

'+ (c*+m -C* )- /+l/2,aV J+l' l' K J' l' k )

2AxAy

(ebD, t+i * 

+

'

t+1

,,

i

t+l

2Bj,i,k

+

i \r-t-i

/ *

j,i,k

65
+ I* *~>* /I*

j+l,i,k ~ Cj-l,i,k+l ~ ~/-l

+ -

( \£+l i * * * * 

Vlj,i,k-\/2\ i'i+l- k + CJ,i+lk-l ~ Cj,i-\,k ~ Cj,i-\,k-\ (A4)

where the dispersion coefficient terms ebD and eD are given by:

v2

\D(i\j+\/2,i,k

(btfz(j+l/2,i,k)

IH+i/2,a

(A5)

-V2,i,k f

\hn\ H I/?/? \ OCl\j-\/2,i,k \D(l\j-\/2,i,k

ebD

-aa

= a - a

\j+l/2,i,k

 \

- aT (]f , rv(^); 

b4x(j-l/2,i,k)b4z(j-l/2,i,k) !J L^-J: ^ J-:-jL

(A7)

(A8)

(A9)

(A10)

/ , n \

\ebD\ V yy)j,i+i/2

=a

,k

a v

(A12)

66
-a<*TH (k)

(ebD\ = (aL(k) - aT )^zl^d^ (A14)

V yx >j,i-\/2,k V L(k) TH(k)j \ h I v > D

ebD\ *>J.i - aT 5 r

+aTv (k+l/2)

4y(j,i, (XTv (k-U2) + aL(k-l/2)

-aa

where all terms are at time level t+l. The flux terms normal to the finite-difference block (or cell 

faces) are known directly from the solution to the flow equation. For the horizontal terms, bq, the 

volumetric flux per unit width is known, whereas the vertical flux is specific discharge, or 

volumetric flux per unit area.

To compute the dispersion coefficients at a block interface, the flux must be computed at 

this location. The flux normal to the block face is known from the finite-difference solution of the 

flow equation. However, the other two components must be interpolated from nearby values. 

Horizontal fluxes are averaged for x fluxes at y block interfaces, and vice versa, by:

67
L

t}4xi x(j,i+V2,k)

, _ b4y(j,i-l/2,k) + b4y(j,i+\/2,k) + b4y( j+\,i-\/2,k) + b4y(j+\,i+l/2,k) t>4y(j+l/2,i,k)- ^ 

Horizontal fluxes per unit width are averaged and normalized by corresponding layer thicknesses 

to convert to fluxes per unit area at layer interfaces:

_ bclx(j-\/2,i,k-\) +bclx(j+\/2,i,k-\) + . bclx(j-\/2,i,k) +bc - lx(j+\/2,i,k) f \^n\ (A27)

_ -t A( j i/^.,i,K,-t-is____-t A( y-ri/z,,t> A,-riy , "x(j \/2,i,k)____-^A; JTI/^.,I,K.J /AOQ\ ^f;,a+l/2J - T7 + 77 CAZ8)

_ b4y(j,i-\/2,k+\) +b(ly(j,i+V2,k+l)

 

b4y(j,i-l/2,k) + b4y(j,i+l/2,k)

Vertical fluxes per unit area are averaged and multiplied by corresponding layer thicknesses to 

convert to fluxes per unit width at row and column interfaces:

'

b4z(j,i+l/2,k) = ~~£~[4z(j,i,k-l/2) + Qz(j,i,k+l/2)\ + '4 '

i / / i k r

(j,i,k-l/2) + Qz(j,i,k+\/2)\ + ' ' \4z(j,i-\,k-\/2)

7 i k F 1 7+1 i k

b4z(j+l/2,i,k) = ~

i / i / ^ r i

J + "[^zO-1,1^-1/2) + ^z(;-l,i,Jfc+l/2) J68
The magnitudes of fluid flux at block interfaces are:

i/2

l/2 r 2 / \2 \2 1 \bq\ 1/0; = (bq) ,. . , /0 ,, + (for) ... . « ,, + (bq} ,. . . /0 ,,

J I / / I /v Jt' IV * / Yl 1 1 1 /V K \ \ * / Vr 7 / I / y * 1 \ *- J 7( 1 1 1 /V fc' ) I 1 If* *-* *** * I * I /»* * » ^j **/ /l/j* -1' ^*> **/ <*!/»* *-/ -^i"' / I

./? ' L ^*7 ' 7/ ^^^' '/ ^17' (/ J

(A39)

fAs before, all terms that are determined from the flow solution are at time level t+l.

69
APPENDIX B: DATA INPUT INSTRUCTIONS FOR MOC3D

MODFLOW Name File

Transport simulation is activated by including a record in the MODFLOW name file using 

the file type (Ftype) "CONC" to link to the transport name file. The transport name file specifies 

the files to be used when simulating solute transport in conjunction with a simulation of groundwater flow using MODFLOW. The transport name file works in the same way as the MODFLOW 

name file.

MODFLOW Source and Sink Packages

Except for recharge, concentrations associated with fluid sources (C') are read as auxiliary 

parameters in the MODFLOW source package. The source concentration is read from a new 

column appended to the end of each line of the data file describing a fluid sink/source (see 

documentation for revised MODFLOW model; Harbaugh and McDonald, 1996a and 1996b). For 

example, concentrations associated with well nodes should be appended to the line in the WEL 

Package where the well's location and pumping rate are defined. These concentrations will be read 

if the auxiliary parameter "CONCENTRATION" (or "CONC") appears on the first line of the well 

input data file. The concentration in recharge is defined separately, as described in following 

section "Source Concentration in Recharge File."

To simulate solute transport the MODFLOW option enabling storage of cell-by-cell flow 

rates for each fluid source or sink is required in all fluid packages except recharge. The key word 

"CBCALLOCATE" (or "CBC") must appear on the first line of each input data file for a fluid 

package (see Harbaugh and McDonald, 1996a and 1996b).

MOC3D Input Data Files

All input variables are read using free formats, except as specifically indicated. In free 

format, variables are separated by one or more spaces or by a comma and optionally one or more 

spaces. Blank spaces are not read as zeros.

70
MOC3D Transport Name File (CONC)

FOR EACH SIMULATION:

1. Data: FTYPE NUNIT FNAME

The name file consists of records defining the names and units numbers of the files. Each 

"record" consists of a separate line of data. There must be a record for the listing file and for the 

main MOC3D input file.

The listing (or output) file ("CLST") must be the first record. The other files may be in any 

order. Each record can be no more than 79 characters.

FTYPE The file type, which may be one of the following character strings:

CLST MOC3D listing file (separate from the MODFLOW listing file) [required].

MO C Main MOC3D input data file [required].

CR CH Concentrations in recharge [optional].

CNCA Separate output file containing concentration data in ASCII (text-only) format. 

Frequency and format of printing controlled by NPNTCL and ICONFM 

[optional].

CNCB Separate output file containing concentration data in binary format [optional].

VELA Separate output file with velocity data in ASCII format. Frequency and format 

of printing controlled by NPNTVL and IVELFM [optional].

VELB Separate output file with velocity data in binary format [optional].

PR TA Separate output file with particle locations printed in ASCII format. Frequency 

and format of printing controlled by NPNTPL [optional].

PR TB Separate output file with particle locations printed in binary format [optional]. 

OBS Observation wells input file [optional].

DA TA For formatted files such as those required by the OBS package and for array 

data separate from the main MOC3D input data file [optional].

DATA(BINAR Y) For formatted input/output files [optional].

NUNIT The FORTRAN unit number used to read from and write to files. Any legal unit 

number other than 97, 98, and 99 (which are reserved by MODFLOW) can be 

used provided that it is not previously specified in the MODFLOW name file.

FNAME The name of the file.

71
Main MOC3D Package Input (MOC)

Input for the method-of-characteristics (MOC3D) solute-transport package is read from the 

unit specified in the transport name file. The input consists of 18 separate records or data sets, as 

described in detail below. These data are used to specify information about the transport subgrid, 

physical and chemical transport parameters, numerical solution variables, and output formats. 

Output file controls for the MOC3D package are specified in the transport name file, described 

previously.

FOR EACH SIMULATION:

1. Data: HEDMOC A two-line character-string title describing the 

simulation (80 text characters per line).

2. Data: HEDMOC (continued)

3. Data: ISLAYl ISLAY2 ISROWl ISROW2 ISCOLl ISCOL2

ISLAY1 

ISLAY2 

ISROWl 

ISROW2 

ISCOLl 

ISCOL2

Notes:

Transport may be simulated within a subgrid, which is a "window" within the primary 

MODFLOW grid used to simulate flow. Within the subgrid, the row and column spacing must be 

uniform, but thickness can vary from cell to cell and layer to layer. However, as discussed in the 

section reviewing MOC3D assumptions, the range in thickness values (or product of thickness and 

porosity) should be as small as possible.

ISLAYl ISLAY2 ISROWl ISROW2 ISCOLl

Number of first (uppermost) layer for transport.

Last layer for transport.

First row for transport.

Last row for transport.

First column for transport.

Last column for transport.

4. Data: NODISP DECAY DIFFUS

NODISP Flag for no dispersion (set NODISP=1 if no dispersion in problem; this will reduce 

storage allocation).

DECAY First-order decay rate [1/T] (DECAY=0.0 indicates no decay occurs). 

DIFFUS Effective molecular diffusion coefficient [L2/T].

Notes:

The decay rate (A) is related to the half life (?1/2) of a constituent by A = (In 2)/?1/2 . 

The effective molecular diffusion coefficient (Dm) includes the effect of tortuosity.

72
5. Data: NPMAX NPTPND

NPMAX Maximum number of particles available for particle tracking of advective transport 

in MOC3D. If set to zero, the model will calculate NPMAX according to the 

following equation:

NPMAX = 2xNPTPNDxNSROWxNSCOLxNSLAY.

NPTPND Initial number of particles per cell in transport simulation (that is, at t = 0.0). Valid 

options for default geometry of particle placement include 1, 2, 3, or 4 for onedimensional transport simulation; 1,4, 9, or 16 for two-dimensional transport 

simulation; and 1, 8, or 27 for three-dimensional transport simulation. The user 

can also customize initial placement of particles by specifying NPTPND as a 

negative number, in which case the minus sign is recognized as a flag to 

indicate custom placement is desired. In this case, the user must input local 

particle coordinates as described below.

IF NPTPND IS NEGATIVE IN SIGN:

6. Data: PNEWL PNEWR PNEWC

PNEWL Relative position in the layer (z) direction for initial placement of particle within any 

finite-difference cell.

PNEWR Relative position in the row (y) direction for initial placement of particle. 

PNEWC Relative position in the column (x) direction for initial placement of particle.

Notes:

The three new (or initial) particle coordinates are entered sequentially for each of the 

NPTPND particles. Each line contains the three relative local coordinates for the new particles, in 

order of layer, row, and column. There must be NPTPND lines of data, one for each particle. The 

local coordinate system range is from -0.5 to 0.5, and represents the relative distance within the cell 

about the node location at the center of the cell, so that the node is located at 0.0 in each direction.

FOR EACH SIMULATION:

7. Data: CELDIS FZERO INTRPL

CELDIS Maximum fraction of cell dimension that particle may move in one step (typically,

0.5 < CELDIS < 1.0). 

FZERO If the fraction of active cells having no particles exceeds FZERO, the program will

automatically regenerate an initial particle distribution before continuing the simulation

(typically, 0.01 < FZERO < 0.05). 

INTRPL Flag for interpolation scheme used to estimate velocity of particles. The default

(INTRPL=1) will use a linear interpolation routine; if INTRPL=2, a scheme will be

implemented that uses bilinear interpolation in the row and column (/' and i)

directions only (linear interpolation will still be applied in the k, or layer, direction).

(See section "Discussion Choosing appropriate interpolation scheme.")

73
FOR EACH SIMULATION:

8. Data: NPNTCL ICONFM NPNTVL IVELFM NPNTDL IDSPFM NPNTPL

NPNTCL Flag for printing concentration data. If NPNTCL=-2, concentration data will be

printed at the end of every stress period; if NPNTCL=-1, data will be printed at the 

end of every flow time step; if NPNTCL=0, data will be printed at the end of the 

simulation; if NPNTCL=N>0, data will be printed every Nth particle moves, and at 

the end of the simulation. Initial concentrations are always printed.

ICONFM Specification for format of concentration data in main output file (see Table 17 and 

MODFLOW documentation on array-reading utility modules).

NPNTVL Flag for printing velocity data. If NPNTVL=-1, velocity data will be printed at the end 

of every stress period; if NPNTVL=0, data will be printed at the end of the 

simulation; if NPNTVL=N>0, data will be printed every Nth flow time steps, and 

at the end of the simulation.

IVELFM Specification for format of velocity data, if being printed in main output file (see Table 

17).

NPNTDL Flag for printing dispersion equation coefficients that include cell dimension factors 

(see section "Program Segments"). If NPNTDL=-2, coefficients will be printed at 

the end of every stress period; if NPNTDL=-1, coefficients will be printed at the 

end of the simulation; if NPNTDL=0, coefficients will not be printed; if 

NPNTDL=N>0, coefficients will be printed every Nth flow time step.

IDSPFM Specification for format of dispersion equation coefficients (see Table 17).

NPNTPL Flag for printing particle locations in a separate output file (only used if file types 

"PRTA" or "PRTB" appear in the MOC3D name file). If neither "PRTA" or 

"PRTB" is entered in the name file, NPNTPL will be read but ignored (so you must 

always have some value specified here). If either "PRTA" or "PRTB" is entered in 

the name file, initial particle locations will be printed to the separate file first, 

followed by particle data at intervals determined by the value of NPNTPL. If 

NPNTPL=-2, particle data will be printed at the end of every stress period; if 

NPNTPL=-1, data will be printed at the end of every flow time step; if 

NPNTPL=0, data will be printed at the end of the simulation; if NPNTPL=N>0, 

data will be printed every Nth particle moves, and at the end of the simulation.

Table 17. Formats associated with MOC3D print flags

Print flag

0

i

2

3

4

5

6

Format

10G11.4

11G10.3

9G13.6

15F7.1

15F7.2

15F7.3

15F7.4

Print flag

7

8

9

10

11

12

13

Format

20F5.0

20F5.1

20F5.2

20F5.3

20F5.4

10G11.4

10F6.0

Print flag

14

15

16

17

18

Format

10F6.1

10F6.2

10F6.3

10F6.4

10F6.5

74
FOR EACH SIMULATION:

9. Data: CNOFLO Concentration associated with inactive cells of subgrid (used for

output purposes only).

FOR EACH LAYER OF THE TRANSPORT SUBGRID:

10. Data: CINT(NSCOL,NSROW) Initial concentration.

Module: U2DREL*

FOR EACH SIMULATION, ONLY IF TRANSPORT SUBGRID DIMENSIONS ARE 

SMALLER THAN FLOW GRID DIMENSIONS:

11. Data: CINFL (ICINFL) C to be associated with fluid inflow across the

boundary of the subgrid.

Module: UlDREL*

Notes:

The model assumes that the concentration outside of the subgrid is the same within each 

layer, so only one value of CINFL is specified for each layer within and adjacent to the subgrid. That 

is, the size of the array (ICINFL) is determined by the position of the subgrid with respect to the 

entire (primary) MODFLOW grid. If the transport subgrid has the same dimensions as the flow grid, 

this parameter should not be included in the input data set. If the subgrid and flow grid have the 

same number of layers, but the subgrid has fewer rows or fewer columns, ICINFL=NSLAY. Values 

are also required if there is a flow layer above the subgrid and/or below the subgrid. The order of 

input is: C" for first (uppermost) transport layer (if required); C' for each successive (deeper) 

transport layer (if required); C" for layer above subgrid (if required); and C" for layer below 

subgrid (if required).

FOR EACH SIMULATION

12 . Data: NZONES Number of zone codes among fixed-head nodes in transport subgrid. 

IF NZONES > 0:

Data: IZONE ZONCON

I ZONE Value identifying a particular zone.

ZONCON Source concentration associated with nodes in the zone defined by IZONE above.

Notes:

Zones are defined within the IBOUND array in the BAS Package of MODFLOW by 

specifying unique negative values for fixed-head nodes to be associated with separate fluid source 

concentrations. Each zone is defined by a unique value of IZONE and a concentration associated 

with it (ZONCON). There must be NZONES lines of data, one for each zone. Note that values of 

IZONE in this list must be negative for consistency with the definitions of fixed-head nodes in the 

IBOUND array in the BAS Package. If a negative value of IBOUND is defined in the BAS package 

but is not assigned a concentration value here, MOC3D will assume that the source concentrations 

associated with those nodes equal 0.0.

Module is a standard MODFLOW input/output module.

75
FOR EACH LAYER OF THE TRANSPORT SUBGRID:

13. Data: IGENPT (NSCOL,NSROW) Flag to treat fluid sources and sinks as

either "strong" or "weak."

Module: U2DINT*

Notes:

Where fluid source is "strong," new particles are added to replace old particles as they are 

advected out of that cell. Where a fluid sink is "strong," particles are removed after they enter that 

cell and their effect accounted for. Where sources or sinks are weak, particles are neither added nor 

removed, and the source/sink effects are incorporated directly into appropriate changes in particle 

positions and concentrations. If IGENPT=0, the node will be considered a weak source or sink; if 

IGENPT=1, it will be a strong source or sink. See section on "Special Problems" and discussion by 

Konikow and Bredehoeft (1978).

IF NODISP * I (If dispersion is included in simulation):

14. Data: ALONG (NSLAY) Longitudinal dispersivity.

Module: UlDREL*

15 . Data: ATRANH (NSLAY) Horizontal transverse dispersivity.

Module: UlDREL*

16 . Data: ATRANV (NSLAY) Vertical transverse dispersivity.

Module: UlDREL*

FOR EACH SIMULATION:

17 . Data: RF (NSLAY) Retardation factor (RF=1 indicates no retardation).

Module: UlDREL*

Notes:

If RF=0.0 in input, the code automatically resets it as RF=1.0 to indicate no retardation.

FOR EACH LAYER OF TRANSPORT SUBGRID:

18a. Data: THCK (NSCOL,NSROW) Cell thickness.

Module: U2DREL*

18b. Data: FOR (NSCOL, NSROW) Cell porosity.

Module: U2DREL*

Notes:

The thickness and porosity are input as separate arrays for each layer of the transport 

subgrid. The sequence used in data set 18 is to first define the thickness of the first layer of the 

transport subgrid, and then define the porosity of that same layer. Next, that sequence is repeated for 

all succeeding layers. The product of thickness and porosity should not be allowed to vary greatly 

among cells in the transport subgrid.

Module is a standard MODFLOW input/output module.

76
Source Concentration in Recharge File (CRCH)

Concentrations in recharge, if the recharge package is used, are read from a separate unit 

specified in the MOC3D name file. This is defined using the file type (Ftype) "CRCH."

FOR EACH STRESS PERIOD, IF RECHARGE PACKAGE USED:

1. Data: INCRCH Flag to reuse or read new recharge concentrations.

Notes:

Read new recharge concentrations if INCRCH > 0. Reuse recharge concentrations from the 

last stress period if INCRCH < 0.

2 . Dat a: CRECH (NSCOL, NSROW) Source concentration associated with fluid

entering the aquifer in recharge.

Module: U2DREL*

Observation Well File (OBS)

Nodes of the transport subgrid can be designated as "observation wells." At each such 

node, the time, head, and concentration after each move increment will be written to a separate 

output file to facilitate graphical postprocessing of the calculated data. The input file for specifying 

observation wells is read if the file type (Ftype) "OBS" is included in the MOC3D name file.

FOR EACH SIMULATION, IF OBS PACKAGE USED:

1. Data: NUMOBS IOBSFL

NUMOBS Number of observation wells.

IOBSFL If IOBSFL = 0, well data are saved in NUMOBS separate files. If IOBSFL>0, all 

observation well data will be written to one file, and the file name and unit 

number used for this file will be that of the first observation well in the list.

FOR EACH OBSERVATION WELL:

2. Data: LAYER ROW COLUMN UNIT

LAYER Layer of observation well node.

ROW Row of observation well node.

COLUMN Column of observation well node.

UNIT Unit number for output file.

Notes:

If NUMOBS>1 and IOBSFL = 0, you must specify a unique unit number for each observation well and match those unit numbers to DATA file types and file names in the MOC3D name file. 

If IOBSFL>0, you must specify a unique unit number for the first observation well and match that 

unit number to a DATA file type and file name in the MOC3D name file.

Module is a standard MODFLOW input/output module.

77
APPENDIX C: ANNOTATED EXAMPLE INPUT DATA SET FOR 

SAMPLE PROBLEM

This example input data set is the one used to generate the solution for the base case in the 

one-dimensional steady-state flow problem. Parameter values are indicated in Table 11 and 

selected results are shown in fig. 18. Several of the following data files (mod/low.nam, 

bas95.dat, bcfll.dat, and sipl9.dat) are those required for MODFLOW-96, and their formats are 

described by Harbaugh and McDonald (1996a).

In several of the data files shown below, the right side of some data lines includes a semicolon followed by text that describes the parameters for which values are given. These comments 

(including the semi-colon) are not read by the program because in free format the code will only 

read the proper number of variables and ignore any subsequent information on that line. This style 

of commenting data files is optional, but users may find it helpful when viewing the content of data 

files.

Following (enclosed in a border) are the contents of the MODFLOW name file for the 

sample problem; explanations are noted outside of border:

Filename: mod flow, nam

list

bas

bcf

sip

cone

T

i

16

95

11

19

33

T

2

flow. out

bas95.dat

bcfll.dat

sipl9.dat

moc .nam

t

3

4 Designates main output file for MODFLOW

4- Basic input data for MODFLOW

4 Block-centered flow package

4 Input for numerical solution of flow equation

4 Transport name file (turns transport "on")

1 Ftype (that is, the type of file)

2 Unit number

* File name (name chosen to reflect contents of file)

78
Following (enclosed in a border) are the contents of the basic package input data file for the 

MODFLOW simulation of the sample problem; explanations are noted outside of border:

File name: bas95. dat

Finite: Compare

NLAY

FREE

-1 1

1 1

1 1

1 1

1 1

0

12.1

120.0

end

1

0

95

1 1

1 1

1 1

1 1

1 1

00

95

to Wexler program and MOC3D

NROW

1

1

NCOL

122

NPER

1

BAS Input

ITMUNI

1

; IAPART, ISTRT

1(2513)

1

1

1

1

1

1 1

1 1

1 1

1 1

1 1

111

111

111

111

111

111

111

111

111

111

1 1

1 1

1 1

1 1

1 1

1

1

1

1

1

1.0(122F5.0)

1 1. ; PERLEN, NSTP,

1 1

1 1

1 1

1 1

1 1

;

1

3

1 1

1 1

1 1

1 1

1 1

HNOFLO

HEAD

; IBOUND

11111

11111

11111

11111

1 -2

TSMULT

<- 1

<- 1

^- 2

<- 3

<- 4

<r- 5

<- 5

<- 5

<- 5

<- 5

<- 5

^6

<- 7

<- 7

^- 8

<- 9

Two header lines of comments. For convenience and clarity, the second line is used to label names of 

parameters on subsequent line of file.

Flow grid dimensions, number of periods, and time units.

Options line (new in MODFLOW-96)

Flags for buffer array and drawdown calculations.

IBOUND identifiers (first line) and array

Head value assigned to inactive cells

Initial head information

MODFLOW time-step information

Final comment line

Following (enclosed in a border) are the contents of the block-centered flow package input data 

file; explanations are noted outside of border:

File name: bcf 11 . dat

1 0 0.0

0

0

0

0

0

0 0.0 0 0

1.0

0.1

0.1

0.01

; ISS, flags

; LAYCON

; TRPY

; DELR

; DELC

; TRAN

BCF Input <- 1

<-3

1 Flag for steady-state flow, flag for cell-by-cell flow terms, five flags related to wetting

2 Layer type

3 Anisotropy factor

4 Grid spacing information

^ Transmissivity data

79
Following (enclosed in a border) are the contents of the strongly implicit procedure package input 

data file; explanations are noted outside of border:

File name: sipl 9. da t

500 5

1. 0.0000001

; MXITER,NPARM SIP Input 

0 0.001 ; ACCL,ERR,IPCALC,WSEED

<- 1

Maximum iterations, number of iteration parameters 

Acceleration parameter, head change criterion, flag for seed, seed

Following (enclosed in a border) are the contents of the MOC3D name file for the sample problem; 

explanations are noted outside of border:

File name: moc. nam

cist

moc

obs

data

cnca

cncb

T

i

97

96

44

45

22

23

T

2

moc . out

moc96 . dat

obs 4 4 .dat

obs . out

cone . txt

cone .bin

t

3

1 Ftype

2 Unit number

3 File name

< Designates main output file for MOC3D 

< Main input data file for MOC3D 

< Input data file for observation wells 

< Output file for observation well data 

< Separate output file for cone, data (ASCII) 

< Separate output file for cone, data (binary)

80
Following (enclosed in a border) are the contents of the main input data file for the MOC3D 

simulation for the sample problem; selected explanations are noted outside of border:

File name: moc9 6. da t

One-dimensional , Steady Flow, No

ISLAY1 ISLAY2 ISROWl

1

0

2000

0.5 0

000-1000

0.0

0

0

2

-1

-2

0

0

0

0

0

0

0

1 1

0.0 0.0 ;

3

.05 2 ;

Decay, Low Dispersion: BASE CASE MOC3D Input

ISROW2 ISCOL1 ISCOL2

1 2 121

NODISP, DECAY, DIFFUS

NPMAX, NPTPND

CELDIS, FZERO, INTRPL

; NPNTCL, ICONFM, NPNTVL, IVELFM, NPNTDL, IDSPFM, NPRTPL

; CNOFLO

0.0 (122F3.0)

1.

1.0

0.0

0

0.1

0.1

0.1

1.0

1.0

0.1

; initial concentration

; C ' inflow

; NZONES to follow

; I ZONE, ZONCON

; I ZONE, ZONCON

; IGENPT

; longitudinal disp.

; transverse disp. horiz.

; transverse disp. vert.

; retardation factor

; thickness

; porosity

<- 1

<- 1

<- 2

<- 3

^ _ A

A __ A

<- 5

< 6

<- 7

<- 7

<- 7

<- 8

Two header lines of comments. For convenience and clarity, the second line is used to label names of 

parameters on subsequent line of file.

Indices for transport subgrid

Flag for no dispersion, decay rate, diffusion coefficient 

Particle information for advective transport 

Print flags 

" Value of concentration associated with inactive cells

' Concentrations associated with fixed-head nodes (fixed head nodes are defined in the IBOUND array in the 

MODFLOWBAS package)

8 Flag for "strong" sources or sinks

Following (enclosed in a border) are the contents of the observation well input data file for the 

sample problem; explanations are noted outside of border:

File name: obs 4 4 . da t

3 1

112

1 1 42

1 1 112

45

;NUMOBS IOBSFL Observation well data

;layer, row, column, unit number

;layer, row, column

;layer, row, column

<- 1

Number of observation wells, flag to print to one file or separate files

Node location and unit number for output file (linked to the Ftype DATA in MOC3D name file)

81
APPENDIX D: SELECTED OUTPUT FOR SAMPLE PROBLEM

This example output was generated from the input data sets listed in Appendix C for the 

base case of the one-dimensional steady-state flow problem. The line spacing and font sizes of the 

output files have been modified in places to enhance the clarity of reproduction in this report. 

Some repetitive lines of output have been deleted where indicated by an ellipsis (...).

Some brief annotations were added in a few places in this sample output listing to help the 

reader understand the purpose of various sections of output. These annotations are written in bold 

italics to clarify that they are not part of the output file.

Following are the contents of the MOC3D main output file for the sample problem.

FILE INFORMATION

U.S. GEOLOGICAL SURVEY

METHOD-OF-CHARACTERISTICS SOLUTE TRANSPORT MODEL 

MOC3D (Version 1.0) 11/08/96

HOC BASIC INPUT READ FROM UNIT 

LISTING FILE: moc.out UNIT 97

OPENING moc96.dat

FILE TYPE: HOC UNIT 96

OPENING obs44.dat

FILE TYPE: OBS UNIT 44

OPENING obs.out

FILE TYPE: DATA UNIT 45

OPENING conc.txt

FILE TYPE: CNCA UNIT 22

OPENING cone.bin

FILE TYPE: CNCB UNIT 23

MOC BASIC INPUT READ FROM UNIT 96

2 TITLE LINES:

One-dimensional, Steady Flow, No Decay, Low Dispersion: BASE CASE MOC3D Input 

ISLAY1 ISLAY2 ISROWl ISROW2 ISCOLl ISCOL2

PROBLEM DESCRIPTORS, INCLUDING GRID CHARACTERISTICS AND PARTICLE INFORMATION:

MAPPING OF SOLUTE TRANSPORT SUBGRID IN FLOW GRID:

FIRST LAYER FOR SOLUTE TRANSPORT = 1 LAST LAYER FOR SOLUTE TRANSPORT = 1 

FIRST ROW FOR SOLUTE TRANSPORT = 1 LAST ROW FOR SOLUTE TRANSPORT = 1 

FIRST COLUMN FOR SOLUTE TRANSPORT= 2 LAST COLUMN FOR SOLUTE TRANSPORT = 121

UNIFORM DELCOL AND DELROW IN SUBGRID FOR SOLUTE TRANSPORT

NO. OF LAYERS = 1 NO. OF ROWS =

NO SOLUTE DECAY

NO MOLECULAR DIFFUSION

MAXIMUM NUMBER OF PARTICLES (NPMAX) =

NO. OF COLUMNS = 120

2000

82
14006 ELEMENTS IN X ARRAY ARE USED BY MOC 

12 ELEMENTS IN X ARRAY ARE USED BY OBS

NUMBER OF PARTICLES INITIALLY IN EACH ACTIVE CELL (NPTPND) = 3 

PARTICLE MAP ("o" indicates particle location; shown as

fractions of cell distances relative to node location)

o o o

-1/3 0 1/3

INITIAL RELATIVE PARTICLE COORDINATES

1 0.00000 0.00000 -0.33333

2 0.00000 0.00000 0.00000

3 0.00000 0.00000 0.33333

CELDIS= 

FZERO =

0.500

0.050

INTRPL= 2; BILINEAR INTERPOLATION SCHEME

NPNTCL= 0: CONCENTRATIONS WILL BE WRITTEN AT THE END OF THE SIMULATION 

MODFLOW FORMAT SPECIFIER FOR CONCENTRATION DATA: ICONFM= 0

NPNTVL= 0: VELOCITIES WILL BE WRITTEN AT THE END OF THE SIMULATION 

MODFLOW FORMAT SPECIFIER FOR VELOCITY DATA: IVELFM= -1

NPNTDL= 0: DISP. COEFFICIENTS WILL NOT BE WRITTEN 

NPNTPL= 0: PARTICLE LOCATIONS WILL NOT BE WRITTEN

CONCENTRATION WILL BE SET TO O.OOOOOE+00 AT ALL NO-FLOW NODES (IBOUND=0). 

INITIAL CONCENTRATION = 0.OOOOOOOE+00 FOR LAYER 1

OUTPUT 

CONTROL

VALUES OF C 1 REQUIRED FOR SUBGRID BOUNDARY ARRAY = 1 

ONE FOR EACH LAYER IN TRANSPORT SUBGRID

ORDER OF C 1 VALUES: FIRST LAYER IN SUBGRID, EACH SUBSEQUENT LAYER, 

LAYER ABOVE SUBGRID, LAYER BELOW SUBGRID:

SUBGRID BOUNDARY ARRAY = 1.000000

NUMBER OF ZONES FOR CONCENTRATIONS AT FIXED HEAD CELLS =

ZONE FLAG = -1 

ZONE FLAG = -2

INFLOW CONCENTRATION = l.OOOOE+00 

INFLOW CONCENTRATION = O.OOOOE+00

SINK-SOURCE FLAG = 0 FOR LAYER

LONGITUDNL. DISPERSIVITY = 0.1000000

HORIZ. TRANSVERSE DISP. = 0.1000000

VERT. TRANSVERSE DISP. = 0.1000000

RETARDATION FACTOR = 1.000000

INITIAL THICKNESS = 1.000000 FOR LAYER

INITIAL POROSITY = 0.1000000 FOR LAYER

INITIAL AND 

BOUNDARY 

CONDITIONS 

FOR SOLUTE

83
COORDINATES FOR 3 OBSERVATION WELLS:

WELL # LAYER ROW COLUMN UNIT

1 1 1 2 45

2 1 1 42 45

3 1 1 112 45 

ALL OBSERVATION WELL DATA WILL BE WRITTEN ON UNIT 45

CONCENTRATION DATA WILL BE SAVED ON UNIT 22 IN ASCII FORMAT 

CONCENTRATION DATA WILL BE SAVED ON UNIT 23 IN BINARY FORMAT

TOTAL NUMBER OF PARTICLES GENERATED (GENPT) = 360 

TOTAL NUMBER OF ACTIVE NODES (NACTIV) = 120

MAX. NUMBER OF CELLS THAT CAN BE VOID OF PARTICLES (NZCRIT) = 6 

(IF NZCRIT EXCEEDED, PARTICLES ARE REGENERATED)

CALCULATED VELOCITIES (INCLUDING EFFECTS OF RETARDATION, IF PRESENT):

EFFECTIVE MEAN SOLUTE VELOCITIES IN COLUMN DIRECTION

AT NODES

1

VELOCITY (COL) IN LAYER 1 AT END OF TIME STEP 1 IN STRESS PERIOD 1

123456789... 

1 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02 ...

Ill 112 113 114 115 116 117 118 119 120 

1 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02 9.917E-02

EFFECTIVE MEAN SOLUTE VELOCITIES IN ROW DIRECTION

AT NODES

1

VELOCITY (ROW) IN LAYER 1 AT END OF TIME STEP 1 IN STRESS PERIOD 1

123456789... 

1 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 ...

Ill 112 113 114 115 116 117 118 119 120 

1 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00

EFFECTIVE MEAN SOLUTE VELOCITIES IN LAYER DIRECTION

AT NODES

1

VELOCITY (LAYER) IN LAYER 1 AT END OF TIME STEP 1 IN STRESS PERIOD 1

123456789... 

1 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 ...

Ill 112 113 114 115 116 117 118 119 120 

1 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00

84
STABILITY CRITERIA M.O.C.

MAXIMUM FLUID VELOCITIES: C-VEL = 9.92E-02 R-VEL = l.OOE-20 L-VEL = l.OOE-18 

MINIMUM TIME TO TRAVEL THCK = l.OOE+18

TIMV = 5.04E-01 NTIMV = 239

MAX. C-VEL. IS CONSTRAINT AND OCCURS BETWEEN NODES ( 2, 1, 1) AND ( 1, 1, 1) 

TIMD = 5.04E-01 NTIMD = 239

THERE ARE NO FLUID SOURCES IN THE TRANSPORT SUBGRID

NUMBER OF MOVES FOR ALL STABILITY CRITERIA: 

CELDIS DISPERSION INJECTION 

239 239 1

CELDIS IS LIMITING 

DISPERSION IS LIMITING

NO. OF PARTICLE MOVES REQUIRED TO COMPLETE THIS TIME STEP = 239 

MOVE TIME STEP (TIMV)= 5.020920634270E-01

(NUMERICAL SOLUTION TO TRANSPORT EQUATION STARTS AT THIS POINT)

NP 

NP 

NP 

NP 

NP 

NP

NP 

NP 

NP 

NP 

NP 

NP

360 AT START 

360 AT START 

360 AT START 

360 AT START 

360 AT START 

360 AT START

360 AT START 

360 AT START 

360 AT START 

360 AT START 

360 AT START 

360 AT START

OF MOVE 

OF MOVE 

OF MOVE 

OF MOVE 

OF MOVE 

OF MOVE

OF MOVE 

OF MOVE 

OF MOVE 

OF MOVE 

OF MOVE 

OF MOVE

IMOV 

IMOV 

IMOV 

IMOV 

IMOV 

IMOV

IMOV 

IMOV 

IMOV 

IMOV 

IMOV 

IMOV

234

235

236

237

238

239

ONE LINE PRINTED 

FOR EACH MOVE TO 

TRACK PROGRESS 

AND NUMBER OF 

ACTIVE PARTICLES

85
SOLUTE BUDGET AND MASS BALANCE FOR TRANSPORT SUBGRID

VALUES CALCULATED AT END OF: 

STRESS PERIOD 1 

FLOW TIME STEP 1 

TRANSPORT TIME INCREMENT 239

ELAPSED TIME = 1.2000E+02

OUT OF 

OUT OF 

OUT OF

1

1

239

CHEMICAL MASS IN STORAGE:

INITIAL: MASS DISSOLVED = O.OOOOE+00 

PRESENT: MASS DISSOLVED = 1.1341E-01

CHANGE IN MASS STORED = -1.1341E-01

MASS SORBED 

MASS SORBED

O.OOOOE+00 

O.OOOOE+00

CUMULATIVE SOLUTE MASS (L**3) (M/VOL)

IN:

DECAY = O.OOOOE+00

CONSTANT HEAD = O.OOOOE+00

SUBGRID BOUNDARY = 1.190IE-01

RECHARGE = O.OOOOE+00

WELLS = O.OOOOE+00

RIVERS = O.OOOOE+00

DRAINS = O.OOOOE+00

GENL. HEAD-DEP. BDYS. = O.OOOOE+00

EVAPOTRANSPIRATION = O.OOOOE+00

TOTAL IN = 1.1901E-01

OUT:

DECAY = O.OOOOE+00

CONSTANT HEAD = O.OOOOE+00

SUBGRID BOUNDARY = -5.6659E-03

RECHARGE = O.OOOOE+00

WELLS = O.OOOOE+00

RIVERS = O.OOOOE+00

DRAINS = O.OOOOE+00

GENL. HEAD-DEP. BDYS. = O.OOOOE+00

EVAPOTRANSPIRATION = O.OOOOE+00

TOTAL OUT = -5.6659E-03

SOURCE-TERM DECAY = O.OOOOE+00

ITEMIZED 

BUDGETS FOR 

SOLUTE FLUXES

RESIDUAL 

PERCENT DISCREPANCY

-6.7927E-05

-0.5708E-01 RELATIVE TO MASS FLUX IN

86
Following (enclosed in a border) are the abridged contents of the observation well output file for 

the sample problem. This output file was generated using the option to write all observation well 

data to a single file (IOBSFL = 1).

File name: obs. ou t

"OBSERVATION WELL DATA"

"TIME, THEN

"

5

1

1

2

2

3

1

1

1

1

TIME:

.0209E-01

.0042E+00

.5063E+00

.0084E+00

.5105E+00

.0126E+00

.1849E+02

.1900E+02

.1950E+02

.2000E+02

HEAD AND CONC.

1

1

1

1

1

1

1

1

1

1

H & C AT

.190E+01

.190E+01

.190E+01

.190E+01

.190E+01

.190E+01

.190E+01

.190E+01

.190E+01

.190E+01

1

2

6

5

7

7

8

1

1

1

1

FOR EACH

, 1, 2

.503E-01

.539E-01

.994E-01

.914E-01

.747E-01

.578E-01

.OOOE+00

.OOOE+00

.OOOE+00

.OOOE+00

OBS. WELL AT NODE (K,I,J)"

H & C AT

7.934E+00

7.934E+00

7.934E+00

7.934E+00

7.934E+00

7.934E+00

7.934E+00

7.934E+00

7.934E+00

7.934E+00

1,

0.

0.

0.

0.

0.

0.

1.

1.

1.

1.

1, 42

OOOE+00

OOOE+00

OOOE+00

OOOE+00

OOOE+00

OOOE+00

OOOE+00

OOOE+00

OOOE+00

OOOE+00

H

9.

9.

9.

9.

9.

9.

9.

9.

9.

9.

& C AT

917E-01

917E-01

917E-01

917E-01

917E-01

917E-01

917E-01

917E-01

917E-01

917E-01

1, 1,112 "

0. OOOE+00

0. OOOE+00

0. OOOE+00

0. OOOE+00

0. OOOE+00

0. OOOE+00

6.835E-01

6.900E-01

7.045E-01

7.125E-01

Following (enclosed in a border) are the partial contents of the separate output file for concentration 

in ASCII format. Initial concentrations are abridged; complete set of final concentrations are 

shown.

File name: cone. txt

CONCENTRATIONS AT NODES IN SUBGRID. IMOV= 0, NSTP= 0, NPER= 1, SUMTCH=O.OOOOE+00

SUBGRID LAYER 1 

O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00

O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00 O.OOOOE+00

CONCENTRATIONS AT 

SUBGRID LAYER

NODES IN SUBGRID. IMOV= 239, NSTP= 1, NPER= 1, SUMTCH=1.2000E+02

.OOOOE+00

.OOOOE+00

.OOOOE+00

.OOOOE+00

.OOOOE+00

.9999E-01

.9990E-01

9.9904E-01

9.9330E-01

9.6710E-01 9

, OOOOE+00 

.OOOOE+00 

.OOOOE+00 

.OOOOE+00 

.OOOOE+00 

.9999E-01 

.9988E-01 

.9882E-01 

.9201E-01 

.6207E-01

8.8536E-01 8.7219E-01 

7.1254E-01 6.8920E-01

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

9.9999E-01

9.9984E-01

9.9854E-01

9.9050E-01

9.5647E-01

8.5826E-01

6.6520E-01

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

9998E-01

9.9980E-01

9.9821E-01

9.8877E-01

9.5035E-01

8.4369E-01

6.4163E-01

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

9998E-01

9975E-01

9.9782E-01

9.8678E-01

9.4358E-01

8.2811E-01

6.1868E-01

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

9997E-01

9.9968E-01

9.9735E-01

9.8449E-01

9.3594E-01

8.1104E-01

5.9552E-01

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

9.9997E-01

9.9960E-01

9.9679E-01

9.8182E-01

9.2738E-01 

7.9266E-01 

5.7147E-01

.OOOOE+00

.OOOOE+00

.OOOOE+00

.OOOOE+00

.OOOOE+00

9.9996E-01

9.9950E-01

9.9611E-01

9.7876E-01

9.1810E-01

7.7364E-01

5.4711E-01

.OOOOE+00

.OOOOE+00

.OOOOE+00

.OOOOE+00

.OOOOE+00

9.9994E-01

9.9938E-01

9.9532E-01

9.7534E-01

9.0818E-01

7.5433E-01

5.2527E-01

OOOOE+00

OOOOE+00

OOOOE+00

OOOOE+00

9999E-01

9993E-01

9922E-01

9439E-01

9.7149E-01

8.9737E-01

7.3423E-01

5.1194E-01

87