PB85-214443
PLUME 3D: THREE-DIMENSIONAL PLUMES IN UNIFORM GROUND
WATER FLOW
Oklahoma State University
Stillwater, OK
Jun 85
U.S. DEPARTMENT OF COMMERCE
National Technical Information Service
-------
EPA/600/2-85/067
June 1985
PLUME3D :
THREE-DIMENSIONAL PLUMES IN UNIFORM GROUND WATER FLOW
Jan Wagner and Stephanie A. Watts
School of Chemical Engineering
Oklahoma State University
' Stillwater, Oklahoma 74078
Douglas C. Kent
Department of Geology
Oklahoma State University
Stillwater, Oklahoma 74078
CR811142
Project Officer
Carl G. Enfield
Robert S. Kerr Environmental Research Laboratory
U. S. Environmental Protection Agency
Ada, Oklahoma 74820
ROBERT S. KERR ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ADA, OK 74820
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1. REPORT NO. 2. 3. RECIF
EPA/ 600/2-85/067 52? r
4 Lr V* —
4. TITLE AND SUBTITLE 5. REPO
PLUME 3D: THREE-DIMENSIONAL PLUMES IN UNIFORM J'
GROUM) WATER FLOW. 6. PERF
7. AUTHOR(S) 8. PERF
Jan Wagner, Stephanie Watts, Douglas Kent
9. PERFORMING ORGANIZATION NAME AND ADDRESS 10, PRO
Oklahoma State University
Stillwatrr OK 7/1070
12. SPONSORING AGENCY NAME AND ADDRESS 13 TYP
Robert S. Kerr Environmental Research Laboratory Fine
Office of Research and Development 14. SPO
U.S. Environmental Protection Agency
Ada, OK .74820 EPA/
"lENT'S ACCESSION NO.
5 2' k •'' •'; > /AS
RT DATE
jne 1985
ORMING ORGANIZATION CODE
ORMING ORGANIZATION REPORT NO.
GRAM ELEMENT NO.
ABRD1A
«»»xxK3««woxx«. Uoop. Agr.
CR-811142
E OF REPORT AND PERIOD COVERED
il 9/83 - 2/85
MSORING AGENCY CODE
'600/15
15. SUPPLEMENTARY NOTES
Carl G. Enfield, Project Officer
16. ABSTRACT
A _ closed-form analytical solution for three-dimensional plumes was incorporate
. in an interactive computer program. The assumption of an infinite aquifer depth
and uniform source mass rate and source location was overcome by using the
principal of superposition in space and time. The source code was written in a
InDTDrt °f FORTRAN 77 and can be compiled with FORTRAN IV, FORTRAN 66 as well as
hURTRAH 77. As a result, the code is nearly independent of hardware and operatina
system.
17. KEY WORDS AND DOCUMENT' ANALYSIS
a. DESCRIPTORS b. IDENTIFIERS/OPEN ENDE
18. DISTRIBUTION STATEMENT 19. SECURITY CLASS (This
UNCLASSIFIED
RELEASE TO PUBLIC 20. SECURITY CLASS (This
UNCLASSIFIED
D TERMS c. COSATI Field/Group
Report! 21 . NO. OF PAGES
82
yagel 22. PRICE
EPA Form 2220-1 (R«v. 4-77) . PREVIOUS EDITION is OBSOLETE ±
-------
DISCLAIMER
The information in this document has been funded wholly or in part
by the United States Environmental Protection Agency under assistance
agreement number CR-811142 to Oklahoma State University. It has been
subject to the Agency's peer and administrative review, and it has been
approved for publication as an EPA document.
-------
FOREWORD
EPA is charged by Congress to protect the Nation's land, air, and water
systems. Under a mandate of national environmental laws focused on air and
water quality,, solid waste management and the control of toxic substances,
pesticides, noise, and radiation, the Agency strives to formulate and imple-
ment actions which lead to a compatible balance between human activities and
the ability of natural systems to support and nurture life.
The Robert S. Kerr Environmental Research Laboratory is the Agency's
center of expertise for investigation of the soil and subsurface.environment.
Personnel at the Laboratory are responsible for management of research pro-
grams to: (a) determine the fate, transport and transformation rates of
pollutants in the soil, the unsaturated zone .and the saturated zones of the
subsurface environment; (b) define the processes to be used in characterizing
the soil and subsurface environment as a receptor of pollutants; (c) develop
techniques for predicting the effect of pollutants on ground water, soil and
indigenous organisms; and (d) define and demonstrate the applicability and
limitations of using natural processes, indigenous to the soil and subsurface
environment, for the protection of this resource.
This project was initiated to develop an interactive computer model which.
could be utilized to predict'toxic chemical fate in homogeneous aquifers. The
model should be useful in making comparisons, between chemicals, for idealized
homogeneous aquifers. This model is not intended for addressing site specific
problems where there is significant heterogeneity in the aquifer.
Clinton W. Hall /
Director
Robert S. Kerr Environmental
Research Laboratory
-------
ABSTRACT
A closed-form analytical solution for two dimensional plumes was
incorporated in an interactive computer program. The assumption of an infinite
aquifer depth and uniform source mass rate and source location was overcome
by using the principal of superposition in space and time. The source code
was written in a subset of FORTRAN 77 and can be compiled with- FORTRAN IV,
FORTRAN 66 as well as FORTRAN 77. As a result, the code is nearly independent
of hardware and operating system.
-------
TABLE OF CONTENTS
INTRODUCTION 1
SECTION I - MATHEMATICAL DEVELOPMENT 2
Analytical Solution . 3
Assumptions and Limitations •. 7
Superposition 8
SECTION II - COMPUTER PROGRAM ' 16
Basic Input Data 16
Edit Commands 20
SECTION III - APPLICATIONS .23
Site Location .23
Hydrogeology 23
Plating-Waste Contamination 24
Analytical Model '. -. 25
REFERENCES 32
APPENDIX'A - Input Data Dialog and Results for Example Problem A-l
APPENDIX.B - Listing'of Source Code for PLUME3D R-l
APPENDIX C - Solute Transport in Uniform Ground-Uater Flow
Mathematical Development C-l'
-------
LIST OF TABLES
Table 1 - Edit Commands 22
Table 2 - Input Data Required for the Analytical Three-Dimensional
PI ume Model < 26
Table 3 - Typical Values for Aquifer Properties 28
Table 4 - Model Input for Example Problem 1 ....30
-------
LIST OF FIGURES
Figure 1 - Decomposition of a variable rate using superposition
in time 10
Figure-2 - Use of image sources to account for aquifers of
finite depth 13
Figure 3 - Superposition in space to account for barrier boundaries.' 14
VI 1
-------
INTRODUCTION
This document describes a mathematical model and the associated computer
program which can be used to estimate concentration distributions in'a
leachate plume which emanates from a point source. The model includes both
linear adsorption and first-order reactions.
The use of the computer program is fairly simple, but represents only one
tool which can aid in the analysis and understanding of ground-water
contamination problems. The user must select the appropriate tools for the
problem at. hand, based.on a sound, understanding of the principles of ground-
water hydrology, the physical problem, and the assumptions and limitations of
the mathematical model!
-------
SECTION I
MATHEMATICAL DEVELOPMENT
The differential equation describing the conservation of mass of a
component in a saturated, homogeneous aquifer with uniform, steady flow in the
x-direction can be written as
where
o
C = component mass per unit volume of fluid phase M/L
D = dispersion coefficient in x-direction . L^/t
A
D = dispersion coefficient in y-direction L^/t
p
0 = dispersion coefficient in z-direction L /t
R^ = retardation coefficient
V = average interstitial velocity in x-direction L/t
x,y,z = rectangular coordinates L
X = first-order decay constant 1/t
The retardation coefficient accounts for partitioning of the component between
the fluid and solid phases using a linear adsorption isotherm, and is defined
as
-------
where
Pg = bulk density of the aquifer M/L
0 = effective porosity
M/M
KH = distribution constant for a linear adsorption isotherm Q
• . M/LJ
Analytical Solution
A closed-form analytical solution to Equation 1 can be obtained by making
a change of variables. Let
t=.t/Rd (3)
and
X = x - - V*T (4)
Then Equation 1 is transformed to
222
14+D -^4+D 14-R.AC' (5)
x 7? y ay2 z sz2 d
with boundary conditions
C(X,y,z,o) = 0 (6a)
C(X,y,±-,T) = 0 (6b)
C(X,±-,z,T) = 0 (6c)
-------
C(±-,y,z,T) = 0 (6d)
Equation 5 is a special form of the analogous equations in heat conduction,
for which solutions are given by Carslaw and Jaeger (1959). The solution for
an instantaneous point source of strength M0 at the origin is
M / 2 2 2
C =
i -4077- 4-7- 40
r3 tJ D D D x y Z
x y z
(7)
In terns of the untransformed variables,
MQ - / (x - V t/Rd) 2
'' ° ^ .3 t3 0 0 D" 6Xt 4°* "Rd ' *V^
x y z x
z
4Dz t/Rd-
- U ' •(8)
The solution for a continuous point source is obtained by integrating
Equation 8 with respect to time and letting C0Q = dM0/dt, or
/ ~K
I I v _ U 1- /t^
- (x - V*t/R)2
C =
r
/ C0
J
'c Rfl 7-^ / ° ^ t- 4D t/R
89 7 7T3 D D D J X X d
x y z o
2
J^ ± ^i dt
4D t/Rd " 4Dz t/R
-------
Equation 9 can be rearranged slightly to
*
1 V x
exp -
c =
8e / 3 n n n ^
TT D D D o
x y z
A
Q t"3/2 exp
Rd R U2t
4DxRd
dt (10)
where
and
(11)
4D R
vl/2
' - U = V
(12)
Turner (1972) gives the solution to Equation 10 as
Cc =
/ * \
fl V .x
2 in
\ x /
/ D D
y 2
r' erf
x
l/U2t
2l°xRd
_ RdR
2\ Dxt
2\l/2.
,2\l/2-
(13)
-------
which can be simplified somewhat to yield
C =
8ir9R
Rd R + Ut
exp [^) erfc.[£
1 x/ r / R, D t
.
d x
-*M-^
A
erfc
Rd R - Ut
(14a)
for x > 0.
For x < 0, only the first term, exp (V x/Dx), has the sign of x altered,
because the value of the integral in Equation 10 is the same for both positive
and negative values of x. Therefore
cr _ x_
C SuSR / D DZ
e*p,4r
R + Utx
erfc
exp I- ^ rr-1 erfc
- Ut
v
(14b)
for x < 0.
-------
The steady-state solution for a continuous point source is found from
Equation. 8 as (Hunt, 1978)
Cc.-- /o"Cidt •
and the limit of Equation 12 as t •»• « is
c o / *
V X
Equations 14 and 16 describe the transient and steady-state concentration
distributions arising from a continuous point source in an infinite aquifer
with uniform ground-water flow.
Assumptions and Limitations
Equations 14 and. 16 can be used, to calculate the concentrations in a
leachate plume under the following assumptions and limitations:
1. The ground-water flow regime is completely saturated.
2. All aquifer properties are constant and uniform throughout the
aquifer.
3. The ground-water flow is horizontal, continuous, and uniform
throughout the aquifer.
4. The aquifer is infinite in extent.
5. The leachate source is a point located at the origin of the
coordinate system.
6. The mass flow rate of the source is constant.
7. At zero time the concentration of leachate in the aquifer is zero.
The assumptions of an infinite aquifer depth and a uniform source mass
rate can be overcome by using the principles of superposition in space and
-------
time, respectively (Walton, 1962). Both of these provisions have been
incorporated in the computer program described in the next section.
Superposition
The differential equation describing component mass concentration in a
porous medium, Equation 1, is a linear partial differential equation. The
principal of superposition can be used directly to solve complex ground-water
contamination problems in terms of the simplier solutions described above.
Unfortunately,'the scattered applications of this principle are not explained
in any single reference. Some texts indicate that superposition means that
any sum of solutions is also a solution.. Superposition is commonly used to
generate a linear no-flow boundary condition through the use of "image wells"
or to simulate multiple sources and sinks (Walton; 1962, 1970). The principle
of superposition is also complicated by referring to the "Duhamel theorem,"
the "Faltung integral," and/or "convolution integrals." These terms often
have no apparent physical interpretation. For the purposes of this report,
"superposition in space" will refer to the approximation of sources of finite
area or.volume as the sum of a finite number of point sources or the
generation of no-flow boundaries using image wells. "Superposition in time"
will refer to the approximation of a variable source rate of contamination as
the sum of a finite number of constant source rates distributed in time.
The three-dimensional solutions presented above can be used to simulate
aquifers of finite width or depth or sources of finite volume. Applications
of this type require a thorough understanding of the physical interpretation
of the principal of superposition.
However, some applications are relatively straight forward, and the
computer program provides for the approximation of a non-uniform source rate
-------
using superposition in time. Multiple sources and aquifers of finite
thickness are also included using superposition in space.
Consider the variable source of contamination shown in Figure 1. The
solutions of the governing differential equation presented in this report are
of the form
C(x,z,t) = CQQ' f(x,z,t) = Q1 f(x,z,t) (17)
where Q1 is the source mass rate per unit length. The principle of
superposition in time can be written for any position as
n
C(x,z,t) = E Q- f(x.,z,t-) (18)
^ 1
Now, the variable rate schedule shown in Figure la can be decomposed into a
series of positive and negative mass rates as shown in Figure Ib. The
concentration at a point x,y,z at the end of the simulation, ts, can be
evaluated as
• •
C(x,y,z,t) = QJ f(x,y,z,t1) - Q^ f(x,y,z,t2)
• •
+ Q£ f(x,y,z,t2) - Q£ f(x,y,z,t3)
+ Q3 f(x,y,z,t3) - Q^ f(x,y,z,t4)
+ Q^ f(x,y,z,t4) (19;
-------
o
•
UJ QS
h- • ,
CC • ,
CO 1
CO
i
63
0 Q'
oT
i Q'4
CO
CO .
*-t
-Qi
\
• I
0
.
-
•
r~
TIME, t
(a)
.
•
.
•
i
~-t4—
* ^
""• *y
*
p* \*£ •-
PERIOD OF
SIMULATION
TIME, t
(b)
Figure 1. Decomposition of a variable source rate using superoosition
in time.
-------
In general terms
C(x,y,z,ts) = £ (Qj - Qj_i) f(x,y,z,ti) (20)
•
with Q0 = 0
Note the time corresponding to a given source rate, t^, is the period
beginning with the start of the given rate to the end of the simulation
period; time is not the duration of a given rate. For ease of application,
Equation 20 can be rewritten as
C(x,y,z,t) = £ (Q1, - Q' ,) f(x,y,z,t -t, ,) (21)
b I/—1 N~l J N •" 1
K — 1
•where t^.} .is the time corresponding to the end of mass rate Q^.j or the
beginning of rate Q^ with Q0 = 0 and tQ = 0.
A continuous non-uniform rate schedule may be approximated as closely as
desired by increasing the number, of discrete rates in the source rate .
schedule. In theory an infinite number of discrete rates would be required.
An understanding of the physical problem and the assumptions incorporated in
the mathematical model are the best guidelines for decomposing a continuous
non-uniform source of contamination.
The influence of geohydrologic boundaries on the movement of a tracer is
similar to the influence of these boundaries on the drawdown response of an
aquifer to pumping. The applications of image well theory described by Walton
(1962, 1970) can be extended to the horizontally-averaged solution to the
solute transport problem considered in this report. The following discussion
11
-------
parallels Walton's examples of the use of image wells to account for barrier
boundaries.
Consider the contaminant plume which would exist if the aquifer were of
infinite depth as shown in Figure 2a. If the contaminant plume was to
intersect an impermeable base of the aquifer as shown in Figure 2b, the
vertical concentration gradient must change since there can be no transport of
mass across the boundary as a result of dispersion. In mathematical terms
Dz!=0
at z = B. Now, if an imaginary, or image, source were placed across the
boundary at a distance equal to the depth of the aquifer, as shown in
Figure 2c, this source would create a concentration gradient from the boundary
to the image, water table equal to the concentration gradient from the boundary
to the real water table. A "concentration divide" would be established at
boundary, and the no-transport boundary condition (3C/3Z = 0) would be
satisfied. . •
The imaginary system of a contaminant source and its image in an aquifer
of infinite depth satisfies the boundary conditions dictated by the finite
depth system. The resultant concentration distribution is the sum of
concentrations in both the real and image systems as shown in Figure 2d.
In theory an infinite number of image systems may be required. For
example, if the plume in the infinite system intersects the water table in the
image system a second no-transport boundary is encountered as shown in Figure
3. This boundary can be handled by introducing another image system across
the imaginary boundary and equidistant from the first image system. This
process of adding image systems could be repeated indefinitely. In practice
12
-------
WATER TABLE
WATER TABLE
(a)
(b)
(d)
Figure 2. Use of image sources to account for aquifers of finite depth.
13
-------
INFINITE
AQUIFER
c(x, z)
IMAGE 1st
CONCENTRATION IMAGE SYSTEM
[WATER TABLE
2nd
IMAGE SYSTEM
2d
c(x, z)
oc(x, 2d-z)
[MAG_E
WATER TABLE
4d
J^M AGIN A R Y _
BOUNDARY
oc(x, 2d+z)
\\\\\\\ V V\\VV\ \\\ \\V\\
c(x, 4d-z)
6d
c'(x, z) = c(x,
c(x,
f ffi t f / / iff / f / i
oc(x, 6d-z)
n=1
, 2nd-z)-»-c(x, 2nd+z)]
Figure 3. Superposition in space to account for barrier boundaries.
-------
only a few image systems are required. The computer program automatically
introduces an appropriate number of image systems.
15
-------
SECTION II
COMPUTER PROGRAM
The computer program evaluates the analytical solution of the
differential equation describing concentration distributions in a three-
dimensional plume with uniform ground-water flow. The program has been
designed for interactive use and requires input data under two modes of
operation — "Basic Input Data" and "Edit."
Basic Input Data
Basic input data are required to initiate a new problem using the PLUME3D
program. The user is prompted for the required data through a series of input
commands described below. Numeric data may be entered through the keyboard
with or without decimal points and multiple data entries should be separated
by comma(s). The first basic input command is:
ENTER TITLE
?
Any valid keyboard characters can be used. The first 60 characters will be
retained for further problem identification.
The next three input commands define the units for all variables used in
the calculations. Any consistent set of units may be used.
ENTER UNITS FOR LENGTH (2 CHARACTERS)
?
Any valid keyboard characters can be used. The first two characters will be
retained for identifying the units of the length dimension which may be
required for other input data or output listings.
ENTER UNITS FOR TIME (2 CHARACTERS)
7
16
-------
Any valid keyboard characters can be used. The first two characters will be
retained for identifying the units of the time dimension which may be required
for other input data or output listings.
ENTER UNITS FOR CONCENTRATION (6 CHARACTERS)
?
The first six characters of any valid keyboard entries will be retained for
identifying the concentration units for data input and output.
The remaining input commands are used to initiali-ze all variables for a
given problem. They include both aquifer and contaminant parameters. Input
data errors which may interrupt the computational sequence are detected by the
program, and a command is issued to reenter the data for the appropriate
variable.
ENTER SATURATED THICKNESS, (0. FOR INFINITE THICKNESS), L
?
The saturated thickness must be entered in the units requested with dimensions
of L. If a zero or negative value is entered, the calculations will be
carried out assuming an aquifer of infinite depth. The program automatically
includes up to 20 image wells for aquifers of finite depth.
ENTER AQUIFER POROSITY
?
Enter the volume void fraction.
ENTER SEEPAGE VELOCITY, L/t
?
The seepage, or interstitial, velocity must be entered with dimensions of L/t
in the units requested. Numerical values must be greater than zero.
ENTER RETARDATION COEFFICIENT
?
The retardation coefficient includes the effects of absorption of the tracer
on the solid matrix. The numerical value must be greater than 1.0, or equal
to 1.0 if absorption is neglected.
17
-------
ENTER X DISPERSION COEFFICIENT, SQ L/t
?
Dispersion coefficients have dimensions of L2/t and must be entered in the
units requested. Numerical values must be greater than zero. The next two
commands will ask for the Y and Z dispersion coefficients, respectively. They
also have dimensions of L2/T and must be entered in the units requested.
Numerical values must be greater than zero.
ENTER Y DISPERSION COEFFICIENT, SQ L/t
?
ENTER Z DISPERSION COEFFICIENT, SQ L/t
?
The subsequent command .is:
ENTER DECAY CONSTANT, 1/t
?
The first order decay constant has dimensions of 1/t and must be entered in
the units requested. The decay constant must be greater than, or equal to,
zero.
SELECT TRANSIENT OR STEADY-STATE SOLUTION
TR FOR TRANSIENT SOLUTION
SS FOR STEADY-STATE SOLUTION
?
Selection of the transient solution also allows the approximation of a
nonuniform rate schedule by a series of uniform rates. Approximation is
accomplished through superposition of a series of uniform rates. If steady-
state solution is chosen, the steady state concentration will be evaluated.
ENTER THE NUMBER OF SOURCES (MAXIMUM OF N)
?
The number of sources of contaminant should be entered. The value entered
must be greater than zero.
MASS RATES HAVE UNITS OF (M/L3) (L3/t)
TIME HAS UNITS OF t
18
-------
This statement reminds the user of the units that will be used for mass rates
and for time. All mass rates and time values entered must be in these units.
The next three commands will be repeated for each source.
ENTER X, Y, AND Z COORDINATES OF SOURCE I (L)
?
The input units for the coordinates must be in the units requested. The Z-
coordinate must be greater than or equal to zero.
If the transient solution was chosen the following two commands will be
issued.
ENTER THE NUMBER OF RATES FOR SOURCE I (MAXIMUM OF N)
?
The number of uniform rates used to approximate a nonuniform rate schedule for
this source is entered. The value must be greater than zero.
SOURCE I, RATE J STARTS AT TIME t
ENTER MASS RATE AND ENDING TIME
?
The source mass rate is entered in units of concentration times the volumetric
rate. Note the actual source concentration and rate are not required, but the
units must be consistent. The time units must also be consistent.
If the steady-state solution has been selected, the following command
will be entered instead of the two previously listed commands.
ENTER STEADY-STATE MASS RATE I
?
The next three basic input commands are used to define the matrix of
observation points, or coordinates at which concentration will be evaluated.
ENTER XFIRST, XLAST, DELTAX (L)
? ? ?
.,.,.
The input units for the coordinates must also be in the units requested. A
zero entry for DELTAX will result in a single X-coordinate observation.
19
-------
Results of calculations for multiple X-coordinates will be listed from XFIRST
to XLAST.
ENTER YFIRST, YLAST, DELTAY (L)
? 7 ?
• » • ^ •
Any of the numerical values used to define the Y-coordinates of
observation points may be positive or negative.
ENTER ZFIRST, ZLAST, DELTAZ (L)
?,?,?
Both ZLAST and ZFIRST must be greater than or equal to zero and less than or
equal to the saturated thickness.
ENTER PLANE FOR SECTIONING AQUIFER (XY, XZ, OR YZ)
?
The selection of a-particular plane determines the presentation of the output
of the 'program. Concentrations at the specified coordinates in the selected
plane will be printed for a constant value of the third coordinate.
ENTER TFIRST, TLAST, DELTAT (t)
? ? ?
. • > • 9 •
The beginning value and ending value of the time interval of contaminant
transport being modeled is entered. Both TFIRST and TLAST must be positive
values in the units requested. A zero entry for DELTAT will result in model
output at a single value of time.
Edit Commands
Once the basic input data have been entered, the problem as currently
defined is listed and the program enters the "edit" mode. The edit commands
are listed in Table 1 and are also listed the first time the program enters
the edit mode. The request for information is:
ENTER NEXT COMMAND ?
20
-------
One of the reponses from Table 1 should be given. If the response is
incorrect or improperly formulated the statement
ERROR IN LAST COMMAND — REENTER ?
is issued. Error messages for invalid numerical data will be issued as
described under the Basic Input Commands. The request for information will be
repeated until one of the responses Mil, LI, RN, NP, or DN is entered.
MU will list the table of edit commands.
LI will list the problem as currently defined.
RN will initiate the calculation of concentration-s and print the results.
NP will request a complete new problem .using the "Basic Input Data"
dialog.
DN will terminate the program.
Although many tests for valid input data and properly formulated edit
commands, have been embedded in the program, the user is 'encouraged to correct
"keyboard errors" before the data are transmitted. These precautions will
serve to minimize the frustration of program termination as a result of fatal
errors during execution of the numerical computations.
21
-------
Table 1
EDIT COMMANDS
Command Variable changed/Execution
ST . Saturated thickness
PO Porosity
VX Seepage Velocity
RD Retardation Coefficient
DE Decay Constant
DX . X-Dispersion Coefficient
DY Y-Dispersion Coefficient .
DZ Z-Dispersion Coefficient
RT Source-Rate Schedule
OB Observation Points
XC X-Coordinates
ZC Z-Coordinates
YC Y-Coordinates
TC Observation Times
AS Aquifer Sectioning
CS Change Solution/Sources
MU Menu of Edit Commands
LI List Input Data
RN Run
NP New Problem
DN Done
22
-------
SECTION III
APPLICATIONS
The example problems presented in this document are based on the
dispersal of chromium discharged to the ground water in southeastern Nassau
County, New York. The hydrogeology and history of contamination have been
documented by Perlmutter and Lieber (1970) and will be briefly summarized in
the following paragraphs.
Site Location
The general area of the documented case history of the concentrations of
contaminants in ground-water is in southeastern Nassau County, Long Island,
New .York (Figure 4). The detailed study area included an .industrial park at
South Fanningdale. During World War II, the 'industrial park was occupied by
an aircraft company whose cadmium and chromium enriched metal plating waste
was the source of much of the heavy-metal contamination in a shallow glacial
aquifer.•
Hydrogeology
The upper glacial aquifer which is addressed in this document extends
from the water table, at depths ranging from 0 to 15 feet, below the ground
surface, to the top of the deeper Magothy aquifer, at depths of 80 to 140 feet
below the ground surface. The upper unit consists of beds and lenses of fine
to coarse sand and gravel. In some parts of the aquifer, thin lenses of fine
to medium sand, as well as some silt, are interbedded with the coarse
23
-------
material. Data from scattered well borings indicate that the lower 8 to 10
feet of the upper unit may consist of silty and sandy clays.
The principal direction of ground-water flow in the upper glacial unit is
from north to south. The regional hydraulic gradient is approxmately 0.0025
o
ft/ft and the average coefficient of permeability is assumed to be 1*600 gal'/day/ftc
with an average total porosity of 0.35 (Perlmutter and Lieber, 1970). The
ground-water movement is horizontal throughout the area of interest, with the
exception of local recharge and discharge areas where vertical or oblique flow
may be predominant.
Plating-Waste Contamination
The chromium contamination in the upper glacial unit was derived
principally from the disposal of metal-plating and anodizing waste water to
unlined basins. During World War II, and.for several, years thereafter,
essentially untreated plating-water effluent was recharged to the- aquifer
through these disposal basins. Perlmutter and Lieber (1970) estimated that
during the early 1940's as much as 200,000 to 300,000 gallons of effluent
containing approximately 52 pounds of chromium was discharged daily into the
upper glacial unit. At the end of World War II, the amount of plating-waste
effluent was decreased substantially and a chromium removal unit was installed
in 1949.
By 1949, or approximately 9 years after the start of plating-waste
disposal, a plume of contaminated water had-migrated southward approximately
3,900 feet with a maximum width of about 850 feet. Maximum observed chromium
concentration was approximately 40 mg/1.
Since 1949 the chromium concentrations decreased substantially in nearly
all sections of the sections of the plume. The maximum observed chromium
24
-------
concentration was about 10 mg/1 in 1962. However, the plume of contaminated
water had migrated further southward beyond Massapequa Creek, a small stream
which serves as a natural drain for part of the contaminated ground water.
The plume was approximately 4,300 feet long with a maximum width of almost
1000 feet in 1962.
Analytical Model
The data required for the analytical two-dimensional plume model are
listed in Table 2. The following paragraphs will discuss the estimation of
•these model parameters based on the limited field data, judgement, and
experience.
All input data for the computer program must be in consistent units. For
the example problems, the following system of units will be selected: length
in feet, time in days, and concentration in mg/1.
' Assuming that the water table is located approximately at the ground
surface, the average saturated thickness of the upper glacial unit can be
estimated from the depth to the top of the deeper Magothy aquifer. Thus, the
average saturated thickness, St, will be taken as:
st - 8Q * 14° - 110 ft
The superficial, or Darcy, velocity can be estimated from the average
coefficient of permeability, K, and the regional hydraulic gradient, dh/dx.
From Darcy's law
V = - K.
dx
25
-------
TABLE 2
Input Data Required for the Analytical
Three-Dimensional Plume Model
Title - Units for length, time, and concentration
Saturated thickness (for aquifer of finite depth)
Effective porosity
Ground water interstitial velocity
Retardation coefficient
Longitudinal, dispersion coefficient
Transverse dispersion coefficient
Vertical dispersion coefficient
First-order decay constant
Type of solution (transient or steady-state)
Number of sources
Location and .rate schedules for each source
Coordinates of observation points
Observation times (for transient solution)
26
-------
where the minus sign indicates that the flow is in the direction of decreasing
hydraulic head. Thus
v _ 1600 gal ft3- 0.0025 ft
~ TF" 7.48 gal ft
day ft 3
and • ' .
V = 0.52 ft/day
The average interstitial, or pore, velocity, V*, can be estimated by assuming
that the area! porosity is equal to the volume porosity, 0, and
Using the estimated value of effective porosity
or
= 0.52 ft/day
OB"
V* = 1.5 ft/day
The longitudinal, transverse and vertical dispersion coefficients are the
most- difficult parameters to estimate. Typical values for dispersivities, as
well as hydraulic conductivity, are summarized in Table 3. Based on the
results of a numerical model for the same site, Pinder (1973) estimated
longitudinal and transverse dispersivities as ctx = 69.9 feet and ay = 14.0
feet, respectively. Considering the interbedding of finer materials in the
27
-------
TABLE 3
Typical Values of Aquifer Properties
(after Yeh, 1981)
Parameter
Bulk density, lb/ft3
Effective porosity
Hydraulic Conductivity,
gal/day/ft2
Dispersivity, ft
Longitudinal
Transverse
Vertical
Clay
87.36 - 137.2
0.03 - 0.05
0.01 - 0.1
0.1 - 1.0
0.01 - 0.1
0.01 - 0.1
Material
Silt
80.50 - 112.3
0.05 - 0.10
1 - 10
1 - 10
0.1 - 1.0
0.1 - 1.0
Sand
73.63 - 98.59
0.10 - 0.30
100 - 100,000
10 - 100
1.0 - 10
1 - 10
28
-------
shallow upper unit together with the typical values of longitudinal and
vertical dispersivities in Table 3, the vertical dispersivity will be assumed
to be approximately two orders of magnitude smaller than the longitudinal
dispersivity. Thus oz = 0.7 ft, and the dispersion coefficients are evaluated
as
Dv = a, V* = (69.9 ft)(1.5 ft/day) = 105 ft2/day
A I—
D = a V* = (14.0 ft)(1.5 ft/day) = 21.0 ft2/day
and
DZ = az V* = (0.7 -ft) ('1.5 ft/day) = 1.05'ft2/day.
Chromium is believed to be a conservative material in this system
(Perlmutter and Lieber, 1970), and both adsorption and chemical/biological
decay can be neglected. Therefore, Kd = 0 (or Rd = 1) and x = 0.
The last set of model input data to be specified is the source mass/rate
schedule. A steady mass rate of 52 Ib/day of chromium will be assumed from
the time of initial injection to mid-1949, or approximately 2800/days. At
this time, the actual extent of the contamination was estimated; and the
chromium removal unit was installed. The actual volumetric rate and
concentration of the source of chromium contamination are not required.
However, the units of the source term must be consistent with other parameters
in the model.
The mass rate of chromium can be converted to units of concentration
times volume rate per unit width of aquifer as follows:
29
-------
TABLE 4
Model Input Data for Example Problem 1.-
Title: Hexavalent Chromium Plume—Example 1
Units for length: ft
Units for time: dy
Units for concentration: mg/1
Saturated thickness: 110 ft
Effective porosity: 0.35
Interstitial velocity: l.b ft/dy
Retardation coefficient: 1.0
x-dispersion coefficient 105 ft/dy
y-dispersion coefficient 21.0 ft^/dy
o
z-dispersion coefficient 1.05 ft/dy
First-order decay constant: . 0/dy
Source/rate schedule
Number of sources: 1
Location of source: x=0, y=0, z=0
Number of rates: 1
Mass rate and time: 833,586 (mg/1)(ft3/dy) from 0 to 2800 dy
30
-------
52 Ib 454 x 1Q3 mg ft3 Q,, ,-QC / n .., w^3.. >
Ifay Tb 28.321 liter = 833>586 (mg/^ter)(ft /day)
for 2800 days.
The model input data for this example problem are summarized in Table 4.
The input data dialog for the example problem and the model results are
included in Appendix A.
This example problem is not intended to serve as a verification of a
model. The example illustrates the type of input data required, and some
methods for estimating input data. The input data dialog and model results
can be used to partially test the model code. The analytical model and the
computer program are tools which can aid in the analysis of ground-water
contamination problems-. The user must select the best tools for the problem
at hand based on a sound .understanding of the principles-of ground-water
hydrology, the physical problems, and the limitations of the mathematical
model(s).
31
-------
REFERENCES
Abramowitz, M. and I. A. Stegan, 1966, Handbook of Mathematical Functions with
Formulas, Graphs, and Mathematical Tables, National Bureau of Standards
Applied Mathematics Series 55, U. S. Department of Commerce, 1046 pp.
Cars!aw, H. S. and J. C. Jaeger, 1959, Conduction of Heat in Solids, Oxford
University Press, Oxford-, England, 510 pp.
Hunt, B., 1978, "Dispersive Sources in Uniform Ground-Water Flow, Jour.
Hydraul. Div.. ASCE, 104, No. HY1, January 1978, pp. 75-85.
Perlmutter, N. M. and M. Lieber, 1970, "Dispersal of Plating Wastes and Sewage
Contaminants in Ground Water and Surface Water, South Farmingdale-
Massapequa Area, Nassau County, New York," U. S. Geological Survey Water-
Supply Paper 1879-G, pp. G1-G67.
Pinder, G. F., 1973, "A Galerkin Finite Element Simulation of Ground-water
Contamination on Long Island," Water Resources Research, Vol. 9, No. 6,
pp. 1657-1669.
Turner, G. A., 1972, Heat and Concentration Waves, Academic Press, New York,
New-York, 233 pp.
Walton, W. C., 1962, "Selected Analytical Methods for Well and Aquifer
Evaluation," Bulletin 49, Illinois State Water Survey* Urbana, Illinois,
81 pp.
Walton, W. C., 1970, Groundwater Resource Evaluation, Mc-Graw-Hil1, New York,
New York, 664 pp.
Yen, G. T., 1981, "AT123D: Analytical Transient One-, Two-, and Three-
Dimensional Simulation of Waste Transport in the Aquifer System,"
Publication No. 1439, Environmental Sciences Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee, 83 pp.
32
-------
APPENDIX A
Input Data Dialog and. Results for Example Problem
A-l
-------
ENTER TITLE
7HEXAVALENT CHROMIUM PLUME — EXAMPLE 1
'ENTER UNITS FOR LENGTH (2 CHARACTERS)
?FT
ENTER UNITS FOR TIME (2 CHARACTERS)
?DY
ENTER UNITS FOR CONCENTRATION (6 CHARACTERS)
7MG/L
ENTER SATURATED THICKNESS (0 FOR INFINITE THICKNESS), FT
7110.
ENTER AQUIFER POROSITY
?0.35
ENTER SEEPAGE VELOCITY', FT/DY
ENTER RETARDATION COEFFICIENT
ENTER X DISPERSION COEFFICIENT, SQ FT/DY
?105.
ENTER Y DISPERSION COEFFICIENT, SQ FT/DY
721.
ENTER Z DISPERSION COEFFICIENT, SQ FT/DY
71.05
ENTER DECAY CONSTANT, 1/DY
70.
SELECT TRANSIENT OR STEADY-STATE SOLUTION
TR FOR TRANSIENT SOLUTION
SS FOR STEADY-STATE SOLUTION
7TR
A-2
-------
ENTER THE NUMBER OF SOURCES (MAXIMUM OF 10 )
MASS RATES HAVE UNITS OF (MG/L ) (CU FT/DY)
TIME HAS UNITS OF DY
ENTER X, Y, AND Z COORDINATES OF SOURCE 1 (FT)
?,?,?0. ,0. ,0.
ENTER THE NUMBER RATES FOR SOURCE 1 (MAXIMUM 0F 10)
SOURCE 1, RATE 1 STARTS AT 0.0 DY
ENTER MASS RATE AND ENDING TIME
?,7833586.,2800.
ENTER XIRST, XLAST, DELTAX (FT)
7, 7,7600.,3600.,600.
ENTER YFIRST, YLAST, DELTAY (FT)
.?,?,?450.,-450.,150.
ENTER ZFIRST, ZLAST, DELTAZ (FT)
?,?,?0. ,110.,55.
ENTER PLANE FOR SECTIONING AQUIFER (XY,XZ,OR YZ)
?XY
ENTER TFIRST, TLAST, DELTAT (DY)
?,?,72800.,2800. ,0.
A-3
-------
PLUME3D
VERSION 2.01
PAGE 1
HEXAVALENT CHROMIUM PLUME — EXAMPLE 1
SATURATED THICKNESS, (FT)
SEEPAGE VELOCITY, (FT/DY)
X DISPERSION COEFFICIENT (FT**2/DY)
Y DISPERSION COEFFICIENT (FT**2/DY)
Z DISPERSION COEFFICIENT (FT**2/DY)
POROSITY
110.0000
1.5000
105.0000
21.0000
1.0500
.3500
RETARDATION COEFFICIENT
FIRST ORDER DECAY CONSTANT (1/DY)
1.0000
0.0000
SOURCE/RATE SCHEDULE (MG/L )(CU FT/DY)
SOURCE . RATE MASS
NO X (FT) Y (FT) Z (FT) NO . RATE
TIME (DY)
START END
0 .00
0.00
0.00
833586 .00
0.00 2800.00
OBSERVATION POINTS (FT), AND TIMES (.DY)
XFIRST =
YFIRST =
ZFIRST =
TFIRST =
600.00
450.00
0.00
2800.00
XLAST =
YLAST =
ZLAST =
TLAST =
3600.00
-450.00
110.00
2800.00
DELX =
DELY =
DELZ =
DELT =
600.0000
150.0000
55.0000
0.0000
AQUIFER SECTIONED IN XY PLANE
A-4
-------
MENU OF EDIT COMMANDS
SATURATED THICKNESS ST
POROSITY PO
SEEPAGE VELOCITY . VX
RETARDATION COEFFICIENT RD
X DISPERSION COEFFICIENT DX
Y DISPERSION COEFFICIENT DY
Z DISPERSION COEFFICIENT DZ
DECAY CONSTANT . DE
SOURCE/RATE SCHEDULE RT
.CHANGE SOLUTION/SOURCES CS
OBSERVATION POINTS OB
X COORDINATES XC
Y COORDINATES YC'
Z COORDINATES ZC
OBSERVATION TIMES TC
AQUIFER SECTIONING AS
NEW PROBLEM NP
MENU OF COMMANDS MU
LIST INPUT DATA LI
RUN CALCULATIONS RN
DONE DN
ENTER NEXT COMMAND
?RN
A-5
-------
PLUME3D
VERSION 2.01
PAGE 2
HEXAVALENT CHROMIUM PLUME — EXAMPLE 1
CONCENTRATION DISTRIBUTION AT 2800.00. DY (MG/L )
Z = 0.00 FT
* X(FT)
*
Y(FT) *
600.00
1200.00
1800.00 2400.00
3000.00
3600.00
450.00
300.00
150.00
0.00
-150.00
-300.00
-450.00
1
10
62
134
62
10
1
.1622
.5229
.9100
.5398
.9100
.5229
.1622
3
16
46
67
46
16
3
.7737
.8523
.6486
.2738
.6486
.8523
.7737
6
17
35
44
35
17
6
.0164
.7165
.3420
.8561
.3420
.7165
.0164
7
16
28
33
28
16
7
.2392
.6967
.0852
.•5146
.0852
.6967
.2392
7
14
22
25
22
14
7
.3914
.7097
.4262
.8517
.4262
.7097
.3914
6.2020
11 .3227
16.3152
18.4413
16.3152
11.3227
. 6.2920
A-6
-------
PLUMESD .
VERSION 2.01
PAGE 3
HEXAVALENT CHROMIUM PLUME — EXAMPLE 1
CONCENTRATION DISTRIBUTION AT 2800.00 DY (MG/L )
Z = 55.00 FT
*
* X(FT)
* 600.00 1200.00 1800.00 2400.00 3000.00 3600.00
Y(FT) *
450.00
300.00
150.00
' 0.00
-150.00
-300.00
-450.00
2
12
21
12
2
.4395
.9774
.3888
.5268
.3888
.9774
.4395
1
7
18
26
18
7
1
.8379
.3832
.7138
.0413
.7138
.3832
.8379
3-.
9.
19.
24.
19.
9.
' 3.
5409
9691
2646
1684
2646
9691
5409
.4
10
18
21
18
10
4
.8545
.9487
.1440
.5383 •
.1440
.9487
.8545
5
10
16
. 18
16
10
5
.4148
.6514.
.1159
.5286
.1159
.6514'.
.4148
4.8044
8.7202
12.5181
14.1310
12.5181
8.7202
4.8044
A-7
-------
PLUME3D
VERSION 2.01
PAGE 4 .
HEXAVALENT CHROMIUM PLUME — EXAMPLE 1
CONCENTRATION DISTRIBUTION AT 2800.00 DY (MG/L )
Z = 110.00 FT
* X(FT)
*
Y(FT) *
*
450.00
300.00
150.00
0.00
-150.00
-300.00
-450.00
600.00
.0755
.3219
.8561
1.2145
.8561
.3219
.0755
1200.00
.5147
1.6706
3.5798
4.6664 -
3.5798
1.6706
.5147
1800.00
1.4813
3.7892
6.8411
8.3725
6.8411
3.7892
1.4813
2400.00
2.6624
5.7406
9.2287
10.8380
9.2287
5.7406
2.6624
3000.00
3.5125
6.7651
10.0941
11.5486
10.0941
6.7651
3.5125
3600.00
3.4326
6.1702
8.8018
9.9142
8.8018
6.1702
3.4326
ENTER NEXT COMMAND
?YC
ENTER YFIRST, YLAST, DELTAY (FT)
? , ?,?0.,0.,0.
ENTER NEXT COMMAND
?ZC
ENTER ZFIRST, ZLAST, DELTAZ (FT)
ENTER NEXT COMMAND
?AS
ENTER PLANE FOR SECTIONING AQUIFER (XY,XZ,OR YZ)
?XZ
ENTER NEXT COMMAND
A-8
-------
PLUME3D
VERSION 2.01
PAGE 5
HEXAVALENT CHROMIUM PLUME — EXAMPLE 1
SATURATED THICKNESS, (FT)
SEEPAGE VELOCITY, (FT/DY)
X DISPERSION COEFFICIENT (FT**2/DY)
Y DISPERSION COEFFICIENT (FT**2/DY)
Z DISPERSION COEFFICIENT (FT**2/DY)
POROSITY
110.0000
1.5000
105.0000
21.0000
1.0500
.3500
RETARDATION COEFFICIENT
FIRST ORDER DECAY CONSTANT (1/DY)
1.0000
0.0000
SOURCE/RATE SCHEDULE (MG/L )(CU FT/DY)
SOURCE , RATE MASS
NO X (FT) Y- (FT) . Z (FT) MO RATE
TIME (DY)
START' . END
0.00
0.00
0.00
833586.00
0.03 2800.00
OBSERVATION POINTS (FT), AND TIMES (DY)
XFIRST =
YFIRST =
ZFIRST =
TFIRST =
600.00
0.00
0.00
2800.00
XLAST =
YLAST =
ZLAST =
TLAST =
3600.00
0.00
110.00
2800.00
DELX =
DELY =
DELZ =
DELT =
600 .0000
0 .0000
20.0000
0.0000
AQUIFER SECTIONED IN XZ PLANE
ENTER NEXT COMMAND
?RN
A-9
-------
PLUME3D
VERSION 2.01
PAGE 6
HEXAVALENT CHROMIUM PLUME — EXAMPLE 1
CONCENTRATION DISTRIBUTION AT 2800.00 DY (MG/L. .)
Y = 0.00 FT
* X(FT)
*
Z(FT) *
*
0.00
20.00
40 .00
60.00
80 .00
100.00
110.00
600.00
134.5398
101.2264
47.1345
16.1375
4.7086
1.5140
1.2145
1200.00
67.2738
58.9650
40.1774
22 .0086
10.3607
5.2569
4.. 6664
1800.00
44.8561
41.2153
32.1269
21.6384
13 .3152
8.9261
8.3725
2400.00
33.5146
31.5259
26.3553.
19.9389
14.3904
11.2473
10.8380
3000 .00
25.8517
24.6661
21.5287
17.5144
13.9192
11.8248
. 11 .5486
3600 .00
18.4413
17.7508
15.9102
13 .5252
11.3590
10.0831
9.9142
ENTER NEXT COMMAND
?DN
STOP
A-10
-------
APPENDIX B
Listing of Source Code for PLUME3D
B-l
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
PLUME3D
VERSION 2.02
THREE-DIMENSIONAL PLUMES IN UNIFORM GROUND-WATER FLOW
JAN WAGNER
SCHOOL OF CHEMICAL ENGINEERING
OKLAHOMA STATE UNIVERSITY
STILLWATER, OK 74078
PHONE (405) 624-5280
MARCH, 1984
REVISIONS:
2.01
2.02
23 NOV 84
9 DEC 84
DIMENSION COL(7).CON(7),D(3),DEL(3),XF(3),XL(3),XS(10),YS(10),
1ZS(10)
DIMENSION IC(21),IS(4),KSEC(3),LBL(2.6),NP(3).NR(10),TITLE(30)
REAL LAMBDA
INTEGER TITLE.
COMMON/IO/NI.NO
COMMON/RATE/QC10,12),T(10,12)
COMMON/PHYPRO/ALPHA.BETA,GAMMA,DX,LAMBDA,PE,RD,V
DATA IC/'NP','ST','PO','VX','RD','DX','DY'.'DZ','DE'
1'YC'.'ZC','AS'.'TC','LI','RN','OB','RT','MU','DN'/
DATA IS/'R','M','A','D'/
DATA KAR1,KAR2.KAR3/'X','Y','Z'l
DATA KSEC/'XY'.'XZ','YZ'/
DATA LBL/' ','(C'.' '.'ON',' ','TI','
1' '.')'/
DATA IY/'Y'/
DATA KSOL1,KSOL2/'TR','SS'/
'NU' , '
READ DEVICE: NI
NI = 1
N0=1
WRITE DEVICE: NO
MAXIMUM NUMBER OF PRINTED COLUMNS PER PAGE IS SET TO MAXCOL
DIMENSION COL(MAXCOL).CON(MAXCOL)
MAXCOL = 7 '
MAXIMUM NUMBER OF PRINTED ROWS PER PAGE IS SET TO MAXROW
MAXROW = 40
MAXIMUM'NUMBER OF SOURCES IS SET TO MAXSOR
DIMENSION XS(MAXSOR),.YS(MAXSOR),ZS(MAXSOR),NR(MAXSOR)
MAXSOR = 10
MAXIMUM NUMBER OF SOURCE RATES FOR SUPERPOSITION IN TIME
IS SET TO MAXRT
COMMON/RATE/ 0(MAXSOR.MAXRT+2),T(MAXSOR.MAXRT+2)
MAXRT = 10
MAXIMUM NUMBER OF IMAGE WELLS FOR.SUPERPOSITION IN SPACE
IS SET TO MAXIMG
MAXIMG = 20
INITIALIZE PROGRAM FLOW CONTROL VARIABLES
10 IEDIT = 1
KNTL = 1
NPAGE = 1
SECTION I -- BASIC INPUT DATA
READ TITLE
WRITE(NO,15)
15 FORMAT(1H1.2X,'ENTER TITLE',/'
READ(NI,25) (TITLE(I). 1=1,30)
25 FORMAT(30A2)
PL3D001
PL3D002
PL3D003
PL3D004
PL3D005
PL3D006
PL3D007
PL3D008
PL3D009
PL3D010
PL3D011
PL3D012
PL3D013
PL3D014
PL3D015
PL3D016
PL3D017
PL3D018
PL3D019
PL3D020
PL3D021
PL3D022
PL3D023
PL3D024
PL3D025
PL3D026
PL3D027
PL3D028
PL3D029
PL3D030
PL3D031
PL3D032
.PL3D033
PL3D034
PL3D035
PL3D036
PL3D037
PL3D038
PL3D039
PL3DO40
PL3D04-1
PL3D042
PL3D043
PL3D044
PL3D045
PL3DO46
PL3D047
PL3D048
PL3D049
PL3D050
PL3D051
PL3D052
PL3D053
PL3D054
PL3D055
PL3D056
PL3D057
PL3D058
PL3D059
PL3D060
PL3D061
PL3D062
PL3D063
PL3D064
PL3D065
PL3D066
PL3D067
PL3D068
PL3D069
PL3D070
5-2
-------
DEFINE UNITS
WRITE(NO,35)
35 FORMATOX, 'ENTER UNITS FOR LENGTH (2 CHARACTERS)',/,' ?')
READ(NI,45) IL
45 FORMAT(A2)
WRITE(NO,55)
55 FORMATOX,'ENTER UNITS FOR TIME (2 CHARACTERS)',/,' ?')
READ(NI,45) IT
WRITE(NO,65)
65 FORMATOX,'ENTER UNITS FOR CONCENTRATION (6 CHARACTERS)',/,' '
READ(NI.75) IM1.IM2.IM3
75 FORMAT(3A2)
ENTER .DATA FOR FIRST PROBLEM
SATURATED THICKNESS
80 IMAGE = MAXIMG/2
84 WRITE(NO,85) IL
85 FORMATOX,'ENTER SATURATED THICKNESS (0 FOR INFINITE THICKNESS),
1A2,/, ' ?' )
READ(NI,95,ERR=84) ST
95 FORMAT(F10.0)
IF(ST.GT.O.O) GO TO 100
IMAGE = 1 .
ST = 1.OE32
100 CONTINUE
GO TO '( 120. 110) . IEDIT
110 .IF(-XF(3) .LE.ST.AND.XLO) .LE.ST) GO TO 3000
WRITE(NO,115) ST;IL
115 FORMATOX.'RANGE OF Z-COORDINATES IS OUTSIDE UPPER LIMIT OF',/,
13X,'SATURATED THICKNESS = '.F10.4.A3).
GO TO 59O
'POROSITY
120 WRITE(NO.125)
125 FORMATOX. 'ENTER AQUIFER POROSITY',/,' ?')
130 READ(NI,95,ERR=120) P
IFIP.GT.O.O.AND.P.LT.1.0) GO TO 150
WRITEfNO,145)
145 FORMATOX,'POROSITY MUST BE GREATER THAN ZERO',
1' AND LESS THAN ONE -- REENTER',/,' . ?')
GO TO 130
150 GO TO (160,3000).IEDIT
SEEPAGE VELOCITY
160 WRITE(NO,.165) IL.IT
165 FORMAT(3X, 'ENTER SEEPAGE VELOCITY, 'A2. '/' .A2./. ' ?')
170 READ(NI,95,ERR=160) V
IF(V.GT.O.O) GO TO 180
WRITE(NO,175)
175 FORMATOX,'SEEPAGE VELOCITY MUST BE GREATER THAN ZERO',
1' -- REENTER',/,' ?' )
GO TO 170
180 GO TO (190,3000),IEDIT
•RETARDATION COEFFICIENT
190 WRITE(NO,195)
195 FORMATOX,'ENTER RETARDATION COEFFICIENT',/,' ?')
200 READ(NI,95,ERR=190) RD
IF(RD.GE.1.0) GO TO 21O
WRITE(NO,205)
205 FORMATOX,'RETARDATION COEFFICIENT MUST BE GREATER THAN OR',
1' EQUAL TO ONE',/,3X.' -- REENTER',/,' ?')
GO TO 200
210 GO TO (220,3000),IEDIT
X DISPERSION COEFFICIENT
220 WRITE(NO,225) IL.IT
225 FORMATOX,'ENTER X DISPERSION COEFFICIENT. SO ',A2,
PL3D071
PL3D072
PL3D073
PL3D074
PL3D075
PL3D076
PL3D077
PL3D078
PL3D079
PL3D080
) PL3D081
PL3D082
PL3D083
PL3D084
PL3D085
PL3D086
PL3D087
PL3D088
PL3DO89
' PL3D090
PL3D091
PL3D092
PL3D093
PL3D094
PL3D095
PL3D096
PL3D097
PL3D098
PL3D099
PL3D100
PL3D101
PL3D102
PL3D1C3
PL3D104
PL3D105
PL3D106
PL3D107
PL3D108
• PL3D109
PL3D110
PL3D111
PL3D1 12
PL3D113
PL3D114
. PL3D115
PL3D116
PL3D1 17
PL3D1 18
PL3D119
PL3D120
PL3D121
PL3D122
PL3D123
PL3D124
PL3D125
PL3D126
PL3D127
PL3D128
PL3D129
PL30130
PL3D131
PL3D132
PL3D133
PL3D134
PL3D135
PL3D136
PL3D137
PL3D138
PL3D139
PL3D140
3-3
-------
DISPERSION COEFFICIENT MUST BE GREATER THAN ZERO'
/, ' ?' )
230 READ(NI,95.ERR=220) DX
IF(DX.GT.0.0) GO TO 240
WRITE(NO,235)
235 FORMATOX,'X DISPERSION COEFFICIENT MUST BE GREATER THAN ZERO'
1' -- REENTER'./,' ?')
GO TO 230
240 GO TO (250,3000).IEDIT
C
C Y DISPERSION COEFFICIENT
250 WRITE(NO,255) IL.IT
255 FORMATOX,'ENTER Y DISPERSION COEFFICIENT, SO '.A2,
260 READ(Ni,95,ERR=250) DY
IF(DY.GT.0.0) GO TO 270
"WRITE(N0.265)
265 FORMATOX, 'Y
1' -- REENTER'
GO TO 260
270 GO TO (280.3000),I EDIT
C
C Z DISPERSION COEFFICIENT
280 WRITE(NO,285) IL.IT
285 FORMATOX, 'ENTER Z DISPERSION COEFFICIENT, SO
290 READ(NI,95,ERR=280) DZ
IF(DZ.GT.0.0) GO TO 300
•WRITE(NO,295)
. 295 FORMATOX. 'Z
1' -- REENTER'
GO TO 290
300 GO TO (310,3000).IEDIT
C
C FIRST-ORDER DECAY CONSTANT
310 WRITE(NO,315) IT
315 FORMATOX, 'ENTER DECAY CONSTANT,
READ(NI,95,ERR=310) DECAY
GO TO (320,3000).IEDIT
DISPERSION COEFFICIENT MUST. BE GREATER THAN ZERO'
/. ' ?' )
DEFINE LOCATIONS AND RATES OF SOURCES
INITIALIZE SOURCE/RATE ARRAYS
320 MAXRT2 = MAXRT + 2
DO 330 1=1.MAXSOR
DO 330 d=1.MAXRT2 •
0(1.J) = 0.0
T(I.J) = 0.0
330 CONTINUE
JFLOW = 3
340 WRITE(NO,345)
345 FORMATOX,'SELECT TRANSIENT OR STEADY-STATE SOLUTION',/,
16X,'TR FOR TRANSIENT SOLUTION',/.
26X.-SS FOR STEADY-STATE SOLUTION'./,' ?')
350 REAO(NI,45) KSOL
IF(KSOL.EQ.KSOLI) JFLOW=1
IF(KSOL.EO.KSOL2) JFLOW=2
GO TO (370,370,360), JFLOW
360 WRITE(NO,365)
365 FORMATOX,'ERROR IN SELECTION -- REENTER'/,' ?')
GO TO 350
370 WRITE(NO,375) MAXSOR
375 FORMATOX.'ENTER THE NUMBER OF SOURCES (MAXIMUM OF',13,'
1' ?')
380 REAO(NI,385,ERR=370) FDUM
385 FORMAT(FIO.O)
NS=FDUM
IF(NS.GT.O.AND.NS.LE.MAXSOR) GO TO 400
WRITE(N0.395) MAXSOR
395 FORMATOX.'NUMBER OF SOURCES MUST BE GREATER THAN ZERO '
1'AND LESS THAN',13,' -- REENTER',/,' ?')
PL3D141
PL3D142
PL3D143
PL3D144
PL3D145
PL3D146
PL3D147
PL3D148
PL3D149
PL3D150
PL3D151
PL3D152
PL3D153
PL3D154
PL3D155
PL3D156
PL3D157
PL3D158
PL3D159
PL3D160
PL3D161
PL3D162
PL3D163
PL3D164
PL3D165
. PL3D166
PL3D167
PL3D168
PL3D169
PL3D170
PL3D171
PL3D172
PL3D173
PL3D174
PL3D175
PL3D176
PL3D177
PL3D178
PL3D179
• PL3D180
PL3D181
PL3D182
PL3D183
PL3D184
PL3D185
PL3D186
PL3D187
PL3D188
PL3D189
PL3D190
PL3D191
PL3D192
PL3D193
PL3D194
PL3D195
PL3D196
PL3D197
PL3D198
PL3D199
PL3D200
PL3D201
PL3D202
PL3D203
PL3D204
PL3D205
PL3D206
PL3D207
PL3D208
PL3D209
PL3D210
B-4
-------
GO TO 380
400 WRITE (NO, 405) IM1 , IM2 , IM3 , IL, IT , IT
405 FORMATOX, 'MASS RATES HAVE UNITS OF (',3A2,') (CU '.A2.'/',A2.
1 ' )' ,/.3X, 'TIME HAS UNITS OF ',A2,/)
DO 540 1=1 ,NS
414 WRITE(NO,415) I.IL
415 FORMATOX, 'ENTER X. Y, AND Z COORDINATES OF SOURCE', 12,
1' (',A2. ')'./,' ?.?,?')
READ(NI,425,ERR=414) XS( I ) . YS( I ) , ZS( I )
425 FORMAT(3F10.0)
430 IF(ZS(I) .GE.O.O.AND.ZS(I) .LE.ST) GO TO 440
434 WRITE(NO,435) ST.IL
435 FORMATOX, ' Z-COORDINATE MUST BE GREATER THAN OR EQUAL TO ZERO',
1' AND' ,/,3X, 'LESS THAN OR EQUAL TO SATURATED THICKNESS (',
2F10.4.A3, ' )' ,/.3X, ' -- REENTER',/,' ?')
READ(NI,95,ERR=434) ZS(I)
GO TO 430
440 IF(JFLOW.E0.2) GO TO 530
0(1, 1 ) = 0.0
T(I, 1) = 0.0
450 WRITE(NO,455) I.MAXRT
455 FORMATOX, 'ENTER THE NUMBER RATES FOR SOURCE', 12,
1' (MAXIMUM OF' , 13, ')'./.' ?')
460 READ(NI , 465 , ERR=450) FDUM
465 FORMAT(FIO.O)
NR(I)=FDUM
IF(NR(I) .GT.O.AND.NR(I) .LE.MAXRT j GO TO 480
WRITE(NO,475) MAXRT
475 FORMATOX, 'NUMBER OF RATES MUST BE 'GREATER THAN ZERO AND ',
1'LE-SS THAN'. 13,' -- REENTER',/,' ?')
GO TO 460
480 CONTINUE
NRT = NR(I)
DO 52O J=1 .NRT
M = J + 1 ' .
WRITE(N0.485) I , J , T ( I . M- 1 ) , IT
FORMATOX, 'SOURCE
484
485
RATE '.12,' STARTS AT ' . F8 . 1 , A3 ,•/..
1 3X,'ENTER MASS RATE AND ENDING TIME ',/,' ?,?')
READ(NI,495,ERR=484) Q(I,M),T(I.M)
495 FORMAT(2F10.0)
500 IF(T(I,M).GT.T(I,M-1)) GO TO 510
504 WRITE(NO,505)
505 FORMATOX,'ENDING TIME MUST BE GREATER THAN STARTING TIME '
1 ' -- REENTER',/,'?')
READ(NI,95.ERR=504) T(I,M)
GO TO 500
510 CONTINUE
520 CONTINUE
GO TO 540
530 WRITE(N0.535) I
535 FORMATOX,'ENTER STEADY-STATE MASS RATE'.I2,/,' ?')
READ(NI,95,ERR=530) 0(1,1)
NR(I) =0
540 CONTINUE
IF( IEDIT.EQ.2.AND.JFLOW.EQ. 1 .AND.TF.'LE . 1 .OE-06) GO TO 720
550 GO .TO (560,3000),IEDIT
C
C COORDINATES OF THE OBSERVATION POINTS
560 WRITE(NO,565) IL
565 FORMATOX,'ENTER XIRST, XLAST, DELTAX (',A2,')' ./.' ?,?,?')
READ(NI,575,ERR = 560) XF(1),XL(1),DEL( 1 )
575 FORMATOF10.0)
DEL(1) = ABS(DEL(1))
IF(DEL(1 ) .LE.1.OE-06) XL(1)=XF(1)
GO TO (580.3000),KNTL
580 WRITE(NO,585) IL
585 FORMATOX,'ENTER YFIRST, YLAST, DELTAY
READ(NI,575,ERR=580) XF(2).XL(2),DEL(2)
DEL(2) = ABS(DEL(2))
C.A2, ')',/.
PL3D211
PL3D212
PL3D213
PL3D214
PL3D215
PL3D216
PL3D217
PL3D218
PL3D219
PL3D220
PL3D221
PL-3D222
PL3D223
PL3D224
PL3D225
PL3D226
PL3D227
PL3D228
PL3D229
PL3D230
PL3D231
PL3D232
PL3D233
PL3D234
PL3D235
PL3D236
PL3D237
PL3D238 .
PL3D239
PL3D240
PL3D241
PL3D242
PL3D243
PL3D244
PL3D245
PL3D246
PL3D247
PL3D248
PL3D249
PL3D250
PL3D251
PL3D252
PL3D253
PL3D254
PL3D255
PL3D256
PL3D257
PL3D258
.PL3D259
PL3D260
PL3D261
PL3D262
PL3D263
PL3D264
PL3D265
PL3D266
PL3D267
PL3D268
PL3D269
PL3D270
PL3D271
PL3D272
PL3D273
PL3D274
PL3D275
PL3D276
PL3D277
PL3D278
PL3D279
PL3D280
B-5
-------
IF(DEL(2).LE.1.OE-06) XL(2)=XF(2)
GO TO (590.3000),KNTL
C
590 WRITE(NO,595) IL
595 FORMATOX.'ENTER ZFIRST, ZLAST, DELTAZ ('.A2,')',/,' ?,?,?')
READ(NI,575,ERR=590) XF(3),XL(3),DELO)
DELO) = ABS(DELO))
600 IF(XFO) .GE.O.O.AND.XFO) .LE.ST.AND.DELO) .LE. 1 .OE-06) GO TO 620
IF(XF( 3) .GE.O.O.AND.XFO) .LE.ST) GO TO 610
604 WRITE(N0.605) ST.IL
605 FORMATOX,'ZFIRST MUST BE GREATER THAN OR EQUAL TO ZERO AND',/.
13X,' LESS THAN OR EQUAL TO SATURATED THICKNESS (', F 10. 4 , A3 ,')','/,
2' -- REENTER'./, ' ?' )
READ(NI.95.ERR=604) XF(3)
• GO TO 600
610 IF(XL(3) .GE.0.0.AND.XLO) .LE.ST) GO TO 630
614 WRITE(NO,615) ST.IL
615 FORMATOX.'ZLAST MUST BE GREATER THAN OR EQUAL TO ZERO AND',/,
13X,' LESS THAN OR EQUAL TO SATURATED THICKNESS ('.F10.4,A3,')',/,
23X,' -- REENTER',/,'?')
READ(NI,95,ERR=614) XLO)
GO TO 630
620 XLO) = XF(3)
630 GO TO (640.3000), IEDIT
C
C PLANE FOR SECTIONING'AQUIFER
640 WRITE(NO,645)
645 FORMATOX, 'ENTER PLANE FOR- SECTIONING AQUIFER (XY.XZ.OR YZ)',
650
660
665
•
670
680
690
'
700
READfNI ,45) ISEC
DO 660 LSEC=1 ,3
IF(ISEC.EQ.KS.ECUSEC) ) GO TO 67O
CONTINUE
WRITE(NO,665)
FORMATOX, ' INVALID SECTION -- REENTER',/,'
GO TO 650
GO TO (680.690,700). LSEC
KHAR1 = KAR3
KHAR2 = KAR1
KHAR3 = KAR2
GO TO 710
KHAR1
KHAR2
KHAR3
KAR2
KAR1
KAR3
GO TO 7 1O
KHAR1
KHAR2
KHARS
KAR1
KAR2
KAR3
710 GO TO (720, 3000), IEDIT
OBSERVATION TIMES
720 IF(JFLOW.EQ.2) GO TO 770
724 WRITE(NO,725) IT
725 FORMATOX. 'ENTER TFIRST. TLAST. DELTAT (', A2 .')'./,' ?,?,
730 READ(NI,575.ERR=724) TF.TL.DELT '
DELT = ABS(DELT)
740 IF(TF. GT.O.O. AND. DELT. LE. 1 .OE-06) GO TO 760
IF(TF.GT.O.O) GO TO 750
744 WRITE(NO,745)
745 FORMATOX. 'TFIRST MUST BE GREATER THAN ZERO -- REENTER',/.'
READ(NI.95.ERR=744) TF
GO TO 740
750 IF(TL. GT.O.O) GO TO 770
754 WRITE(NO,755)
755 FORMATOX, 'TLAST MUST BE GREATER THAN ZERO -- REENTER',/.'
READ(NI ,95,ERR=754) TL
GO TO 750
760 TL = TF
770 GO TO ( 1000, 780), IEDIT
780 IF(JFLOW.EQ.2) WRITE(NO , 785 )
PL3D281
PL3D282
PL3D283
PL3D284
PL3D285
PL3D286
PL3D287
PL3D288
PL3D289
PL3D290
PL3D291
.PL3D292
PL3D293
PL3D294
PL3D295
PL3D296
PL3D297
PL3D298
PL3D299
PL3D300
PL3D301
PL3D302
PL3D303
.PL3D304
PL3D305
PL3D306
PL3D307
PL3D308
PL3D309
PL3D310
PL3D311
PL3D312-
PL3D313
PL3D314
PL3D315
PL3D316
PL3D317
PL3D318
PL3D319
PL3D32O
PL3D321
PL3D322
PL3D323
.PL3D324
PL3D325
PL3D326
PL3D327
PL3D328
PL3D329
PL3D330
PL3D331
PL3D332
PL3D333
PL3D334
PL3D335
PL3D336
PL3D337
PL3D338
PL3D339
PL3D340
PL3D341
PL3D342
PL3D343
PL3D344
PL3D345
PL3D346
PL3D347
PL3D348
PL3D349
PL3D350
B-6
-------
785 FORMAT(3X,'TIME IS NOT A PARAMETER IN STEADY-STATE SOLUTION')
GO TO 3000
LIST PROBLEM DEFINITION
1000 WRITE(NO,1005) NPAGE,(TITLE(I),1=1,30)
1005 FORMAT(1H1,/,3X,'PLUME3D',/,3X,'VERSION 2.02'.
1/.3X.'PAGE ',I3,///.3X,30A2.///)
NPAGE = NPAGE +1
IF(ST.LE.0.9E32) WRITE(NO,1015) IL.ST
1015 FORMAT(1H0.2X,'SATURATED THICKNESS, (',A2,') '.26X.F10.4)
WRITE(NO,1025) IL,IT,V,IL,IT,DX,IL.IT,DY,IL,IT,DZ,P
1025 FORMAT(3X, 'SEEPAGE VELOCITY, (',A2, '/',A2, ' ) ' ,27X,F10.4,/,
13X,'X DISPERSION COEFFICIENT (',A2.'**2/',A2,') ',15X,F10.4,/,
23X,'Y DISPERSION COEFFICIENT (',A2,'-*2/',A2.') ',15X,F10.4./.
33X.'Z DISPERSION COEFFICIENT (',A2,'**2/',A2,') '.15X.F10.4,/,
43X,'POROSITY '.44X.F10.4)
WRITE(NO,1035) RD,IT,DECAY
1035 FORMAT(//,3X,'RETARDATION COEFFICIENT',30X,F10.4,/,
13X,'FIRST ORDER DECAY CONSTANT (1/',A2,')',20X,F10.4)
GO TO (1070,1040), JFLOW
1040 WRITEtNO,1045) IL.IL,IL,IM1.IM2,IMS.IL,IT
1045 FORMAT(//,3X,'STEADY-STATE SOURCE RATES'.//.
13X.'SOURCE',6X,'X'.11X,'Y',11X,'Z',17X,'RATE',/.
25X.'NO',6X,'(',A2,')',8X.'(',A2,')',8X.'(',A2,')',6X,'(',3A2,
3')(CU ',A2,'/',A2,')',/)
DO 1O60 I=1.NS
WRITE(NO, 1055) 'I,XS(I ) ,YS(I),ZS(I).0(1,1)
1055 FORMAT(5X.I2.F10.2,2X,F10.2.2X,F1O.2.6X.F16.4)
1060 CONTINUE
GO TO 11 10
1070 WRITEtNO,1075) IM1,IM2,IM3,IL,IT,IT,IL,IL,IL .
1075 FORMAT(//,3X. 'SOURCE/RATE SCHEDULE (•' .3A2, ' )(CU '.A2,'/'.A2.
1')'//, 15X, 'SOURCE' , 13X. 'RATE' ,4X. 'MASS',8X, 'TIME (',A2.')',
Y (',A2,')
,ZS(II
Z ( ' ,A2, ' )
NO',5X.'RATE',
2/.3X,'NO X (',A2,')
35X . ' START.' . 6X , ' END ' , / )
DO 1100 1=1.NS
WRITEtNO,1085) I,XS(I),YS(I
1085 FORMAT(/,3X.I2.3F9.2)
NRT = NR(I)
DO 1100 J=1,NRT
M = J + 1
WRITE(NO,1095) J,Q(I,M),T(I,M-1),T(I,M)
1095 FORMAT(34X. I2..F12.2.2F9.2)
1100 CONTINUE
1110 WRITE(NO,1115) IL,IT,XF(1),XL(1),DEL(1),XF(2),XL(2),DEH2)
1 XF(3).XL(3),DEL(3)
'OBSERVATION POINTS (',A2
'.F10.2.5X,'XLAST
',F10.2,5X,'YLAST
1,F10.2,5X,'ZLAST
IF(JFLOW.EO.1) WRITE(NO,1125) TF.TL.DELT
1125 FORMAT(/,5X,'TFIRST =',F10.2,5X,'TLAST =
WRITE(NO,1135) KSEC(LSEC)
1135 FORMAT(//,3X,'AQUIFER SECTIONED IN ',A2,
GO TO 3000
1115 FORMAT(//,3X,
15X,'XFIRST
25X.'YFIRST
35X.'ZFIRST
1 .F10.2.5X,
' .F10.2.5X.
.F10.2.5X.
AND TIMES (',A2,
DELX ='.F10.4,/,
DELY
'DELZ
F10.4./,
F10.4)
.F10.2.5X,'DELT
PLANE')
«**«« SECTION II -- NUMERICAL EVALUATION OF CONCENTRATION AT
SPECIFIED GRID COORDINATES
NUMBER OF OBSERVATION POINTS IN EACH COORDINATE DIRECTION
2000 CONTINUE
DO 2020 L=1,3
NP(L) = 1
DEL(L) * ABS(DELU))
IF(DEL(L).LE.1.OE-03) GO TO 2020
DIF = XL(L) - XF(L)
PL3D351
PL3D352
PL3D353
PL3D354
PL3D355
PL3D356
PL3D357
PL3D358
PL3D359
'PL3D360
PL3D361
PL3D362
PL3D363
PL3D364
PL3D365
PL3D366
PL3D367
PL3D36B
PL3D369
PL3D370
PL3D371
PL3D372
PL3D373
PL3D374
PL3D375
PL3D376
PL3D.377
PL3D378
PL3D379
PL3D380
PL3D381
PL3D382
PL3D383
PL3D3S4
PL3D385
PL3D386
PL3D387
PL3D388
PL3D389
PL3D39C
PL3D391
PL3D392
PL3D393
PL3D394
PL3D395
PL3D396
PL3D397
PL3D398
- PL3D399
PL3D400
PL3D401
PL3D402
F10.4) PL3D403
PL3D404
PL3D405
PL3D40S
PL3D407
PL3D408
PL3D409
PL3D410
PL3D411
PL3D412
PL3D413
PL3D414
PL3D415
PL3D416
PL3D417
PL3D418
PL3D419
PL3D420
3-7
-------
IF(ABS(DIF).LE.1.OE-03) GO TO 2020
IF(DIF.LE.O.O) DEL(L)=-DEL(L)
NPTS = ABS(DIF/DEL(D)
REM = DIP - DEL(L)*FLOAT(NPTS)
NPTS = NPTS + 1
NP(L) = NPTS
IF(ABS(REM).LT.1.OE-03) GO TO 2020
NP(L) = NP(L) + 1
2020 CONTINUE
GO TO (2040,2060,2080),LSEC
2040 MAXSC = NP(3)
MAXRW = NP(2)
MAXCL = NP(1)
GO TO 2100
2060 MAXSC = NP(2)
MAXRW = NP(3)
MAXCL = NP(1)
GO TO 2100
2080 MAXSC = NP(1)
MAXRW = NP(3)
MAXCL = NP(2)
2100 CONTINUE
TIME COORDINATES
NTIME = 1
IF(DELT.LE.1.OE-O6) GO TO 2110
NTIME = ABS(TL-TF)/DELT + 1.0
IF(TF.GT.TL) DELT=-DELT
2110 TSOL = TF
MTIME = NTIME
2120
2140
2160
2180
2200
2220
2240
DAMK = DX*DECAY-RD/(V-V)
ALPHA=SORT(1.0+4.O'DAMK)
PE=V/DX
BETA = DX/DY
GAMMA = DX/DZ
LAMBDA = 1 .O/(25. 132741*P«SORT(DY«DZ) '.
DO 2660 NT=1,NTIME
NSEC =' 1
LPRT = 1
LP = 1
NCFLG = 1
NROW1 = 1
NROW2 = MAXROW
IF(NROW2.GT.MAXRW) NROW2=MAXRW
DO 2580 NROW=NROW1.NROW2
GO TO (2180,2220,2200),NCFLG
NCOL1 = 1
NCOL2 = MAXCOL
IF(NCOL2.GT.MAXCL) NCOL2=MAXCL
NCOL = MAXCOL
IF(NCOL2.EQ.MAXCL) NCOL=NCOL2-NCOL1+1
GO TO (2240,2260,2280), LSEC
1X1 = NCOL1
1X2 = NCOL2
JY1 = NROW
JY2 = NROW
KZ1 = NSEC
KZ2 = NSEC
GO TO 2290
2260
1X1
1X2
JY1
JY2
KZ1
KZ2
NCOL1
NCOL2
NSEC
NSEC
NROW
NROW
PL3D421
PL3D422
PL3D423
PL3D424
PL3D425
PL3D426
PL3D427
PL3D428
PL3D429
PL3D430
.PL3D431
PL3D432
PL3D433
PL3D434
PL3D435
PL3D436
PL3D437
PL3D438
PL3D439 .
PL3D440
PL3D441
PL3D442
PL3D443
PL3D444
PL3D445
PL3D446
PL3D447
PL3D448
PL3D449
PL3D450
PL3D451
PL3D452
PL3D453
PL3D454
PL3D455
PL3D456
PL3D457
PL3D45S
PL3D459
PL3D460
PL3D461
PL3D462
PL3D463
PL3D464
PL3D465
PL3D466
PL3D467
PL3D468
PL3D469
PL3D470
PL3D471
PL3D472
PL3D473
PL3D474
PL3D475
PL3D476
PL3D477
PL3D478
PL3D479
PL3D480
PL30481
PL3D482
PL3D483
PL3D484
PL3D485
PL3D486
PL3D487
PL3D488
PL3D489
PL3D490
B-8
-------
2280
2290
2300
GO TO 2290
1X1 = MSEC
1X2 = NSEC
JY1 = NCOL1
JY2 = NCOL2
KZ1 = NROW
KZ2 = NROW
CONTINUE
DO 2300 L =
CON(L) *
CONTINUE
1.MAXCOL
0.0
DO 2440 N=1,NS
D( 1) = ST - ZS(N)
IF(ST.GE.0.9E32) D(1)=0.0
D(2) = 2S(N)
D(3) = D(1)
COEF =1.O
IF(D(1).LT.1.OE-03.0R.D(2).LT.1.0E-03) COEF=2.0
DO 2440 1 = 1X1 .1X2
X = XF(1) + FLOAT(I-1)*DEL(1)
IF(I.EQ.NP(1)) X = XL(1)
XXS = X - XS(N) .
PEX = PE*XXS
DO 244O J=JY1.JY2
Y » XF(2) + FLOAT(d-1)«DEL(2)
IF(d.EO.NP(2)) Y=XL(2)
YYS = Y - YS(N)
PEY = PE*YYS
DO 2430 K=KZ1,KZ2
L = I-IX1 + J-JY1 ••- K-KZ1 + 1
.IF(CON(L).LT.0.0) GO TO 2430
Z = XF(3) + FLOAT(K-1)»DEL(3)
• IF(K.EO.NP(3)) Z=XL(3)
ZM = Z - ZS(N)
' IF(ABS(XXS).LT.1.0.AND.ABS(YYS).LT.1.0.AND.
1 AES(ZM).LT.1.0) GO TO 2330
PEZ = PE-ZM
' CALL SOL3D(C.PEX,PEY,PEZ,TSOL,N.NR(N))
CXYZT = 'COEF*C
IF(IMAGE.EQ.1) GO TO 2325
DO 2320 LM=1,2 •
ZM = ( (-1 .0)«*(LM+1 ))*ZM
IF(D(LM).LT.1.OE-03) GO TO 2320
ZIMAGE = 2.0«D(LM) - ZM
PEZ = PE'ZIMAGE
CALL SOL3D(C,PEX,PEY,PEZ,TSOL,N,NR(N))
CXYZT = CXYZT + COEF*C
DO 2310 IM = 1.IMAGE
ZIMAGE = (2.0«D(LM)+ZM) + 2.0»FLOAT(IM)-D(LM+1)
1 + FLOAT(2«IM-2)*D(LM)
PEZ = PE'ZIMAGE
CALL SOL3D(C,PEX,PEY,PEZ,TSOL,N,NR(N))
IF(C.LT.1.OE-06) GO TO 2312
CXYZT = CXYZT + COEF'C
2310 CONTINUE
2312 CONTINUE
DO 2314 IM=1,IMAGE
ZIMAGE = (2.0*D(LM)-ZM) + 2.0*FLOAT(IM)-D(LM+1)
1 + FLOAT(2»IM)*D(LM)
PEZ = PE*ZIMAGE
CALL SOL3D(C,PEX,PEY,PEZ.TSOL,N,NR(N))
IF(C.LT.1.OE-06) GO TO 2320
CXYZT = CXYZT + COEF*C
2314 CONTINUE
WRITE(NO,2315) MAXIMG.X.Y.Z
2315 FORMAT(3X,'»»*»«* WARNING -- SOLUTION DID NOT',
1 ' CONVERGE USING' ,/,9X.12, ' IMAGE WELLS AT X =' ,
PL3D491
PL3D492
PL3D493
PL3D494
PL3D495
PL3D496
PL3D497
PL3D498
PL3D499
PL3D500
PL3D501
PL3D502
PL3D503
PL3D504
PL3D505
PL3D506
PL3D507
PL3D508
PL3D509
PL3D510
PL3D511
PL3D512
PL3D513
PL3D514
PL3D515
PL3D516
PL3D517
PL3D518
PL3D519
PL3D520
PL3D521.
PL3D522
PL3D523
PL3D524
PL3D525
PL3D526
PL3D527
PL3D528
PL3D529
PL3D530
PL3D531
PL3D532
PL3D533
PL3D534
PL3D535
PL3D536
PL3D53T
PL3D538
PL3D539
PL3D540
PL3D541
PL3D542
PL3D543
PL3D544
PL3D545
PL3D546
PL3D547
PL3D548
PL3D549
PL3D550.
PL3D551
PL3D552
PL3D553
PL3D554
PL3D555
PL3D556
PL3D557
PL3D558
PL3D559
PL3D560
B-9
-------
2 MO.4,' Y =',F10.4,' Z =',F10.4)
2320 CONTINUE
2325 CON(L) = CON(L) + CXYZT
GO TO- 2340
2330 CON(L) = -9.9999
2340 GO TO (2360.2380,2400), LSEC
2360 SEC = Z
ROW = Y
COL(L) = X .
GO TO 2420
2380 SEC = Y
ROW'= Z
COL(L) = X
GO TO 2420
2400 SEC = X
ROW = Z
COL(L) = Y
2420 CONTINUE
2430 CONTINUE
2440 CONTINUE
: PRINT CONCENTRATION DISTRIBUTION
GO TO (2460.2560). LPRT
2460 WRITE(NO,1005) NPAGE, (TITLE(I),I=1.30)
NPAGE = NPAGE + 1
IF(UFLOW.EQ.2) GO TO 2500
WRITE(NO,2465) TSOL,IT.IM1,IM2.IMS.KHAR1,SEC.IL.
1(LBL(LP,L),L=1,6).KHAR2,IL .
2465 FORMAT(13X,'CONCENTRATION DISTRIBUTION AT '.F10.2,
11X.A2,' (',3A2.') ' ,//. 13X.A1, ' =' .F10.2, 1X,A2,3X,6A2.//,
2' «',/,' * ',A1,'(',A2,')')
GO TO 2520
2500 WRITE.(NO,2505) IM1 , IM2, IMS. KHAR 1 . SEC ,'IL , ( LBL( LP . L) . L = .1 , 6 ) ,
1KHAR2.IL
2505 FORMAT('13X, 'CONCENTRATION DISTRIBUTION AT STEADY STATE',
1' ('.3A2,') './/13X.A1.' =' ,F10.2, 1X,A2,3X,6A2,//.
2' »''./.' « ' , A1 , ' ( ' . A2, ' ) ' )
2520 CONTINUE
WRITE(NO,2525) (COL(L).L=1.NCOL)
2525 FORMAT;' •',4x,7Fio.2)
WRITE(NO,2545') KHAR3.IL
2545 FORMAT(1X.A1 ,"(' ,A2. ' ) *',/,9X,'«")
2560 WRITE(NO,2565) ROW,(CON(L),L=1,NCOL)
2565 FORMAT(2X,F8.2,7F10.4)
LPRT = 2
2580 CONTINUE
IF(NROW2.EO.MAXRW) GO TO 2600
NROW1 = NROW1 + MAXROW
NROW2 = NROW2 + MAXROW
LPRT = 1
LP = 2
NCFLG = 2
GO TO 2160
2600 I.F(NCOL2.EO.MAXCL) GO TO 2620
NCOL1 = NCOL1 + MAXCOL
NCOL2 = NCOL2 + MAXCOL
LPRT = 1
LP = 2
NCFLG = 3
GO TO 2140
2S20 IF(NSEC.EQ.MAXSC) GO TO 2640
NSEC = NSEC + 1
GO TO 2120
2640 CONTINUE
TSOL = TSOL + DELT
IF(NT.EQ.MTIME) TSOL=TL
2660 CONTINUE
PL3D561
PL3D562
PL3D563
PL3D564
PL3D565
PL3D566
PL3D567
PL3D568
PL3D569
PL3D570
PL3D571
PL3D572
PL3D573
PL3D574
PL3D575
PL3D576
PL3D577
PL3D578
PL3D579
PL3D580
PL3D581
PL3D582
PL3D583-
PL3D584
PL3D585
PL3D586
PL3D587
PL3D588
PL3D589
PL3D590
P.L3D591
PL3D592
PL3D593
PL3D594
PL3D595
PL3D596
PL3D597
PL3D598
PL3D599
PL3D600
PL3D601
PL3D602
PL3D603
PL3D604
PL3D605
PL3D606
PL3D607
PL3D608
PL3D609
PL3D610
PL3D611
PL3D612
PL3D613
PL3D614
PL3D615
PL3D616
PL3D617
PL3D618
PL3D619
PL3D620
PL3D621
PL3D622
PL3D623
PL3D624
PL3D625
PL3D626
PL3D627
PL3D628
PL3D629
PL3D630
B-10
-------
C »«««« SECTION III -- PROBLEM REDEFINITION AND CONTROL OF EXECUTION PL3D631
C PL3D632
C PL3D633
3000 CONTINUE PL3D634
IF(IEDIT.E0.2) GO TO 3010 PL3D635
WRITE(NO,3705) PL3D636
IEDIT = 2 . PL3D637
3010 KNTL =2. PL3D638
WRITE(NO,3015) • PL3D639
3015 FORMAT*//,3X,'ENTER NEXT COMMAND',/,' ?') PL3D640
3020 READ(NI,3025) NEXT ' PL3D641
3025 FORMATU2) PL3D642
C . ' PL3D643
DO 3030 1=1,21 PL3D644
IF(NEXT.EO.IC(I)) GO TO 3040 PL3D645
3030 CONTINUE PL3D646
WRITE(NO,3035) . PL3D647
3O35 FORMAT(3X,'ERROR IN LAST COMMAND -- REENTER',/,' ?') PL3D648
GO TO 302O PL3D649
3040 GO TO (10,BO,120,160,190,220,250,280.310.320,560,580,590,640. PL3D650
1 720,1000,2000,3050,3060,3700,4000),! PL3D651
C PL3D652
C NEW SET OF X, Y, AND Z OBSERVATIONS PL3D653
3050 KNTL = 1 • PL3D65d
GO TO 560. PL3D655
C ' ' . . PL3D656
C PL3D657
C . ' PL3D658
C NEW SOURCE/RATE SCHEDULE ' . PL3D659
3060 WRITE(NO,3065) . PL3D660
3065 FORMAT(3X,'ADD(A).DELETE(D).MODIFY(M) A SOURCE OR RETURN(R)' PL3D661
1 ' TO EDIT ?' ). . - PL3D662
3O7O READtNI ,3075)- ISK PL3DS63
3075 FORMAT!A1) . ' • PL3D664
DO 3080 1=1,4 ' PL3D6S5
IFdSK'.EQ.ISd) ) GO TO 3090 • PL3D666
3080 CONTINUE ' PL3DG67
WRITE(NO,3085) . PL3D668
3085 FORMAT!3X,'ERROR IN SELECTION -- REENTER ?') . PL3D669
GO TO 3070 PL3D670
3090 GO TO (3000,3100.3450,3490),! PL3D671
C . . PL3D672
C MODIFY SOURCE PL3D673
C PL3D674
3100 WRITE(NO,3105) NS ' • PL3D675
3105 FORMAT(3X,I2,' SOURCES IN CURRENT SCHEDULE',/. PL3D676
13X,'ENTER SOURCE TO MODIFY',/.' ?') PL3D677
READ(NI,465,ERR=3100) 'FDUM PL3D678
JS=FDUM • PL3D679
IF(JS.GT.O.AND.JS.LE.NS) GO TO 322O PL3D680
WRITE(NO,3215) JS . PL3D681
3215 FORMAT(3X,'SOURCE',14,' NOT IN SCHEDULE') PL3D682
GO TO 3060 PL3D683
3220 GO TO (3230.3260),JFLOW PL3D684
3230 WRITE(NO,3235) JS,XS(JS),IL,YS(JS),IL,ZS(JS),IL,IT, PL3D685
1IM1,IM2,IM3,IL,IT PL3D686
.12,': X =' .F8.2.A3, ' . Y =',F8.2,A3, PL3D687
•' ,F8.2,A3,//.3X, 'RATE',7X. 'MASS RATE' , 14X. 'TIME (', PL3D688
2A2, ' )' ./,4X, 'NO' ,3X.'(' .3A2.')(CU ',A2. ' /' ,A2.')', PL3D689
38X.'START'.7X.'END',/) PL3D690
NRT = NR(JS) PL3D691
DO 3250 J=1,NRT PL3D692
M = J + 1 PL3D693
WRITE(NO,3245) J,0(JS,M),T(JS,M-1),T(JS,M) PL3D694
3245 FORMAT(4X.I2,5X.F14.2,7X,F8.3,3X,F8.2) PL3D695
3250 CONTINUE PL3D696
GO TO 3270 PL3D697
3260 WRITE(NO,3265) JS,XS(JS),IL,YS(JS),IL.ZS(JS),IL.0(JS,1). PL3D698
1IM1,IM2,IM3.IL,IT PL3D699
3265 FORMAT(3X,'SOURCE ',12,': X =',F8.2,A3,'. Y =',F8.2,A3, PL3D700
3235 FORMAT(3X,'SOURCE
1 ' . Z
B-ll
-------
1'. Z=',F8.2.A3,/,3X,'STEADY-STATE MASS RATE =',M6.4,
2' (',3A3,')(CU ',A2,'/',A2,')',/)
3270 WRITE(NO,3275)
327S FORMAT(3X.'CHANGE COORDINATES (Y/W)?')
READ(NI,3075) OC
IF(JC.NE.IY) GO TO 3290
3276 WRITE(NO,415) JS.IL
READ(NI,425,ERR=3276) XS(JS).YS(JS),ZS(JS)
3280 IF(ZS(JS).GE..O.O.AND.ZS(JS).LE.ST) GO TO 3290
3284 WRITE(NO,435) ST.IL
READ(NI,95,ERR=3284) ZS(JS)
GO TO 3280
3290 GO TO (3300,3430),JFLOW
C
C TRANSIENT SOURCES
3300 WRITE(NO,3305) JS
3305 FORMAT(3X,'MODIFY RATE SCHEDULE FOR SOURCE',13,' (Y/N) ?')
READ(NI,3075) JY
IF(JY.NE.IY) GO TO 3060 '•
3310 WRITE(NO,3315)
3315 FORMAT(3X,'ENTER RATE TO BE CHANGED',/,
13X,'(ENTER 0 TO CHANGE ALL RATES)',/.' ?')
READ(NI,465,ERR=331O) FDUM
JR=FDUM
IF(JR.LE.O) GO TO 335O
IF(JR,LE.NR(JS)) GO TO 3330
WRITE(NO,3325)-JR
3325 FORMAT(3X,'RATE ',12,' NOT IN CURRENT SCHEDULE')
GO TO 3300 . • ' •
3330 WRITE(NO,3335) JS,JR,T(JS.JR),T(JS.JR+1 ) , IT
3335 FORMAT(3X,'SOURCE ',12,'. RATE ',12,' STARTS AT',F8.2,
1' .AND ENDS AT',F8.2,A3./.3X,'ENTER NEW MASS RATE'./.' ?'
M = JR + 1
READINI.3345.ERR=3330) Q(JS.M)
3345 FORMAT( F'10.0)
GO TO 3300
C
3350 NRT = NR(JS)
DO 33SO d=1.NRT
M = J + 1
O(JS.M) = 0.0
T(JS,M) = 0.0
3360 CONTINUE
3370 WRITE(NO,455) JS.MAXRT
3380 READ(NI,465,ERR=3370) FDUM
NR(dS)=FDUM
IF(NR(JS).GT.O.AND.NR(JS).LE.MAXRT) GO TO 3390
WRITE(NO,475) MAXRT
GO TO 338O
3390 CONTINUE
NRT = NR(dS)
DO 3420 J=1,NRT
M = J + 1
3394 WRITE(NO,485) JS, J.T(JS.M- 1 ), IT
READ(NI,495.ERR=3394) 0(JS.M),T(JS,M)
3400 IF(T(JS.M).GT.T(JS.M-1)) GO TO 3410
3404 WRITE(NO,505)
READ(NI,95,ERR=3404) T(JS,M)
GO TO 3400
3410 CONTINUE
3420 CONTINUE
GO TO 3060
C
C STEADY-STATE SOURCES
3430 WRITE(NO,3435) JS
3435 FORMAT(3X,'CHANGE STEADY-STATE RATE FOR SOURCE
READ(NI,3075) JC
IF(JC.NE.IY) GO TO 3060
3444 WRITE(NO,3445) JS
3445 FORMAT(3X,'ENTER NEW STEADY-STATE MASS RATE FOR SOURCE
,12,' (Y/N) ?')
PL3D701
PL3D702
PL3D703
PL3D704
PL3D705
PL3D706
PL3D707
PL3D708
PL3D709
PL3D710
PL3D71L
PL3D712
PL3D713
PL3D714
PL3D715
PL3D716
PL3D717
PL3D718
PL3D719
PL3D720
PL3D721
PL3D722
PL3D723
PL3D724
PL3D725
PL3D726
PL3D727
PL3D728
PL3D729
PL3D730
PL3D731
PL3D732
PL3D733
PL3D734
PL3D735
PL3D736
PL3D737
PL3D738
PL3D739
PL3D740
PL3D741
PL3D742
PL3D743
PL3D744
PL3D745
PL3D746
PL3D747
PL3D748
PL3D749
PL3D750
PL3D751
PL3D752
PL3D753
PL3D754
PL3D755
PL3D756
PL3D757
PL3D758
PL3D759
PL3D760
PL3D761
PL3D762
PL3D763
PL3D764
PL3D765
PL3D766
PL3D767
PL3D768
PL3D769
PL3D770
B-12
-------
1- ?')
READ(NI,3345,ERR=3444) Q(JS,1)
GO TO 3060
C
C ADD A NEW SOURCE
C
3450 NS = NS + 1
JS = NS
3454 WRITE (NO, 41,5) J.S.IL
READ(NI,425,ERR=3454) XS(dS),YS(JS),ZS(JS)
3460 IF(2S(JS).GE.O.O.AND.ZS(JS).LE.ST) GO TO 3470
3464 WRITE(NO,435) ST.IL
READ(NI,95,ERR=3464) ZS(JS)
GO TO 3460
3470 GO TO (3370,3480).JFLOW
C
C STEADY-STATE SOURCES
3480 WRITE(NO,3485) JS
3485 FORMATOX, 'ENTER STEADY-STATE MASS RATE FOR SOURCE
1/,' ?')
READ(NI,3345.ERR=3480) 0(JS.D
NR(JS) = 0
GO TO- 3060
3490
3495
3500
3505
3515
3520
3530
3535
,6X, 'X ( ' , A2 , ' ) ' , 3X . ' Y ( ' , A2 , ' ) ' . 3X .
I.XS(I), YS(I ) ,ZS(I )
. 3X , F8 . 2 )
DELETE A SOURCE
IF(NS.GT.I) GO TO 35OO
WRITE(NO,3495)
FORMATOX, 'ONLY ONE SOURCE IN SCHEDULE -- CAN NOT DELETE'
GO TO 3060
WRITE (NO, 3505) IL.IL.IL
FORMAT(3X, 'SOURCE'
1'Z C ,A2,. ')',/)
DO 352O 1 = 1 ,NS .
WRITE(N0.3515).
FORMAT(5X, 12 , 3X , F8 . 2 , 3X . F8 . 2 .
CONTINUE
WRITE(NO,3535)
FORMATOX, 'ENTER SOURCE TO DELETE'./,
13X, '(ENTER 0 TO CANCEL)', /,' ?')
READ(NI.465,ERR=3530) FDUM
JS=FDUM
IF(JS.LE.O) GO TO .3060
IF(dS.LE.NS) GO TO 3550
WRITE(NO,3545) JS
FORMATOX, 'SOURCE ',12.' NOT IN CURRENT SCHEDULE')
GO TO 3530
WRITE(NO,3555) JS
FORMATOX, 'DELETE SOURCE
READ(NI,3075) JC
IF(JC.NE.IY) GO TO 3530
NSD = NS - 1
GO TO (3560.3590), JFLOW
: TRANSIENT SOURCES
3560 IF(JS.EO.NS) GO TO 3575
DO 3570 J=JS.NSD
3545
3550
3555
(Y/N)?')
3570
3575
XS(J) = XS(J-M)
YS(J) = YS(J+1)
ZS(J) = ZS(J+1)
NR(J) = NR(J+1 )
NRT = NR(d)
DO 3570 K=1 .NRT
M = K + 1
Q( J,M) = Q(J+1 ,M)
T( J.M) = T(J+1 ,M)
CONTINUE
NRT = NR(NS)
DO 3580 K=1 ,NRT
M = K + 1
PL3D771
PL3D772
PL3D773
PL3D774
PL3D775
PL3D776
PL3D777
PL3D778
PL3D779
PL3D780
PL3D781
PL3D782
PL3D783
PL3D784
PL3D785
PL3D786
PL3D787
PL3D788
PL3D789
PL3D790
PL3D791
PL3D792
PL3D793
PL3D794
PL3D795
PL3D796
PL3D797
PL3D798
PL3D799
PL3D800
PL3D801
PL3D802
PL3D803
PL3D804
PL3D805
PL3D806
PL3D807
PL3D808
PL3D809
PL3D810
PL3D811
PL3D812
PL3D813
PL3D814
PL3D815
PL3D816
PL3D817
PL3D818
PL3D819
PL3D820
PL3D821
PL3D822
PL3D823
PL3D824
PL3D825
PL3D826
PL3D827
PL3D828
PL3D829
PL3D830
PL3D831
PL3D832
PL3D833
PL3D834
PL3D835
PL3D836
PL3D837
PL3D838
PL3D839
PL3D840
B-13
-------
0(NS,M) = 0.0
T(NS,M) = 0.0
3580 CONTINUE
NR(NS) « 0
NS = NSD
GO TO 3060
C
C STEADY-STATE SOURCES
3590 IF(OS.EO,NS) GO TO 3605
DO 3600 d=JS,NSD
0(J,1) " Q(J+1,1)
XS(J) = XS(d+1)
YS(d) = YS(J+1)
ZS(d) = ZS(J+1)
3600 CONTINUE
3605 0(NS. 1.) = 0.0
NS = NSD
GO TO 3060
C
C
C MENU OF EDIT COMMANDS FOR PLUME3D VERSION 2.O2
3700 WRITE(NO.3705)
3705 FORMAT(1H1,/,3X,'MENU OF EDIT COMMANDS',//.
13X
23X
33X
43X
53X
53X
63X
73X
83X
93X
13X, '
GO TO 3000
SATURATED THICKNESS ST
POROSITY PO
SEEPAGE VELOCITY VX
RETARDATION COEFFICIENT RD
X DISPERSION COEFFICIENT DX
Y DISPERSION COEFFICIENT DY
Z DISPERSION COEFFICIENT DZ
DECAY CONSTANT DE
SOURCE/RATE SCHEDULE RT
CHANGE SOLUTION/SOURCES CS
OBSERVATION POINTS OB',/
X COORDINATES XC',/
Y -COORDINATES . YC'./
Z COORDINATES ZC' , /
OBSERVATION TIMES ' TC' ./
AQUIFER SECTIONING AS',/
NEW PROBLEM NP',/
MENU OF COMMANDS MU'./
LIST INPUT' DATA LI ' ,/
RUN CALCULATIONS RN',,
DONE DM')
4000 STOP
END
PL3D841
PL3D842
PL3D843
PL3D844
PL3D845
PL3D846
PL3D847
PL3D848
PL3D849
PL3D850
PL3D851
PL3D852
PL3D853
PL3D854
PL3D855
PL3D856
PL3D857
PL3D858
PL3D859
PL3D860
PL3D861
PL3D862
PL3D863
PL3D864
PL3D865
PL3D866
PL3D867
PL3D868
PL3D869
PL3D870
PL3D871
PL3D872
PL3D873
PL3D874
PL3D875
PL3D876
PL3D877
PL3D878
B-14
-------
FUNCTION ERFCLG(Z)
RATIONAL APPROXIMATION OF THE COMPLIMENTARY ERROR FUNCTION
SEE SECTION 7.1 OF ABRAMOWITZ AND STEGUN (1966)
THE FOLLOWING IDENTITIES ARE USED TO HANDLE NEGATIVE ARGUMENTS
ERFC( Z) = 1 - ERF(Z)
ERF (-Z) = -ERF (Z)
REAL'S COEFLG.DERFC.DI,FX,TERMI,TERMO.SUM,X
COMMON/IO/NI,NO
X = ABS(Z)
IF (X.GT.3.0DOO) GO TO 50
FOR X<3 A RATIONAL APPROXIMATION OF THE COMPLIMENTARY ERROR
FUNCTION IS USED.
DERFC = 1.ODOO/((1.ODOO + 7.05230784D-02*X + 4.22820123D-02"(X"*2)
1 + 9.2705272D-03*(X*»3) + 1.520143D-04*(X**4 )
2 + 2.76572D-04«(X**5) + 4.3OS38D-05*(X««6))«-16)
ERFCLG = DLOG(DERFC)
GO TO'100
FOR X>3 AN ASYMPTOTCI EXPANSION OF THE COMPLIMENTARY ERROR
FUNCTION IS USED.
50 COEFLG = X-X + DLOG(X) + 0.57236494DOO
FX « 2.0DOO*X"X
SUM = 1.ODOO
TERMO = 1.ODOO
DO 60 1=2,50
DI = I
TERMI = -TERMO»(2.ODOO'DI - 3.0DOO)/FX
IF(DABS(TERMI).GT.DABS(TERMO)) GO TO 70
SUM = SUM + TERMI
TEST = TERMI/SUM.
IF(ABSITEST).LT.1.OE-16) GO TO 70
TERMO = TERMI . .
60 CONTINUE
WRITE(N0.65)
65 FORMAT(6X,'*•* WARNING -- ASYMPTOTIC EXPANSION FOR ERFC DID NOT',
1' CONVERGE WITH 50 TERMS IN THE SUMMATION')
70 SUM = DLOG(SUM) - COEFLG
ERFCLG = SUM
100 CONTINUE
FOR Z<0, ERFC(rZ) = 2-ERFC(Z)
IF(Z.LT.O.O) GO TO 200
RETURN
200 ERFC = 2.0 - EXP(ERFCLG)
ERFCLG = ALOG(ERFC)
RETURN
END
EFLG001
EFLG002
EFLG003
EFLG004
EFLG005
EFLG006
EFLG007
EFLG008
'EFLG009
EFLG010
EFLG011
EFLG012
EFLG013
EFLG014
EFLG015
EFLG016
EFLG017
EFLG018
EFLGO19
EFLG020
EFLG021
EFLG022
EFLG023
EFLG024
EFLG025
EFLG026
EFLG027
EFLG028
EFLGO29
EFLG030
EFLG031
EFLG032
EFLG033
EFLG034
EFLG035
EFLG036
EFLG037
EFLG038
EFLG039
EFLG040
EFLG041
EFLG042
EFLG043
EFLG044
EFLG045
EFLGO46
EFLG047
EFLG048
EFLG049
EFLG05O
EFLG051
B-15
-------
c
c
SUBROUTINE SOL3D( C , PEX , PEY , PE2 , TSOL , N , NR )
NUMERICAL EVALUATION OF ANALYTICAL SOLUTION
REAL LAMBDA
COMMON/RATE/Q(10, 12), T( 10, 12)
COMMON/PHYPRO/ALPHA , BETA , GAMMA , DX , LAMBDA , PE , RD , V
PEXYZ « SORT(PEX*PEX + ,BETA*(PEY*PEY)
R = PEXYZ/PE
PEL = PEXYZ«ALPHA
PEX12 = PEX/2.0
PEL12 = PEL/2.0
COEF = LAMBDA/R
MT = NR + 1 .
IF(MT.GT.1) GO TO 10
STEADY-STATE SOLUTION
S = 0(N,1)
IF(CLG.LT.-72.0) CLG = -72.0
IF(CLG.GT. 72.0) CLG * 72.0
C = 2.0*EXP(CLG)
GO TO 50
TRANSIENT SOLUTION
10 C = 0.0
IF(f(N.MT).LT.TSOL) MT=MT+1
DO 40 K=2,MT
IF(T(N,K-1 ) .GT..TSOL) GO TO 50
S = O(N.K) - Q(N,K-1)
SGN. =' 1 .0
IF(S.LT.O.O) SGN =-1.0
SA = ABS(S)
IFCSA.LT.1.OE-3) GO.TO 40
• COEFLG = ALOG(COEF-SA)
PINJ = V«V»(TSOL-T(N,K-1))/(DX»RD)
PINJL = PINJ-ALPHA»ALPHA
TAU = PINJ/(PEXYZ*PEXYZ)
Z1 = 0.'5/SORT(TAU)
22 = O.5«SORT(PINJL)
Z = Z1 + Z2
T1 = ERFCLG(Z)
C1LG = T1 + PEL12 + PEX12 + COEFLG
Z = Z1 - Z2
T2 = .ERFCLG(Z)
C2LG = T2 - PEL12 + PEX12 + COEFLG
IF(C1LG.GT. 72.0) C1LG= 72.0
IF(C1LG.LT.-72.0) C1LG=-72.0
IF(C2LG.GT. 72.0) C2LG= 72.0
IF(C2LG.LT.-72.0) C2LG=-72.0
CK = EXP(C1LG) + EXP(C2LG)
C = C + SGN'CK
40 CONTINUE
50 RETURN
END
GAMMA*( PEZ*PEZ) )
SOL3001
SOL3002
SOL3003
SOL3004
SOL3005
SOL3006
SOL3007
SOL3008
SOL3009
SOL3010
SOL3011
SOL3012
SOL3013
SOL3014
SOL3015
SOL3016
SOL3017
SOL3018
SOL3019
SOL3020
SOL3021
SOL3022
SOL3023
SOL3024
SOL3025
SOL3026
SOL3027
SOL3028
SOL3029
SOL3030
SOL3031
SOL3032
SOL3033
SOL3034
SOL3035
SOL3036
SOL3037
SOL3038
SOL3039
SOL3040
SOL3041
SOL3042
SOL3043
SOL3044
SOL3045
.SOL3046
SOL3047
SOL3048
SOL3049
SOL3050
SOL3051
SOL3052
SOL3053
SOL3054
SOL3055
B-16
-------
APPENDIX C
Solute Transport in Uniform Ground-Water Flow
Mathematical Development
C-l
-------
SOLUTE TRANSPORT IN UNIFORM GROUND-WATER FLOW
MATHEMATICAL DEVELOPMENT
Analytical solute-transport models are highly idealized mathematical
approximations of complex physical phenomena. However, if applied properly
this class of models represents a useful tool for analyzing the transport and
fate of substances under both field and laboratory conditions.
This mathematical derivation of a governing differential equation for
solute transport in uniform ground-water flow is developed in an effort to
show the relationship between physical processes in the aquifer and terms in
the mathematical model. In this context, the "mathematical manipulations" of
the derivation are not as important as the simplifying assumptions which are
required to obtain a mathematical description of the problem.
The resulting differential equation serves as the basis of many
analytical so.l ute-transport models. Solutions of the differential equation
which would provide the actual concentration distributions in time and space
are not addressed in this document. However, the simplifying assumptions
which are incorporated in the formulation of the mathematical model must be
considered in both the application of the model and in the interpretation of
modeling results, regardless of how the model is solved.
C-2
-------
Derivation or the Governing Differential Equation.
Consider a differential element of homogeneous aquifer as shown in
Figure 1. A material balance for any tracer can be written as
Rate of _ Rate of + Rate of-Mass _ Rate of Mass (1)
Mass In ~ Mass Out . Generation ' Accumulation
Each of the terms in this expression are developed in the following
paragraphs.
Rate of Mass In:
The tracer can enter the differential control volume by two mechanisms.
The first is convection, or bulk flow, of tracer in solution. This mechanism
can be expressed as
where Qv, Q., and Q7 are the volumetric flow rates in the x, y and z
A y L
directions, respectively, and C is the concentration of tracer in the fluid.
Tracer can also enter the control volume by molecular diffusion and
hydrodynamic dispersion. For the time being these processes can be combined
as a "dispersion flux" term, q" , with dimensions of mass input per unit time
per unit area. Thus,
AZAZ + q!J| AZAy
represents the rate of mass entering the control volume by "dispersion."
C-3
-------
QYC|
QZC|Z
qz|zAXAY
Qx<=|x
QXC!X
q" I AYAZ
A A
+AX
QYC!Y+AY
Figure 1 - Differential Control Volume for.Mass Balance
C-4
-------
Rate of Mass Out:
Tracer can leave the control volume by the same mechanisms- described for
the input terms. Thus, the rates of mass leaving by convection and dispersion
can be written as . .
and
respectively
Rate of Mass Generation:
Tracer can be generated (or degraded) in the control volume by physical,
chemical, and/or biological reaction. A general relationship can be written
as
AxAyAZ
where rj is the total, or overall, rate of generation per unit volume of
aquifer.
Rate of Mass Accumulation:
The total rate of mass accumulated in the control volume during a
differential period of time is
C-5
-------
where Cy is the total mass of trace per unit volume of aquifer.
Substituting the expressions for the various terms into Equation 1, and
dividing by (AxAyAz) yields
At
Now
Q:
Vx =
and
AX (AyAz) Ay UXAZ) AZ (AxAy)
Ax Ay- AZ
(CT),,A. - (C )
I t^-At I t
(3C)
where V is the superficial, or Darcy, velocity. Taking the limit of
Equation 2 as Ax, Ay, AZ, and At go to zero yields
C-6
-------
C) 3(V C) a(V C) 9q'' 3q" 3q"
x y £ *. y f_ + r =
3y 3z V 9y 3z T
Equation 4 is the differential mass balance for a tracer in a porous medium.
If a Fickian model is assumed-for the dispersion flux, these fluxes can
be expressed in terms of concentration gradients. Neglecting surface
diffusion of tracer which may be adsorbed on the solid matrix,
(5)
where n is the coordinate direction, Rn is a dispersion coefficient, and 0 is
the effective porosity or fractional void volume. The concentration gradient
in the liquid phase is the driving force for mass tranport by dispersion. -The
porosity has been included since the mass balance has been formulated for a
unit volume of aquifer, which includes the solid matrix as well as the fluid-
filled pores. Substituting an expression of the form of Equation 5.. for each
of the dispersion flux terms in Equation 4 yields
3CT a(v c) . a(v c) a(v c)
— L + - * — + - 1 — +
9x '3y 'z
Two assumptions will be made at this point. The first is that the
aquifer is homogeneous, which implies that the porosity and dispersion
coefficients are not functions of position. The second assumption is that the
ground-water flow is uniform and directed along the x-axis, i.e., V =
constant and V = Vz = 0. With these assumptions, Equation 6 reduces to
C-7
-------
T ar a C a C a r
+ r
This equation is a statement of conservation of tracer in homogeneous aquifer
with uniform ground-water flow. The accumulation and reaction terms are
developed in the following paragraphs.
Accumulation. In general, the total mass of tracer per unit volume of
aquifer, Cy, can be distributed as dissolved solute in the fluid-filled pore
volume and as adsorbed solute on the solid matrix, or
mass of tracer _ mass of tracer volume of solution
bulk volume volume of solution bulk volume
+ mass of tracer mass of so.l ids ,„»
mass of. solids bulk volume * '
Equation 8 can be written in terms of aquifer properties as
CT = 0C + PB Cs (9)
where pg is the bulk density of the solid matrix and C$ is the adsorbed mass
concentration (mass of tracer per unit mass of solids). For a homogeneous
aquifer, the accumulation term in Equation 7 can be written as
9CT 9C 9Cs
IT ' e £ + "B *T <10>
In general, the concentration of absorbed solute, Cs, is a function of
the concentration of solute in solution, C, and
C-8
-------
dC c
S dl
dC 8
For a linear adsorption isotherm,
or
Cs = KdC (12)
dCs
dC Kd
where K^ is a constant commonly referred to as an adsorption, or distribution,
coefficient. • The rate of accumulation of tracer per unit volume of aquifer
can be expressed in terms of the concentration in solution by combining
Equations 10, 11 and 13 as follows:
T
The coefficients of 3C/3t are often combined as
where R^ is referred to as a "retardation coefficient." Rewriting Equation 14
as
3CT
C-9
-------
or
3CT
3(t/Rd)
gives some insight into the effect of adsorption on accumulation of tracer in
the homogeneous aquifer. The apparent effect is a distortion or "retardation"
of the time dimension.
Reaction. Only first-order reactions will be considered in formulating
the rate of reaction term, rj, in Equation 7. The kinetic models for reaction
of tracer in solution and adsorbed tracer are
f=-xfc
and
3C,
3t S S
(17)
where \f and xs are the fluid phase and solid phase rate coefficients,
respectively. Equations 16 and 17 have been written for degradation of
tracer, i.e. negative generation. The overall rate of reaction can be written
as
9CT
rT -=HT= - exfC - pBXsCs
Including the linear adsorption isotherm developed above
C-10
-------
Cs = KdC
and
pBKd
(19)
or
rT = - GXyC
(20)
where Xj is an apparent overall first-order rate constant defined as
(21)
For radioactive decay, the rate of reaction is usually expressed in terms
of the "half-life" or the time required for the concentration to be reduced to
one-half of the initial concentration, 11/ . Integrating Equation 16 or 17
from t = 0 to t = tl/
C 12
.dJL
C
or
In C
CQ/2
= - X t
Solving for the rate constant,
C-ll
-------
x = 1n 2
which is an expression for evaluating a first-order rate constant from the
half-life of the reaction.
In the case of radioactive decay, the rate constants are independent of
the phase in which the reaction is occuring, and
X = Xf = Xs
Equation 21 can then be written as
;AT - Rd x
and Equation 19 becomes
PT = - 0 Rd XC (22)
Equation 22 applies to all cases where the first-order reaction rate constants
are the same for both fluid and solid phase reactions.
Differential Equation j_n_ Terms of F1 uid Phase Concentrations.
The differential mass balance for a tracer in a homogeneous aquifer with
uniform ground-water flow, Equation 7, can be written in terms of fluid-phase
concentrations by incorporating Equation 14 for linear adsorption and Equation
19 for first-order reactions. Making these substitutions and rearranging
yields
C-12
-------
• Dz -Z-j. - ATC (23)
where V* is the average interstitial, or pore, velocity defined as
V*=l (24)
Integration of Equation 23 with appropriate initial and boundary
conditions yields the temporal and spacial distribution of a tracer in a
homogeneous aquifer with uniform ground-water flow.
Closed-form analytical solutions to Equation 23 can be obtained by making
a change of variables. Let
T = t/Rd (25)
and
X = x - V*T (26)
Now, C = C(x,y,z,t), and holding y and z constant
Also,
C-13
-------
Substituting Equation 28 into Equation 27 yields
For the second derivative term in x,
Substituting Equations 25, 29 and 30 into Equation 23 yields
3X y 9y
- XTC (31)
T
which is a special form of .the heat conduction equation. Closed-form
analytical solutions for this equation are available .in the literature for a
variety of boundary conditions.
Summary.
Equation 23 provides the basis for many of the analytical solutions for
solute transport in uniform ground-water flow. This equation is a
mathematical model of complex physical phenomena and incorporates many
simplyfying assumptions which are required to obtain a solution to the
problem. Assumptions incorporated in the formulation of the differential
equation are also present in the solution and must be considered in
interpreting any numerical results.
The assumption that the aquifer is homogeneous is seldom satisfied in
practical field problems. Also, the use of an equilibrium adsorption isotherm
implies that adsorption of solute on the solid matrix is both reversible and
C-14
-------
instantaneous. Although these assumptions are seldom met in either field or
laboratory problems, the approximations to the physical system may be
reasonable for initial estimates of concentration distributions.
Solutions to Equation 23 with a continuous source of tracer are often
encountered in the literature. The assumption of a uniform velocity in the x-
direction makes no provision for the effects of a high- volumetric source rate
on the flow-field in the region of the source. For most problems, this
assumption is probably reasonable at moderate distances from either sources or
sinks of fluid, or on a regional basis.
The use of a Fickian model for the dispersive flux is probably the most
frequently misinterpreted or incorrectly applied portion of the mathematical
model. Hydrodynamic dispersion is an observed effect of one or more physical
phenomena which are difficult to define and cannot be measured. A discussion
of the topic is beyond the scope of this brief treatise. However, the
mathematical formulation of the problem as developed in this paper treats
dispersion as a potential flow problem. Analogous mechanism are conduction
in heat transfer and molecular diffusion in mass transfer. Thus, the model
does not distinguish between hydrodynamic dispersion in the direction of
ground-water flow or opposite to the direction of ground-water flow. The
model is a statement of conservation of tracer, and solutions of the governing
differential equation can lead to higher concentrations upgradient of sources
(and thus lower concentrations downgradient) than would be observed in
practice.
Analytical solutions to the solute transport equation present viable and
valuable alternatives for analyzing fairly complex problems, even with the
simplifying assumptions which have been incorporated. Interpretation of the
results of an analytical solution to the problem must be based on an
understanding of both the physical system and the mathematical model.
C-15
------- |