PB85-214450
PLUME 2D: TWO-DIMENSIONAL PLUMES IN UNIFORM GROUND
WATER FLOW
Oklahoma State University .
Stillwater, OK
Jun 85
U.S. DEPARTMENT OF COMMERCE
National Technical Information Service
NTIS
-------
EPA/600/2-85/065
June 1985
PLUME2D .
TWO-DIMENSIONAL PLUMES IN UNIFORM GROUND WATER FLOW
by
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 Instruction! on the reverse before completing)
1. REPORT NO.
EPA/600/2-85/065
2.
3. RECIPIENT'S ACCESSION NO.
c-:";. ••'. -.;; o /AS
4. TITLE AND SUBTITLE
PLUME2D: TWO-DIMENSIONAL PLUMES IN UNIFORM GROUND
WATER FLOW
5. REPORT DATE
June 1985
6. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
8. PERFORMING ORGANIZATION REPORT NO.
Jan Wagner, Stephanie Watts, Douglas Kent
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Oklahoma State University
Stillwater OK 74078
10. PROGRAM ELEMENT NO.
ABRD1A
Coop. Agr.
CR-811142
12. SPONSORING AGENCY NAME AND ADDRESS
Robert S. Kerr Environmental Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Ada, OK 74820
13. TYPE OF REPORT AND PERIOD COVERED
Final 9/83 - 2/85
14. SPONSORING AGENCY CODE
EPA/600/15
15. SUPPLEMENTARY NOTES
Carl G. Enfield, Project Officer
16. ABSTRACT
A closed-form analytical solution for two 'dimensional plumes was incorporated
in an interact! ve. computer program. The assumption of an infinite aq.uifer 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. The model can be solved for either vertically or horizontally averaaed
conditions. 3
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS C. COSATI Field/Group
18. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (This Report/
UNCLASSIFIED
21. NO. OF PAGES
95
20. SECURITY CLASS (This page}
UNCLASSIFIED
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
11
-------
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. The model can be solved for either vertically
or horizontally averaged conditions.
-------
TABLE OF CONTENTS
INTRODUCTION 1
SECTION I - MATHEMATICAL DEVELOPMENT 2
. Vertically-averaged solution.. 7
Horizontally-averaged solution ; 13
Assumptions and Limitations 17
Superposition 17
SECTION II - COMPUTER PROGRAM 25
Basic Input Data • 25
Edit Commands 30
SECTION III - APPLICATIONS .... 32
REFERENCES 40
APPENDIX A - Example Problems , A-l
APPENDIX B - Description of Program PLUME2D B-l
APPENDIX C - Numerical Analysis of the Hantush Well Function C-l
APPENDIX D - Listing of Program PLUME2D .' D-l-
APPENDIX E - Listing of Utility Function Subroutines...... : .• E-l
-------
LIST OF FIGURES
Figure 1 - Coordinate system for vertically averaged solution 6
Figure 2 - Coordinate system for horizontally averaged solution 6
Figure 3 - Decomposition of a variables source rate using
superposition in time 20
Figure 4 - Use of image sources to account for aquifers of
finite depth 22
Figure 5 - Superposition in space to account for barrier boundaries......24
Figure-6 - Results of hexava.lent chromium plume simulation at 3280 days..33
Figure 7 - Results of hexavalent chromium spill simulation at 365. days.'. .36
-------
LIST OF TABLES
Table 1 - Edit Commands..... 31
Table 2 - Comparison of Concentrations Calculated Using Superposition
in Time and an Analytical Solution for an Instantaneous
Line Source. 38
Table 3 - Typical Values for Aquifer Properties 39
vn
-------
INTRODUCTION
Relatively simple analytical methods can often be used to evaluate
ground-water contamination problems, depending upon the complexity of the
system and the availability of field data. • Analytical models can also serve
as valuable tools in developing parameters for more sophisticated numerical
models. Although the numerical evaluation of an analytical solution to a
ground-water probl-em may be mathematically complex, analytical models are well
suited for interactive use on digital computers. Many analytical solutions to
ground-water contamination problems can be coded on programmable hand-held
calculators. In general, very few input parameters are required to define a
given problem and numerical results can be calculated in a few seconds.
•This report presents analytical solutions to two ground-water .pol lution
problems — two-dimensional -plumes in uniform ground-water flow. An
interactive computer code has been developed which enables the user to modify
the definition of a given problem, and thus gain some insight into the effects
of various parameters on the extent of a contaminant plume.
-------
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
r
where
C = component mass per unit of fluid phase ' M/L
•3
Cy = total component mass per unit volume of aquifer M/L°
o
DX = dispersion coefficient in x-direction . L /t
D = dispersion coefficient in y-direction LVt
p
DZ = dispersion coefficient in z-direction L /t
rt = rate of degradation of mass per unit volume
of aquifer ' • M/L3t
.V = Darcy, or seepage, velocity in the x-direction L/t
x,y,z, = rectangular coordinates at the point of interest L
o o
Q = porosity of porous media L°/L
-------
The total mass of a component per unit volume of aquifer is distributed
as dissolved solute in the fluid phase and adsorbed solute on the solid
matrix. Let
Cs = component mass per unit mass of solid M/M
and
pg = bulk density of the aquifer, or the mass of solids
per unit volume of the aquifer M/L .
The total component mass per unit volume of aquifer can be expressed
as
Mass _ Volume of voids Component Mass
Unit Volume ~ Unit. Volume of aquifer Volume of voids
+ Mass of solids Component Mass
- • Unit volume of aquifer Mass of sol-ids
or
CT = G'C .+ PB Cs ' (2)
and, the rate of accumulation of mass in the aquifer becomes
3C
In general, Cs = f(C) and
aC dC „-
at dC at
For a linear equilibrium adsorption isotherm,
dCs = K M/M
d M/L3
-------
where K^ is a distribution constant.
The change in concentration per unit v.olume of porous media, sCy/at, can
be written in terms of fluid phase concentration, C, by substituting Equations
4 and 5 into Equation 3. Therefore,
3CT _ -3C _3C
"at"" e~3t. pBKd at
or
3C
The rate of degradation of component mass per unit volume of porous media
is also distributed between the solid and liquid phases, or
Rate of mass'degraded _ • Rate of mass degraded Volume of fluid
Unit volume of aquifer Unit vo.lume of fluid Volume of aquifer
Rate of mass degraded Mass of sol id •
-.' Unit mass of solid ' Volume of aquifer'
Now, the rate of change in total mass per unit volume of aquifer due to
reaction can be written as
3CT 8C 3Cs
rt - ar = 9lr+ ptf^r
The concentration on the solid, Cs, is related to the concentration in the
liquid, C, through the linear adsorption isotherm assumed previously, and
-------
Assuming first order decay kinetics, the rate of decrease in fluid phase and
solid phase concentrations due to reaction can be expressed as
and
8C
S = XC.
3t S
respectively, where A is a rate constant (1/t), and
rt = (6 + PB"
-------
00
Y O-O
00
GROUND-WATER FLOW
CD
Figure 1. Coordinate system for vertically averaqed solution.
WATER 0
0-6-
TABLE
00
GROUND-WATER FLOW
00
Figure 2. Coordinate system for horizontally averaqed solution.
-------
Equation 13 is a linear partial differential equation which can be
integrated analytically to yield an expression for concentration as a function
of time and position.
Solutions of Equation 13 for two types of ground-water contamination
problems are presented in the following paragraphs. The first is a
vertically-averaged solution which describes a contaminant plume in the x-y
plane (Figure 1). The second is a horizontally-averaged solution, describing
a contaminant plume in the x-z plane (Figure 2).
Vertically-averaged solution. The vertically-averaged solution applies
to a homogeneous aquifer of infinite aerial extent and finite depth. The
contaminant is assumed to be well mixed over the saturated thickness. The
source of contaminant is a vertical line source located at the origin of a
coordinate system in the x-y plane. This conceptual model would apply to an
injection well which fully penetates the saturated zone.
Wilson and Miller (1978) have also applied this solution downstream from
a contaminant source at the surface of the water table. For a relatively thin
saturated zone, vertical dispersion will result in mixing vertically. The
concentration distribution can be considered as being two-dimensional in a
horizontal plane at distances downstream of the source sufficient for the
concentration distribution to become uniform with depth. Mathematically, the
problem is treated as an infinite aquifer with a line source at the origin.
The vertically-averaged formulation of Equation 13 is
Rdf +V*f
3x J 3y
-------
The boundary conditions can be stated mathematically as follows
C(x,y,0) = 0
(15a)
C(x,±-,t) =.0
(15b)
C(±-,y,t) = 0
(15c)
A solution to Equation 14 with the Equation 15 boundary conditions and a
continuous source of strength C0Q' can be written as (Hunt, 1978; Wilson and
Miller, 1978):
0Q' EXP (
C =
4ir0 (D D
x y
.0.5.
W(U,B)
(16)
where
U =
Vxf . °x
H + IT
x/ _ y
4 V t
R.I D
d x
(17)
and
B
1
2
* *
/ V x\
( D J
\ x /
, Dx
D
y
f^) '
\DJ
1C.
4DxRdX"
V*2
(18)
-------
and CQQ' (M/t/L) is the contaminant source rate per unit depth of the
saturated zone.
The function W(U,B) is defined as
» i R2
W(U,B) = / .EXP (-5 -) dc (19)
where c is a dummy integration variable. This function is often referred to
as the "well function for leaky artesian aquifers" (Hantush; 1956, 1964)'.
The steady-state solution of Equations 14 and 15 can be obtained by
noting as t -»• «, U •»• 0 and the well function (Hantush, 1956) can be expressed
as
W(0,B) = 2KQ(B) . •' . (20)
where KQ(B) is the modified Bessel function of the second Id/id of order zero.
At steady-state the vertically-averaged solution can be written as
C0Q' EXP
C= Q-f— K (B) (21)
(D D )U-b °
x y'
The units of the variables in Equations 16 and 21 .can be eliminated by
defining the following dimension.less groups:
Modified Peclet Numbers Convectlve mass transport
Dispersive mass transport
-------
Damkohler Group II • ~ Mass decay rate
K Mass dispersion rate
DRx
Number of Pore Volumes Injected ~ !'!ass transport rate
J Mass accumulation rate
d x
ion, ess Source Te™ - 'rill
(22)
(23)
0(D D )0.5
v x y/
Dimension!ess Concentration
(26)
10
-------
(27)
Note that the number of pore volumes injected can be written as
I =
V2t /vV 2
DxRd"
RdL
(28)
where L is a characteristic length defined as
,2 2 x x 2
L = x + -=- y
(29)
The first group on the right-hand-side of Equation 28 is the Modified
Peclet number
Pe
D ; D AD /
x / y v x '
or
Pe/)
(30)
11
-------
where
The second group on the right-hand-side of Equation 28 is a dimensionless time
variable,
D t
T = -2W- (31)
- RdL2 .
The transient and steady-state solutions to Equations 14 which are given by
Equations 16 and 21 can be written in terms of the dimensionless variables
defined above. The transient solution is
= -EXP( Pe) W(U,B) (32)
and at steady state
Y = |^EXP (\ PeJ KQ(B) (33)
with
12
-------
Pexv2
U --IT^T (34)
and
<35)
The values of dimension!ess concentrations evaluated using Equation 32 or
Equation 33 are val.id for any consistent set of units. Using dimension! ess
variables also tends to "scale" numerical values when working in various
systems of units. ' '
Horizontally-averaged solution. Consider a homogeneous aquifer-with a
continuous line source of infinite length.located at the water table and
normal to the direction of ground-water flow as shown in Figure 2. In other
words the tracer is assumed to be well mixed over the width of the aquifer. A
problem which might fit this conceptual model is seepage from a trench
perpendicular to the direction of ground-water flow.
The horizontally-averaged formulation of Equation 13 is
Rdf +V*f =Dx+Dz-RdXC <36>
13
-------
For an aquifer of infinite depth and a uniform continuous line source, the
appropriate boundary conditions can be written as follows
C(x,z,0) = 0 (37a)
C(±-,z,t) - 0 (37b)
C(x,-,t) = 0 ' (37c)
%C* { v 0 "t* ^
OwlA»U»l»J_f^
TZ ~ U
A solution to Equation.36 with the Equation 37 boundary conditions and a
continuous 1ine source of strength C0Q' is
CQQ'
C = - -r- W(U,B) (38)
2*6
At steady state the horizontally-averaged solution can be written as
C0Q' EXP
where
14
-------
U =
vY2
x /V
Dz \Dx
*?
4 V * t
Rd°x
(40)
V x
\ /V 2
D ID
z \ x
2 -,
(41)
and Q1 (L^/t/L) is the volumetric contaminant source rate per unit width of
the aquifer (or unit length of the line source).
Changing.subscripts, the definition of the dimensionless groups leads to
"•'
and
(42)
r =
Q'
9(0 D )0.5
x z
with
2 2
(43)
(44)
where
(45)
15
-------
By substituting the dimensionless groups described in vertically-averaged
solution and those defined above, Equations 38 through 41 can be written in
terms of dimensionless variables.
The transient solution becomes
Y =-Exp( Pex) W(U«
and at steady state, the horizontally-averaged solution is
KQ(B) (47)
where .
2
Pe
U = - .(48)
and
B= j Pexz(l + 4Dk)7 (49)
The similarity of the solutions of the vertically-averaged and horizontally-
averaged problems facilitates their numerical evaluation using a common
computational algorithm. For the same numerical values of the independent
variables, concentration values.for the horizontally-averaged solutions are
obtained by doubling the vertically-averaged solution values.
16
-------
Assumptions and Limitations
Equations 32-33 and 46-47 can be used to calculate the concentrations in
leachate plumes 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. All ground-water flow is horizontal, continuous, and uniform
throughout the aquifer.
4. The aquifer is infinite in extent for the vertically-averaged
solution, or semi-finite in extent for the horizontally-averaged
solution.
5. The leachate source is a line 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"
17
-------
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.
Both the horizontally-averaged and vertically-averaged solutions can be
used to simulate aquifers of finite width or depth, respectively, or plane
sources of finite.width. Applications of this type require a thorough
understanding of the physical interpretation of the principal;of
superposition.
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 3. The
solutions of the governing differential equation presented in this report are
of the form
C(x,z,t) = CoQ' f(x,z,t) = Q1 f(x,z,t) (50)
•
where Q' is the source mass rate per unit length. The principle of
superposition in time can be written for any position as
18
-------
C(x,z,t) = £ Q.1 f(x,z,t.) (51)
1=1
Now, the variable rate schedule shown in Figure 3a can be decomposed into a
series of positive and negative mass rates as shown in Figure 3b. 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 ffx.y.z.tj) - Qj f(x,y,z,t2)
Q2 f(x,y,z,t2) - Q2' f(x,y,z,t3)
Q3 f(x,y,z,t3) - Q3 f(x,y,z,t4)
^ f(x,y,z,t4) . (52)
In general terms
C(x,y,z,t ) = E (: - Q; , ) f(x,y,z,t.) (53)
s .=1 i 1-1 i
with QQ = 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 53 can be rewritten as
n
C(x,y,z,t ) = £ (Q1 - Q1 ) f(x,y,z,t -t ) (54)
O i _-i ^ N ~ i b N~1
19
-------
•o
uf
cc
CO
co
2
(
Qi
63
61
6;
• p
PERIOD OF
SIMULATION
0 TIME, t
(a)
>o
LJLj"
cr
CO
CO
t
ov.
U3
6"
U2
0"
U1
0"
U4
n%
u 1
Q'o
U2
6',
U3
.
•
"
.
:
•
"
! TIME, t
!
i
Ut4^
^ * „ ^
•" 13 *
*• ^_
12 *
t. _ .. .. .—
(b)
Figure 3. Decomposition of a variable source rate using superposition in time.
.20
-------
where t^ is the time corresponding to the end of mass rate 0^.^ 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
•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 4a. If the contaminant plume was to
intersect an impermeable base of the aquifer as shown in Figure 4b, 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
D -Ir- = 0
z 8z
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 4c, this source would create a concentration gradient from the boundary
to the image water table equal to the concentration gradient from the boundary
21
-------
WATER TABLE
WATER TABLE
(a)
(b)
(d)
Figure 4. Use of image sources to account for aquifers of finite depth.
22
-------
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 4d.
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
5. 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
only a few image systems are required. The computer program automatically
introduces an appropriate number of image systems.
23
-------
no
-pa
INFINITE
AQUIFER
IMAGE
CONCENTRATION
iWATER TABLE
1st
IMAGE SYSTEM
2nd
IMAGE SYSTEM
c(x, z)
2d-rr,
c(x, z)
oc(x, 2d-z)
IMAGE___
WAf ER TABLE
4d
IMAGINARY^
BOUNDARY
oc(x, 2d+z)
vvv vv v\ v w v vv \ vwwvvv
'c(x, 4d-z)
6d
00
c(x, 4d+z)
/ /~r~f~y r f
oc(x, 6d-z)
c'(x, z) = c(x, z)^e(x, 2nd-z)-i-c(x, 2nd*z)]
n=1
Figure 5. Superposition in soace to account for barrier boundaries.
-------
SECTION II
COMPUTER PROGRAM
The computer program evaluates the analytical solutions of the
differential equation describing concentration distributions in two-
dimensional plumes 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 PLUME2D
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 cornma(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 second input command is used to select the vertically-averaged
solution or the horizontally-averaged solution. The command is:
ENTER COORDINATE SYSTEM
XY FOR VERTICALLY-AVERAGED SOLUTION
XZ FOR HORIZONTALLY-AVERAGED SOLUTION
?
Either of the indicated responses is valid.
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)
25
-------
Any valid keyboard characters can be used. The first two characters will be
retained for identifying the units of the length dimensions which may be
required for other input data or output listings.
ENTER UNITS FOR TIME (2 CHARACTERS)
?
Any valid keyboard characters can be used. The first two characters will be
retained for identifying the units of the time dimensions which may be
required for other input data or output listings.
ENTER UNITS FOR CONCENTRATION (6 CHARACTERS)
7
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 initialize 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
?
If horizontally-averaged solution was selected (x-z coordinate system) this
request is issued. 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
7
Enter the volume void fraction.
ENTER SEEPAGE VELOCITY, L/t
7
26
-------
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
7
The retardation coefficient includes the effects of absorption of the tracer
on the solid matrix (see Section I for discussion). The numerical value must
be greater than 1.0, or equal to 1.0 if absorption is neglected.
ENTER X DISPERSION COEFFICIENT, SQ L/t
Dispersion coefficients have dimensions of L/t and must be entered in the
units requested. Numerical values must be greater than zero.
If the X-Y coordinate system has been selected, the next command is:
ENTER Y DISPERSION COEFFICIENT, SQ L/t
If, instead, the X-Z coordinate system has been selected, a command for the Z
dispersion coefficient will be issued.
ENTER Z DISPERSION COEFFICIENT, SQ L/t
The subsequent command will be:
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 (see Section I for
discussion). Approximation is accomplished through superposition of a series
27
-------
of uniform rates. If steady-state solution is chosen,'the steady state
concentration will be evaluated.
ENTER THE NUMBER OF SOURCES (MAXIMUM OF N)
7
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
This statement reminds the user of the units that will be used for mass rates
and-for time. All mass-rate and time values entered must be in these units.
The next series of commands will be repeated for each source.
ENTER X and I 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, instead, the X-Y
coordinate system has been selected, the fol1 owing.command is issued:
ENTER X AND Y COORDINATES OF SOURCE I (L)
?
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.
28
-------
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 two 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.
XFIRST and XLAST can be positive or negative values. A zero entry for DELTAX
will result in a single X-coordinate observation. Results of calculations for
multiple X-coordinates will be listed from. XFIRST to XLAST.
ENTER YFIRST, YLAST, DELTAY (L)
• 9 ' 9 •
Any of the numerical values Jsed to define the Y-coordinates of
observation points may be. positive or negative. If the X-Z coordinate system
has been selected, a command to enter the Z-coordinates, rather than the Y-
coordinates, will be issued.
ENTER ZFIRST, ZLAST, DELTAZ (L)
ENTER TFIRST, TLAST, DELTAT (t)
• i •»•
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.
29
-------
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 ed.it 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?
One of the reponses from Table 2 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 wil1 be
repeated until one of the responses Mil, LI, RN, NP, or DM is 'entered.
ML) will list the table of edit commands.
LI will- list the problem as currently defined.
. RN will initiate-the calculation of concentrations and print the results.
NP will request a complete new problem using the "Basic Input Data"
dialog.
ON will terminate the program.
A listing of the dialog and the results for the example problem discussed
in Section III are included.in Appendix A.
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.
30
-------
Table 1
EDIT COMMANDS
Command Variable changed/Execution
ST Saturated Thickness
PO Porosity
VX • New Seepage Velocity
RD Retardation Coefficient
DE . Decay Constant
DX X-Dispersion Coefficient
DY Y-Dispersion Coefficient
DZ Z-Dispersion Coefficient
RT Source Rate Schedule
08 Observation Points
XC X-Coordinates
YC Y-Coordinates
ZC Z-Coordinates
TC Observation Times
CS Change Solution/Sources
ML) Menu of Edit Commands
LI List input data
RN . Run
NP New Problem
DN Done
31
-------
SECTION III
APPLICATIONS
The case history of ground-water contamination with hexavalent chromium
in South Farmingdale, Nassau County, New York, (Perlmutter and Lieber, 1970),
has been used as an example of the application of the two-dimensional plume
model. The contaminant plume has been modeled numerically by Pinder (1973)
and analytically by Wilson and Miller (1978). Details of the hydrologic
system are described in the above references. A brief summary of the problem
is presented in the following paragraphs.
The aquifer is assumed to have a saturated thickness of 33.52 m with a
porosity of 0.35. Perlmutter and Lieber (1970) estimated the average seepage
velocity to be approximately 0.366 m/dy. Using Pinder's (1973) values of
dispersivity, aY = 21.3 m and a., = 4.27 m, and x and y dispersion coefficients
A J
are
D = (21.3 m) (0.366 m/dy) = 7.70 m2/dy
and
D = (4.27 m) (0.366 m/dy) = 1.56 m2/dy
The source of contamination consisted of three metal-plating-waste disposal
ponds as shown in Figure 6. The mass rate of chromium entering the aquifer
has been estimated at 23.6 kg/dy during the nine year period from 1941 through
1949 (Perlmutter and Lieber, 1970). Chromium is believed to be a conservative
contaminant, thus absorption and degradation can be neglected. The
vertically-averaged model parameters are summarized as follows:
32
-------
OJ
co
400
200
cc
in
t-
LU
-200
AERAL EXTENT OF HEXAVALENT CHROMIUM PLUME
(PERLMUTTER AND LIEBER. 1970)
DISPOSAL
BASINS
400
800
X, METERS
1200
1mg/l
1600
Figure 6. Results of hexavalent chromium plume simulation at 3280 days.
-------
Aquifer porosity 0.35
Seepage velocity 0.366'm/dy
Retardation coefficient 1.0
o •
x-Dispersion coefficient 7.79 nr/dy
y-Dispersion coefficient 1.56 nr/dy
Decay constant 0.0 1/dy
The contaminant source rate is assumed to be constant, and only one rate
period is required. The mass rate can be. converted to units of concentration
times volume rate per unit depth as
23.6 kg 10. mg m 1 ,n/1 , ,, . , 3,, , .
~dy-^ —kg ^3[ 33.52 m = 704 ^/l)
-------
Rate 1: 704 (mg/1) (m3/dy)/m from 0 to 1 day
Rate 2: 0 (mg/1) (m3/dy)/m from 1 to 365 days
Other model parameters are identical to those used in the previous example.
The results of the simulation are summarized in. Figure 7, which shows the
center of mass of the plume moving down-gradient at the seepage velocity and
spreading longitudinally and transversely by diffusion.
The results of the" simulation using superposition in time were compared
with the concentrations calculated using
4*8t(D D )< x - y
x y J •
which is the solution of Equation 14. for an instantaneous line source of
strength CQQ' (m/l3) (L3/L). The' values of concentration and errors in
approximating the instantaneous source through superposition in time are
presented in Table 2. Note that the finite duration of the source results is
si ightly higher concentrations up-gradient from the center of mass than
concentrations down-gradient. For an instantaneous source the concentration
distribution should be symmetrical about a y-z plane through the center of
mass located as x = Vt. A better approximation can be obtained by injecting
the same total mass of contaminant over a shorter period of time; but for
purposes of illustrating superposition in time, the errors in the example
problem are not significant.
The example problems presented above are intended to illustrate the
appl ication of the two-dimensional plume models developed in this report.
These models are tools which can aid in the analysis of ground-water
contamination problems. The user must select the best tool for the problem at
35
-------
oo
en
CO
QC
111
H
Ul
60
40
20
0
-20
-40
-60
0
40
0.05 mg/l
i''' 1 1 L
80 120
X, METERS
160
200
240
Figure 7. Results of hexavalent chromium spill simulation at 365 days.
-------
hand, based on a sound understanding of the principles of ground-water
hydrology, the physicalproblem, and the limitations of the mathematical
model(s).
Perhaps the most difficult step in using any mathematical model is
defining the problem, to be solved. In addition to developing the physical
boundaries of the problem domain, rock and fluid properties must also be
quantified. Typical values of aquifer properties are listed in Table 3, but
the user must accept the responsibility for developing the required model
input data for the specific problem to be solved.
37
-------
y
(meters)
20.0
10.0
0.0
Table 2
COMPARISON OF CONCENTRATIONS CALCULATED
USING SUPERPOSITION IN TIME AND AN
ANALYTICAL SOLUTION FOR AN INSTANTANEOUS
.LINE SOURCE
MODEL PARAMETERS
Aquifer Porosity 0.35
Seepage Velocity 0.366 m/dy
Retardation Coefficient . 1.0
x-Dispersion Coefficient 7.79 m2/dy
y-Dispersion Coefficient 1.56 nr/dy
First-Order Decay Constant 0.0 1/dy
Source Strength . 704.0 (mg/1) (m3/tn)
Concentration at 365 days, mg/1
Superposition
(Equation 54)
% Error
x(meters)
73.59
0.0919
( .0917)
.22
0.0878
( .0877)
.11
0.0771
( .0769)
.26
103.59
0.1165
( .1162)
.26
0.1115
( .1112)
.27
0.0977
( .0975)
.21
133.59
(= vt)
0.1259
( .1258)
.08
0.1206
( .1204)
.17
0.1057
( .1055)
.19
163.59
0.1163
( .1162)
.01
0.1113
( .1112)
.09
0.0975
( .0975)
.0
193.59
0.0916
( .0917)
-.11
0.0876
( .0877)
-.11
0.0768
( .0769)
-.13
38
-------
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 0.01 - 0.1
Material
Clay
87.36 - 137.2
0.03 - 0.05
0.01 - 0.1
0.1 - 1.0
0.01 - 0.1 .
0.1 - 1.0
Silt
80.50 - 112.3
0.05 - 0.10
1 - 10
1 - 10
0.1 - 1.0
1 - 10
Sand
7-3.63 - 98.59
•0.10 - 0.30
100 - 100,000
10 - 100
1.0'- 10. •
39
-------
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.
Hantush, M. S., 1956, "Analys.is of Data from Pumping Tests in Leaky Aquifers,"
Transitions, American Geophysical Union, Vol. 37, No. 6, pp. 702-714.'
Hantush, M. S., 1964, "Hydraulics of Wells," Advances in Hydroscience, Vol. 1,
pp. 281-432.
Hantush, M. S. and C. E. Jacob, 1955, "Non-Steady Radial Flow in an Infinite
Leaky Aquifer," Transitions, American Geophysical Union, Vol. 26., No. 1,
pp. 95-100.
Hunt, B., 1978, "Dispersive Sources in Uniform Ground-Water Flow," Journal of
The Hydraulics Division, ASCE, Vol. 104, No. HY1, 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 Galerkfn Finite Element Simulation 'of Ground-water
Contamination on Lo-ng Island," Water Resources Research, Vol. 9, No. 6,
pp. 1657-1669.
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-Hi11, New York,
New York, 664 pp.
Wilson, J. L. and P. J. Miller, 1978, "Two-Dimensional Plume in Uniform
Ground-Water Flow," Journal of the Hydraulics Division, ASCE, Vol. 104,
No. HY4, pp. 503-514.
Wilson, J. L. and P. J. Miller, 1979, "Two-Dimensional Plume in Uniform
Ground-Water Flow, Discussion," Journal of the Hydraulics Division, ASCE,
Vol. 105, No. HY12. pp. 1567-1570.
Yeh, G. T., 1981, "AT123D: Analytical Transient One-, Two-, and Three-
Dimensional Simulations of Waste Transport in the Aquifer System,"
Publication No. 1439, Environmental Sciences Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee, 83 pp.
40
-------
APPENDIX A
Example Problems
The two example problems presented in the following pages are discussed
in Section III of. this report. The first demonstrates .the application of
PLUME2D to a continuous source of contamination. The second example
approximates an instantaneous source using the principle of superposition in
time as discussed, .in Section I.
A-l
-------
ENTER TITLE .
7HEXAVALENT CHROMIUM 'PLUME
ENTER. COORDINATE SYSTEM
XY FOR VERTICALLY-AVERAGED SOLUTION
XZ FOR HORIZONTALLY-AVERAGED SOLUTION
?XY
ENTER UNITS FOR LENGTH (2 CHARACTERS)
? M
ENTER UNITS FOR TIME (2 CHARACTERS)
?DY
ENTER UNITS FOR CONCENTRATION (6 CHARACTERS)
7MG/L
ENTER AQUIFER POROSITY
?0.35 -
ENTER SEEPAGE VELOCITY, M/DY
70.366
ENTER RETARDATION COEFFICIENT
?1 .0
ENTER X DISPERSION COEFFICIENT, SQ M/DY
?7.79
ENTER Y DISPERSION COEFFICIENT, SQ M/DY
71.56
ENTER DECAY CONSTANT, 1/DY
70.0
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 M/DY)
.TIME HAS UNITS OF DY
ENTER X AND Y COORDINATES .OF SOURCE 1 ( M)
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
?,7704. ,3280.
ENTER XFIRST, XLAST, DELTAX ( M)
?,?,?200. ,1200.,200.
ENTER YFIRST, YLAST, DELTAY '( M)
?,?,?200. ,-200. ,50. ' . •'
ENTER TFIRST, TLAST, DELTAT (DY)
?,?,73280.,0.,0.
A-3
-------
PLUME2D
VERSION 2.01
PAGE 1
HEXAVALENT CHROMIUM PLUME
SEEPAGE VELOCITY, ( M/DY)
X DISPERSION COEFFICIENT ( M**2/DY)
Y DISPERSION COEFFICIENT ( M**2/DY)
POROSITY
.3660
7.7900
1.5600
.3500
RETARDATION COEFFICIENT
FIRST ORDER DECAY CONSTANT (1/DY)
1 .0000
0.0000
SOURCE/RATE SCHEDULE (MG/L )(CU M/DY)
NO
X
SOURCE
M) Y ( M)
RATE
NO
MASS
.RATE
TIME (.DY)
START END
0 .00
0 :00
l •
704.00
0.00 3280.00
OBSERVATION POINTS ( M)
XFIRST =
YFIRST =
20.0 ..00
200.00
XLAST
YLAST
OBSERVATION TIMES (DY)
TFIRST = 3280.00 TLAST
1200.00
-200.00
DELX
DELY
3280.00 DELT =
200 .000U
50.0000
0 .0000
A-4
-------
MENU OF EDIT COMMANDS
RETARDATION COEFFICIENT RD
POROSITY PO
SEEPAGE VELOCITY VX
X DISPERSION COEFFICIENT DX
Y DISPERSION COEFFICIENT DY
DECAY CONSTANT DE
SOURCE RATE SCHEDULE RT
NEW PROBLEM. NP
CHANGE SOLUTION/SOURCES CS
OBSERVATION POINTS OB
X COORDINATES XC
Y COORDINATES YC
MENU OF COMMANDS MU
LIST INPUT DATA LI
RUN CALCULATIONS RN
DONE ' DN
SATURATED THICKNESS ST
OBSERVATION TIMES
TC
ENTER NEXT COMMAND
?RN
A-5
-------
PLUME2D
VERSION 2.01
PAGE 2
HEXAVALENT CHROMIUM PLUME
CONCENTRATION DISTRIBUTION AT
3280.00 DY (MG/L )
* X( M)
*
Y( M) *
*
200 .00
150.00
100.00.
50.00
0.00
-50.00
' '-100 ..00
-150.00 •
-200.00
200.00
.0372
.4289
4.0806
24.5165
51.8245 .
24.5165
4.0806
.4289
.0372
400 .00
.2773
1.8560
8.8387
25.3968
37.0664
25.3968
8.8387
1.8560
.2773
600.00
.8210
3.6177
11 .3609
23 .5539
30.2812
23.5539
11 .3609
3.6177
.8210
800 .00
1.4371
4.8444
11.9818
20 .9946-
25.3930
20.9946
11 .9818
4.8444
1.4371
1000 .00
1.6352
4.7217
10.2348
16.4014
19.2190
16.4014
10.2348
4.7217
1.6352
1200.00
1.13.80
3.0238
6.1201
9.372.1
10.8087
9.3721
6.1201
3.0233
1.1380
ENTER NEXT COMMAND
?DN
STOP
A-6
-------
ENTER TITLE
7ACCIDENTAL HEXAVALENT CHROMIUM SPILL
ENTER COORDINATE SYSTEM
XY FOR VERTICALLY-AVERAGED SOLUTION
XZ FOR HORIZONTALLY-AVERAGED SOLUTION
?XY
ENTER UNITS FOR LENGTH (2 CHARACTERS)
? M . .
ENTER UNITS FOR TIME (2 CHARACTERS)
?DY
ENTER UNITS FOR CONCENTRATION (6 CHARACTERS)
7MG/L
ENTER-AQUIFER POROSITY
?0.35
ENTER SEEPAGE VELOCITY, M/DY
70.366
ENTER RETARDATION COEFFICIENT
ENTER X DISPERSION COEFFICIENT, SQ M/DY
?7.79
ENTER Y DISPERSION COEFFICIENT, SQ M/DY
71.56
ENTER DECAY CONSTANT, 1/DY
70.0
SELECT TRANSIENT OR STEADY-STATE SOLUTION
TR FOR TRANSIENT SOLUTION
SS FOR STEADY-STATE SOLUTION
7TR
A-7
-------
ENTER THE NUMBER OF SOURCES (MAXIMUM OF 10 )
MASS RATES HAVE UNITS OF (MG/L ) (CU M/DY)
TIME HAS UNITS OF DY
ENTER X AND Y COORDINATES OF SOURCE 1 ( M)
ENTER THE NUMBER RATES FOR SOURCE 1 (MAXIMUM 0F 10)
?2
SOURCE 1, RATE 1 STARTS AT 0.0 DY
ENTER MASS RATE AND ENDING TIME
SOURCE 1, RATE 2 STARTS AT' 1.0 DY
ENTER MASS RATE AMD ENDING TIME
?,?0. ,365.
ENTER XFIRST, XLAST, DELTAX ( M)
?,?,773.59,193.59,30.
ENTER YFIRST, YLAST, DELTAY ( M)
?,?,?20.,0.,10.
ENTER TFIRST, TLAST, DELTAT (DY)
?,?,?365.,0. ,0.
A-8
-------
PLUME2D
VERSION 2.01
PAGE 1
ACCIDENTAL HEXAVALENT CHROMIUM SPILL
SEEPAGE VELOCITY, ( M/DY)
X DISPERSION COEFFICIENT ( M**2/DY)
Y DISPERSION COEFFICIENT ( M**2/DY)
POROSITY
.3660
7.7900
1.5600
.3500
•RETARDATION COEFFICIENT
FIRST ORDER DECAY CONSTANT (1/DY)
1 .0000
0.0000
SOURCE/RATE SCHEDULE (MG/L )(CU M/DY)
' SOURCE
NO X ( M) ' Y ( M)
RATE
NO
MASS
RATE
• TIME'(DY)
START END
0.00
0 . 00
1 704.00 0.00 1.00
2 0.00 1.00 365.00
OBSERVATION POINTS ( M)
XFIRST =
YFIRST =
73.59
20.00
OBSERVATION TIMES (DY)
TFIRST = 365 .00
XLAST
YLAST
TLAST =
193.59
0.00
DELX
DELY
365.00 DELT =
30 .0000
10.0000
0.0000
A-9
-------
MENU OF EDIT COMMANDS
RETARDATION COEFFICIENT RD
POROSITY PO
SEEPAGE VELOCITY VX
X DISPERSION COEFFICIENT DX
Y DISPERSION COEFFICIENT DY
DECAY CONSTANT DE
SOURCE RATE SCHEDULE RT
NEW PROBLEM NP
CHANGE SOLUTION/SOURCES CS
OBSERVATION POINTS OB
X COORDINATES XC
Y COORDINATES YC
MENU OF COMMANDS MU
LIST INPUT DATA LI
RUN CALCULATIONS RN
DONE DN
SATURATED THICKNESS ST
OBSERVATION TIMES TC
ENTER NEXT COMMAND
?RN
A-10
-------
PLUME2D
VERSION 2.01
PAGE 2
ACCIDENTAL HEXAVALENT CHROMIUM SPILL
CONCENTRATION DISTRIBUTION AT 365.00 DY (MG/L )
* X( M)
Y(
*
M) *
*
20.00
10.00
0 .00
73.59
.0771
.0879
.0919
103.59
.0977
.1115
.1165
133.59
.1056
.1204
.1260
163.59
.0975
.1113
.1163
193.59
.0768
.0876
.0916
ENTER NEXT COMMAND
?XC
ENTER XFIRST, XLAST, DELTAX ( M)
?,.?,?103.59,163.59,l'5.
ENTER NEXT COMMAND
?RN
A-ll
-------
PLUME2D
VERSION 2.01
PAGE 3
ACCIDENTAL HEXAVALENT CHROMIUM SPILL
CONCENTRATION DISTRIBUTION AT 365.00 DY (MG/L j
* X( M)
Y(
*
M) * .
*
20.00
10.00
0. 00
103.59
.0977
.1115
.1165
118.59
.1036
.1183
.1236
133.59
.1056
.1204
.1260
148.59
.1035
.1181
.1234
163.59
.097.5
.1113
.1163
ENTER NEXT COMMAND
?DN
STOP
A-12
-------
APPENDIX B
Description of Program PLUME2D
Program PLUME2D has been written in an unextended Fortran computer code
in an effort to make the program transportable between computer systems. The
computer code consists of a main program and several function subroutines
which are required to evaluate the Hantush well function. The program has
been documented "internally" through the liberal use of comment statements.
The main program has been divided into three sections. A listing of the
computer code is presented in Appendix D. Section I provides for the "Basic
Input Data" as described in Section II' of this report. The numerical
evaluation of concentration at specified grid coordinates is accomplished in
Section II of the main program which calls subroutine SOL2D,. the code for the
analytical solution of the governing differential equations. Section.Ill
provides for problem redefinition and control of execution under the "Edit"
mode discussed in the body of this report.
Ten function subroutines are used to evaluate the Hantush well function
using the numerical methods described in Appendix C. Listings of the computer
codes are presented in Appendix E. FUNCTION W(U,B) evaluates the Hantush well
function for B < 20. For B > 20, the term EXP(Pex/2) W(U,B) in Equation 45 is
evaluated using FUNCTION WELPRD(U,B,PEX). This procedure is used to avoid
taking the direct product of very large numbers, EXP(Pex/2), and very small
numbers W(U,B), for large values of B.
FUNCTION GAUSS is a 24-point Gauss-Legendre quadrature numerical
integration scheme which is used to evaluate the Hantush well function using
either Equation C-6 or Equation C-7. FUNCTION FUNCTN evaluates the integrand
of Equations C-6 and C-7.
B-l
-------
The six remaining function subroutines are used to evaluate mathematical
functions using rational approximations or polynomial approximations. They
are:
FUNCTION BIO(Z) . Modified. Bessel function of the fi.rst kind
FUNCTION BIOLOG(Z) of order zero and the natural logarithm of
the function.
FUNCTION BKO(Z) Modified Bessel function of the second
FUNCTION BKOLOG(Z) kind of order zero and the natural logarithm
of the function.
FUNCTION EILOG(Z) • Natural logrithm of the exponential
integral.
FUNCTION ERFC(Z) Complimentary error function.
These .six function subroutines are used to support FUNCTION W(U,B) and/or
FUNCTION.WELPRD (U.B.PEX). If system subroutines are available for these
functions they may be substituted for the function subroutines provided with'
Program PLUME2D.
B-2
-------
APPENDIX C
Numerical Evaluation of the Hantush Well Function
The Hantush well function can be defined as
W(U.B) - £i£XP(-e-ld?'. (c-i;
or the reciprocal relation
1 R2
W(U,B) = 2KQ(B) - /2/4uIEXP(- c -f^) df (C-2)
where-c is a dummy integration variable (Hantush, 1964).. Using the identity
/°° fU) d? = /" fU) d? - rQ fU) d? . (C-3;
Equation C-'l can be rewritten as
W(U.B) = /; \ EXP(- 5 - |i) d? - JUQ i (- 5 - |i) dc (C-4)
Now
- JL) dc = 2KQ(B) (C-5)
where KQ- is the modified Bessel function of the second kind of order zero.
Substituting Equation C-5 into Equation C-4, the well function becomes
C-l
-------
W(U,B) = 2KQ(B) - /J jEXP(- C -f|) d? (C-6)
The reciprocal relation, Equation C-2, can also be written in terms of finite
limits. Using the relationship given by Equations C-4 and C-5, the reciprocal
relation can be expressed as
W(U,B) = /* /4U ^EXPf- C -|g) dc (C-7)
For 0 < B < 20, values of W(U,B) for 0 < U < B/2 are obtained from Equation
C-6 by first evaluating the value of the integrand using a 24-poin.t Gauss-
Legendre numerical integration scheme.. For B/2 < U < °°, the reciprocal
relation, Equation C-7', is evaluated using the same numerical integration
scheme.
For 0 < B < 0.1, values of W(U,B) are obtained from the series expansions
presented by Hantush and Jacob (1955). For U < 1
R2
W(U,B) = 2KQ(B) - I
+ EXP(- Jg.) To. 57721566 + ln(U) + EU) + - (1 - ) (C-8)
and for U > 1
W(U,B) = I (B)E (U) - EXP(-U)
01
[(I - -L^} + f^ (^ - -
|_U 36y^ ID tU 4U
(C-9)
where I0 is the modified Bessel function of the first kind of order zero,
is the exponential integral, and 0.57721566 is Euler's constant.
C-2
-------
For B > 20, the third order approximation for W(U,B) presented by Wilson
and Miller (1979) is used to evaluate the well function. The approximation is
W(U,B) - (^)1/2 EXP(-B) [(1 - 4) ERFC (-6) + —L-EXP(-B2)1
& I OB 4B1T1'* J
. (C-10)
where
B =
(4U)1/2
and ERFC is the complimentary error function.
Now, for large positive values of 6,
W(U,B)' * 2 r) EXP(-B) (l.-gf) (C-li;
and an asympotic .expansion for K (B) can be written as
MB) - (^)1/2 EXp(-B) (i - w + —^—9 + - • •) '(c-12)
Thus for B > 20 and B > 7.5 the well function is approximated as
W(U,B) = 2K (B) - (C-13)
Note that this approximation is equivalent to the relationship
W(0,B) = 2KQ(B)
C-3
-------
Evaluations of the Hantush well function using the methods described in
the previous paragraphs have been checked for both accuracy and continuity of
the function between the various approximations. The Gauss-Legendre
quadrature scheme was checked using up to 48 quadrature points. A maximum of
24 quadrature points yielded results accurate to four significant figures in
the mantissa over the entire range of arguments which require numerical
integration. The other approximations for W(U,B) are also accurate to four
significant figures in the mantissa.
C-4
-------
APPENDIX D
Listing of Program PLUME2D
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
PLUME2D
VERSION 2.02
TWO-DIMENSIONAL PLUMES IN UNIFORM GROUND-WATER FLOW
JAN WAGNER
SCHOOL OF CHEMICAL ENGINEERING
OKLAHOMA STATE UNIVERSITY
STILLWATER. OK 74078
PHONE (405) 624-5280
JULY. 1981
REVISIONS:
2.00
2.01
2.02
2.02A
APR 84
24 NOV 84
9 DEC 84
3 MAY 85
DIMENSION TITLE(30).IC(20),XS(10),YS(10),D(3),LBL(2.6),
1 NR(10),IS(4),NP(2),DEL(2).XL(2).XF(2).CON(7),COL(7)
REAL LAMBDA
INTEGER TITLE
COMMON/IO/NI,NO
COMMON/RATE/0(10,12 ) ,T( 10, 12 ) . MT
COMMON/PHYPRO/ALPHA.BETA.DX,LAMBDA,PE,RD,V
DATA IC/'DE'. 'VX','RD' . 'DX', 'DY' .'02' .'PO' ,'OB' , 'XC1 . 'YC
1 'RT' , 'NP' , 'RN' , 'DN' . 'LI ' , ' MU' . .'ST ' . ' CS ' , ' TC ' /
DATA KPR01/'XY'/. KPR02/'XZ'/.. KHARY/'Y'/. KHARZ/'Z'/
DATA NPAGE/1/
DATA KSOL1,KSOL2/'TR'.'SS'/
DATA IS/'R'.'M'.'A'.'D'/
DATA IY/'Y'/
DATA LBL/' ','(C'.' ','ON','
1 ' ',')'/
'ZC'
'TI •
'NU'
'ED'
READ DEVICE : NI
NI»5
NO-6
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).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 MAX IMG
MAXIMG « 20
INITIALIZE PROGRAM FLOW CONTROL VARIABLES
1 IEDIT » 1
KNTL « 1
••• SECTION I -- BASIC INPUT DATA
READ TITLE
WRITE(NO,3)
3 FORMAT(1H1,2X,'ENTER TITLE',/'
READ(NI.S) (TITLE(I). 1-1.30)
5 FORMAT(30A2)
PL2D001
PL2D002
PL2D003
PL2DO04
PL2DO05
PL2DO06
PL2D007
PL2D008
PL2D009
PL2D010
PL2D011
PL2DC12
PL2D013
PL2D014
PL2D015
PL2D016
PL2D017
PL2D018
PL2D019
PL2D020
PL2D021
PL2D022
PL2D023
•PL2D024
PL2D025
PL2D026
PL2D027
PL2D028
PL2D029
PL2D030
PL2D031
PL2D032
PL2D033
PL2D034
PL2DO35
PL2D036
PL2D037
PL2D038
PL2D039
PL2DO4Q
PL2D041
PL2D042
PL2D043
PL2D044
PL2D045
PL2D046
PL2D047
PL2D048
PL2D049
PL2D050
PL2DO51
PL2D052
PL2D053
PL2D054
PL2D055
PL2D056
PL2D057
PL2D058
PL2D059
PL2D060
PL20061
PL2D062
PL2D063
PL2D064
PL2D065
PL2DO66
PL2D067
PL2DO68
PL2DO69
PL2D070
D-l
-------
SELECT VERTICALLY OR HORIZONTALLLY AVERAGED SOLUTION
KFLOW « 3
WRITE(NO,7)
7 FORMAT(3X.'ENTER COORDINATE SYSTEM'./,
16X.-XY FOR VERTICALLY-AVERAGED SOLUTION'./,
26X,'XZ FOR HORIZONTALLY-AVERAGED SOLUTION'./.' ?')
8 READ(NI.S) KFLG
9 FORMATU2)
IF(KFLG.EO.KPROI) KFLOW-1
IF(KFLG.EO.KPR02) KFLOW»2
GO TO (12.12,10). KFLOW
10 WRITE(NO.11)
11 FORMAT(3X,'ERROR IN PROBLEM SELECTION — REENTER'./.' ?')
GO TO B
12 CONTINUE .
KHAR • KHARY
IF(KFLOW.E0.2) KHAR-KHARZ
DEFINE UNITS
WRITE(NO,15)
15 FORMAT(3X,'ENTER UNITS FOR LENGTH (2 CHARACTERS)'./.' ?')
READ(NI,25) IL
25 FORMAT(A2)
WRITE(NO,35)
35 FORMAT(3X.'ENTER UNITS FOR TIME (2 CHARACTERS)'./,' ?')
READ(NI.25) IT
WRITECN0.45)
45 FQRMATox,-ENTER UNITS FOR CONCENTRATION (6 CHARACTERS)'./,' '
READ(NI,26) IM1.IM2.IM3
26 FORMAT(3A2)
ENTER DATA FOR FIRST PROBLEM
SATURATED THICKNESS
IF(KFLOW.EO.1)GO TO 38
30 IMAGE • MAXIMG/2
44 WRITE (NO,46V IL
46 FORMAT(3X,'ENTER SATURATED THICKNESS (0 FOR INFINITE THICKNESS),
1A2./,' ?')
36 READ(NI,37.ERR«44) ST
37 FORMAT(FIO.O)
IF(ST.GT.O.O) GO TO 39
38 IMAGE * 1
ST • 1.OE32
39 CONTINUE
GO TO (50,40O).IEDIT
POROSITY
50 WRITE(NO.SS)
55 FORMAT(3X.'ENTER AQUIFER POROSITY',/.' ?')
READ(NI.56.ERR*50) P
56 FORMAT(FIO.O)
57 IF(P.GT.O.O.AND.P.LT.1.0) GO TO 59
54 WRITE(NO.SS)
58 FORMAT(3X.'POROSITY MUST BE GREATER THAN ZERO'.
1' AND LESS THAN ONE -- REENTER'./,' ?')
READ(NI.56.ERR-54) P
GO TO 57
59 GO TO (6O.40O).IEDIT
SEEPAGE VELOCITY
60 WRITE(N0.65) IL.IT
65 FORMAT(3X.'ENTER SEEPAGE VELOCITY. 'A2.'/',A2 , / . ' ?')
READ(NI.56.ERR-60) V
66 IF(V.GT.O.O) GO TO 69
64 WRITE(N0.67)
67 FORMAT(3X.'SEEPAGE VELOCITY MUST BE GREATER THAN ZERO',
1' -- REENTER',/.' ?')
READ(NI.56.ERR-64) V
GO TO 66
PL2D071
PL2D072
PL2D073
PL2D074
PL2D075
PL2D076
PL2D077
PL2D078
PL2D079
PL2DOBO
PL2D081
PL2D082
PL2D083
PL2D084
PL2D085
PL2D086
PL2D087
PL2D088
PL2D089
PL2D090
PL2D091
PL2D092
PL2D093
PL2D094
PL2D095
PL2D096
PL2D097
) PL2D098
PL2D099
PL2D100
PL2D101
PL2D102
PL2D103
PL2D104
PL2D105
PL2D1O6
PL2D107
' PL2D108
PL2D109
PL2D11O
PL2D111
PL2D1 12
PL2D113
PL2D114
PL2D1 15
PL2D1 16
PL2D1 17
PL2D1 18
PL2D119
PL2D120
PL2D121
PL2D122
PL2D123
PL2D124
PL2D125
PL2D126
PL2D127
PL2D128
PL2D129
PL2D130
PL2D131
PL2D132
PL2D133
PL2D134
PL2D135
PL2D136
PL2D137
PL2D138
PL2D139
PL2D140
D-2
-------
69 GO TO (70,400).IEDIT
RETARDATION COEFFICIENT
70 WRITE(NO,75) -
75 FORMATOX,'ENTER RETARDATION COEFFICIENT',/,' ?')
READ(NI,56,ERR=70) RD.
76 IF(RD.GE.1.0) GO TO 79
74 WRITE(NO,77)
77 FORMAT(3X,'RETARDATION COEFFICIENT MUST BE GREATER THAN OR'.
1' 'EQUAL TO ONE',/,' -- REENTER',/,' ?')
READ(NI,56,ERR«74) RD
GO TO 76
79 GO TO (80,400),IEDIT
X DISPERSION COEFFICIENT
80 WRITE (NO, 81 ) IL.'IT
81 FORMATOX.'ENTER X DISPERSION COEFFICIENT, SO ',A2,
1'/'.A2,/,' ?')
82 READ(NI,56,ERR=80) DX
IF(DX.GT.O.O) GO TO 85
WRITEIN0.83)
83 FORMATOX,'X DISPERSION COEFFICIENT MUST BE GREATER THAN ZERO',
1' -- REENTER'./.' ?' )
GO TO 82
85 GO TO'(86.400).IEDIT
Y OR 2.DISPERSION COEFFICIENT
86 WRITE(NO,87) KHAR.IL.IT • '
87 FORMAT! 3X .'ENTER ' .-A 1 , ' DISPERSION COEFFICIENT, SO '.A2.
1'/' ,A2./, ' ?' )
88 READ(N!,56.ERR=86) DY
IF(DY.GT.'O.O) GO. TO 90
WRITE(N0.89) KHAR
89 FORMATOX .A 1 .' DISPERSION COEFFICI-ENT MUST BE GREATER THAN ZERO'
1' -- REENTER',/.' ?')
GO TO 88 .'
90 GO TO (91.400).IEDIT
FIRST-ORDER DECAY CONSTANT
91 WRITEIN0.95) IT
95 FORMATOX. 'ENTER DECAY CONSTANT ,. 1 /'. A2 ./.' ?')
READINI,56,ERR=91) DECAY
GO TO (1320.400).IEDIT
. SOURCE RATE SCHEDULE.
DEFINE LOCATIONS AND RATES OF SOURCES
INITIALIZE SOURCE/RATE ARRAYS
1320 MAXRT2 = MAXRT .+ 2
DO 1330 I = 1,MAXSOR"
DO 1330 J=1.MAXRT2
0.0
0.0
0(1.J)
T(I.J)
1330 CONTINUE
JFLOW = 3
1340 WRITE(NO,1345)
1345 FORMATOX.'SELECT TRANSIENT OR STEADY-STATE SOLUTION'./.
16X.'TR FOR TRANSIENT SOLUTION',/.
26X.'SS FOR STEADY-STATE SOLUTION',/,' ?')
1350 READ(NI,25) KSOL
IF(KSOL.EO.KSOLI) JFLOW=1
IF(KSOL.EO.KSOL2) JFLOW>»2
GO TO (1370.1370,1360). JFLOW
1360 WRITE(NO.1365)
1365 FORMATOX . 'ERROR IN SELECTION -- REENTER'/.' ?')
GO TO 1350
1370 WRITEtNO,1375) MAXSOR
1375 FORMATOX,'ENTER THE NUMBER OF SOURCES (MAXIMUM OF',13.' )',/.
PL2D141
PL2D142
PL2D143
PL2D144
PL2D145
PL2D146
PL2D147
PL2D148
PL2D149
PL2D150
PL2D151
PL2D152
PL2D153
PL2D154
PL2D155
PL2D156
PL2D157
PL2D158
PL2D159
PL2D160
PL2D161
PL2D162
PL2D163
PL2D164
PL2D1S5
PL2D166
PL2D167
PL2D168
PL2D169
PL2D17O
PL2D171
PL2D172
PL2D173
PL2D17-
PL2D175
•PL2D176
PL2D177
PL2D17S
PL2D179
PL2D130
PL2D181
PL2D-182
PL2D183
PL2D184
PL2D185
PL2D186
PL2D187
PL2D188
PL2D189
PL2D190
PL2D191
PL2D192
PL2D193
PL2D194
PL2D195
PL2D196
PL2D197
PL2D198
PL2D199
PL2D200
PL2D201
PL2D202
PL2D203
PL2D204
PL2D205
PL2D2O6
PL2D207
PL2D208
PL2D209
PL2D210
D-3
-------
1- ?')
1380 READ(NI.1385,ERR=1370) FDUM
1385 FORMAT(FIO.O)
NS»FDUM
IF(NS.GT.O.AND.NS.LE.MAXSOR) GO TO 1400
WRITE(NO,1395) MAXSOR
1395 FORMATOX,'NUMBER OF SOURCES MUST BE GREATER THAN ZERO '
1'AND LESS THAN',13,' -- REENTER',/.' ?')
GO TO 1380
1400 WRITE (NO,1405) IM1.IM2,IMS,IL,IT,IL,IT
1405 FORMATOX,'MASS RATES HAVE 'UNITS OF (',3A2,') (CU ',A2,'/',A2,
1 ' )/'.A2./.3X, 'TIME HAS UNITS OF ',A2./)
DO 1540 1=1,NS
IF(KFLOW.E0.2) GO TO 1414
1406 WRITE(NO,1410) I.IL
1410 FORMATOX.''ENTER x AND Y COORDINATES OF SOURCE',12,
1' C.A2, ')',/,' ?.?')
READ(NI,1425,ERR«1406) XS(I),YS(I)
GO TO 1440 .
1414 WRITE(NO,1415) I.IL . .
1415 FORMATOX,'ENTER X AND Z COORDINATES OF SOURCE',12,
1' (' ,A2, ')'./,' ?.?' )
READfNI . 1425.ERR=1414) XS(I).YS(I )
1425 FORMAT(2F10.0)
1430 IF(YS(I).GE.0.0.AND.YS(I ) .LE.ST) GO TO 1440
1434 WRITE(NO.1435) ST.IL
1435 FORMATOX,'Z-COORDINATE MUST BE GREATER THAN OR EQUAL TO ZERO'
1' AND',/.3X.'LESS THAN OR EQUAL TO SATURATED THICKNESS ('.
2F10.4,A3.')'.V.3X.' -- REENTER'./.' ?')
READINI ,37',ERR=1434) YS(I)
GO. TO 1430
1440 IF(JFLOW.E0.2) GO TO 1530
0( I , 1 ) = 0.0
T(I. 1 ) = 0.0
1450 WRITE(NO.1455) I.MAXRT
1455 FORMATOX. 'ENTER THE NUMBER RATES FOR SOURCE'.12.
1' (MAXIMUM OF'•, 13. ')',/.' ?') •
1460 READfNI.1465.ERR=1450) FDUM
1465 FORMAT!F1O.O)
NRfI)=FDUM
IFfNRf I') .GT.O.ANO.NR( I ) .LE:MAXRT) GO TO 1480
WRITEINO.1475) MAXRT
1475 FORMATOX.'NUMBER OF RATES MUST BE GREATER THAN ZERO AND '.
1'LESS THAN',13.' -- REENTER'./,' ?')
GO TO 146O
1480 CONTINUE
NRT = NRfI)
DO 1520 J»1.NRT
M = J '+ 1
1484 WRITEfNO, 1.485) I , d . T( I . M-1 ) , IT
1485 FORMATOX. 'SOURCE '.12.', RATE '.12.' STARTS AT'.F8. 1 .A3,/ .
1 - 3X.'ENTER MASS RATE AND ENDING TIME ',/,' ?.?')
READfNI,1495,ERR"1484) Q(I,M),T(I,M)
1495 FORMAT(2F10.0)
1500 IFfTfl ,M) .GT.T(I.M-D) GO TO 1510
1504 . WRITEfNO,1505)
1505 FORMATOX.'ENDING TIME MUST BE GREATER THAN STARTING TIME '
1 ' -- REENTER' ./. ' ?' )
READfNI.37,ERR-1504) T(I.M)
GO TO 1500
1510 CONTINUE
1520 CONTINUE
GO TO 1540
1530 WRITEfNO. 1535) I
1535 FORMATOX.'ENTER STEADY-STATE MASS RATE'.I2./.' ?')
READfNI.37,ERR=1530) 0(1,1)
NRfl) = 0
1540 CONTINUE
IF(IEDIT.E0.2.AND.JFLOW.EQ.1.AND.TF.LE.1.OE-06) GO TO 720
123 GO TO (124.400).IEDIT
PL2D211
PL2D212
PL2D213
PL2D214
PL2D215
PL2D216
PL2D217
PL2D218
PL2D219
PL2D220
PL2D221
PL2D222
PL2D223
PL2D224
PL2D225
PL2D226
PL2D227
PL2D228
PL2D229
PL2D230
PL2D231
PL2D232
'PL2D233
PL2D234
PL2D235
PL2D236
PL2D237
PL2D238
PL2D239
PL2D240
PL2D241
PL2D242
PL2Q2-3
PL2D2--
PL2D243
PL2D246
PL2D247
PL2E2-18
PL2D.2J9
PL2D250
PL2D251
PL2D252
PL2D253
PL2D254
PL2D255
PL2D256
PL2D257
PL2D258
PL2D259
PL2D260
PL2D261
PL2D262
PL2D263
PL202S4
PL2D265
PL2D266
PL2D267
PL2D268
PL2D269
PL2D270
PL2D271
PL2D272
PL2D273
PL2D274
PL2D275
PL2D276
PL2D277
PL2D278
PL2D279
PL2D280
D-4
-------
124
125
126
133
134
135
136
137
138
139
140
130
141
142
143
131
144
145
COORDINATES OF THE OBSERVATION POINTS
WRITEfNO,125) IL
FORMATOX,'ENTER XFIRST, XLAST, DELTAX C.A2,')' ,/,' ?,?,?')
READfNI,126,ERR-124) XF(1),XL(1),DELf1)
FORMAT(3F10.0)
DELf1) = ABSfDELf1))
IF(DEL(1).LE.1.OE-06) XL(1)=XF(1)
IF(KNTL.LE.O) GO TO 400
IF(KFLOW.E0.2) GO TO 136
WRITEfNO,135) IL
FORMATOX,'ENTER YFIRST, YLAST, DELTAY (',A2,')',/,' ?.?,?').
READfNI,126,ERR=134) XF(2),XLf2).DELf2)
DEL(2) = ABS(DEL(2))
IF(DELf2).LE.1.OE-06) XL(2)=XF(2)
GO TO 145 '
WRITEfNO,137) IL
FORMATOX,'ENTER ZFIRST, ZLAST, DELTAZ ('.A2,')',/,' ?,?,?').
READfNI,126,ERR=136) XF(2),XL(2).DELf2)
DEL(2) = ABS(DEL(2))
IF(XF(2).GE.0.0.AND.DEL(2).LE.1.OE-06) GO TO 144
IF(XF(2).GE.0.0.AND.XF(2).LE.ST) GO TO 142
WRITEfNO,139)
FORMATOX. 'ZFIRST. MUST BE GREATER THAN OR EQUAL TO ZERO')
IFfIMAGE.GT.-1) WRITE ( NO . 1 40 ) ST.IL
FORMATOX,' AND LESS THAN OR EQUAL TO SATURATED THICKNESS
1 ('.F10.4.A3,')')
WRITEfNO,141)
FORMATOX.' -- REENTER'./,' ?')
READfNI,56,ERR=i30) XF(2)
GO TO 138
IF(XL(2) .GE.0.0.AND.XL(2) .LE.ST) GO TO 145
WRITEfNO-. 143)
FORMATOX. 'ZLAST MUST BE GREATER THAN OR EQUAL TO ZERO-')
IF( IMAGE ..GT . 1 ) WRITE(NO.14O) ST , I L
WRITEfNO.141)
READfNI,56.ERR=131) XL(2)
GO TO 142
XL(2) = XF(2)
GO TO (720.400).IEDIT
OBSERVATION TIMES
720 IFfJFLOW.E0.2) GO TO 770
724 WRITEfNO,725) IT ' '
725 FORMATOX.'ENTER TFIRST, TLAST, DELTAT ('.A2.')',/.' ?.?,?')
730 READ(NI.735,ERR=724) TF.TL.DELT
735 FORMATOF10.0)
. DELT = ABS(OELT)
74O IF(TF.GT.O.O.AND.DELT.LE.1.OE-06C GO TO 760
IF(TF.GT.O.O) GO TO 750
744 WRITEfNO.745)
745 FORMATOX,'TFIRST MUST BE GREATER THAN ZERO -- REENTER'./.' ?'
READfNI,37.ERR=744) TF
GO TO 74O
750 IF(TL.GT.O.O) GO TO 770
754 WRITEfNO.755)
755 FORMATOX.'TLAST MUST BE GREATER THAN ZERO -- REENTER'./.' ?')
READfNI.37,ERR=754) TL '
GO TO 750
760 TL = TF
770 GO TO (146,780),IEDIT
780 IF(JFLOW.E0.2) WRITEfNO,785 )
785 FORMATOX,'TIME IS NOT A PARAMETER IN STEADY-STATE SOLUTION')
GO TO 400
LIST PROBLEM DEFINITION
146 WRITEfNO,147) NPAGE.(TITLEfI),I=1.30)
147 FORMATf1H1./,3X,'PLUME2D'./,3X,'VERSION 2.02',
1/.3X.'PAGE ',I3,///,3X,30A2,///)
NPAGE = NPAGE + 1
PL2D281
PL2D282
PL2D283
PL2D284
PL2D285
PL2D286
PL2D287
PL2D288
PL2D289
PL2D290
PL2D291
PL2D292
PL2D293
PL2D294
PL2D295
PL2D296
PL2D297
PL2D298
PL20299
PL2D300
PL2D301
PL2D302
PL2D303
PL2D3O4
PL2D305
PL2D306
PL2D307
PL2D308
PL2D3O9
PL2D310
PL2D31 1
PL2D312
PL20313
PL2D314
PL2D315
PL2D316
PL2D317
PL20318
PL2D319
PL2D320
PL2D221
PL2D322
PL2D323
PL2D324
. PL2D325
PL2D326
PL2D327
PL2D328
PL2D329
PL2D330
PL2D331
PL2D332
PL2D333
PL2D334
PL2D335
PL2D336
PL2D337
PL2D338
PL2D339
PL2D340
PL2D341
PL2D342
PL2D343
PL2D344
PL2D345
PL2D346
PL2D347
PL2D348
PL2D349
PL2D350
D-5
-------
IFCMAGE.GT. 1 ) WRITE(NO , 148 ) IL.ST
148 FORMAT(1H0.7X,'SATURATED THICKNESS, (',A2.') ',24X
WRITE(NO,149) IL,IT,V,IL,IT,DX,KHAR,IL,IT.DY.P
149 FORMAT(8X,'SEEPAGE VELOCITY. (',A2,'/',A2,') '.25X,
18X,'X DISPERSION COEFFICIENT (' ,A2,'-*2/' ,A2,' )
28X.A1,'. DISPERSION COEFFICIENT ( ' . A2 , ' *«2/ ' . A2 , ' )
38X, 'POROSITY '.42X.F10.4)
WRITE(NO,150) RD,IT,DECAY
150 FORMAT(//,8X,'RETARDATION COEFFICIENT',28X,F10.4./,
18X,'FIRST ORDER DECAY CONSTANT (1/',A2,')'.18X,F10.4)
GO TO (159, 151). OFLOW
151 WRITE(NO. 153) KHAR . IL. IL . IM1 , IM2 , IM3 . l\ , IT..IL
153 FORMAT(//.3X,'STEADY-STATE SOURCE RATES',//,
13X,'SOURCE',6X,'X'.11X.A1,17X,'RATE',/,
25X.'NO',6X,'(',A2.')',8X.'(',A2,')',6X.'(',3A2.
3')(CU ',A2,'/',A2,')/',A2,/)
DO 157 I"1.NS
WRITE (NO, 155). I ,XS( I ) , YS( I) ,0( 1, 1 )
155 FORMAT(5X.I2.F10.2.2X.F10.2.6X.F16.4)
157 CONTINUE
GO TO 171 .
159 WRITE(NO,160) IM1.IM2.IM3,IL,IT,IL.IT,IL.KHAR,IL
160 FORMAT(//,3X.'SOURCE/RATE SCHEDULE ('.3A2,')(CU
1')/'.A2,//,15X,'SOURCE'.13X.'RATE1.4X,'MASS'.8X.'TIME
2/.3X,'NO X l',A2,') ',A1.' (',A2.')'.9X.' t
35X.'START',7X,'END',/)
DO 170 1=1,NS
WRITE(NO,165) I,XS(I ) .YS(I )
165 FORMAT(/'.3X. I2.2F9.2.)
NRT = NR(I )
DO 170 J=1,NRT
M = J * 1 '
' WRITEINO,167) J.0(I.M),T(I.M-1).T(I,M)
167 FORMAT(34X.I 2.F12.2,2F9.2)
170 CONTINUE
1T1 WRITE(NO, 175) IL,XF( 1),XL( 1 ) .DEL( 1 ),KHAR,XF<2),KHAR.XL(2).KHAR,
'1 DEL(2)
175 FORM4T(//.8X.'OBSERVATION POINTS (',A2.')'.//,
112X. 'XFI-RST =' .F10.2.3X. ' XLAST = ' . F 10 . 2 . 3X . ' DELX = ',F10.4.,<
212X,A1,'FIRST ='.F10.2.3X.A1.'LAST = ' .r 10.2,3X, 'DEL' .A 1 , '
1F10.4)
IF(JFLOW.EO.1) WRITEtNO.177) IT.TF.TL,DELT
'177 -FORMAT ( / , 8X , 'OBSERVATION TIMES (', A2 .')',//,
1 12X,'TFIRST ='.F10.2.3X,'TLAST =',F10.2,3X,'DELT' =
180 CONTINUE
GO TO 400
'*- SECTION II -- NUMERICAL EVALUATION OF CONCENTRATION AT
SPECIFIED GRID COORDINATES
NUMBER OF OBSERVATION POINTS IN EACH COORDINATE DIRECTION
2000 CONTINUE
DO 202O L»1.2
NP(L) = 1
OEL(L) = ABSIOEL(L))
IF(DEUL) .LE . 1 .OE-03) GO TO 2020
OIF = XL(L) - XF(L)
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 = OIF - 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
MAXRW = NP(2)
MAXCL « NP(1)
F10.4)
10.4,/,
13X.F10.4,/,
' , 13X.F10.4,/,
)
,2. '/' ,A2.
IE ( ' ,A2, ' ) ' .
5X, 'RATE ' , •
XL(2) .KHAR,
F 1 0 . 4 . / .
. A 1 , ' = ' ,
.F10.4)
IN AT
CTION
PL2D351
PL2D352
PL2D353
PL2D354
PL2D355
PL2D356
PL2D357
PL2D358
PL2D359
PL2D360
PL2D361
PL2D362
PL2D363
PL2D364
PL2D365
PL2D366
PL2D367
PL2D368
PL2D369
PL2D370
PL2D371
PL2D372
PL2D373
PL2D374
PL2D375
PL2D376
PL2D377
PL2D378
PL2D379
PL2C380
PL2D381
PL2C382
PL2D3S3
P'_2D38J
PL2D385
PL2D386
PL2D367
PL2D3S6
PL2D389
PL2D39C
PL2D391
PL2D392
PL2D393
PL2D394
PL2D395
PL2D396
PL2D397
PL2D398
PL2D399
PL2D400
PL2D401
PL2D402
PL2D403
PL2D404
PL2D405
PL2D4Q6
PL2D407
PL2D4Q8
PL2D409
PL2D410
PL2D41 1
PL2D412
PL2D413
PL2D414
PL2D415
PL2D416
PL2D417
PL2D4 18
PL2D419
PL2D420
D-6
-------
TIME COORDINATES
NTIME « 1
IF(DELT.t_E. 1 .OE-06) GO TO 2110
NTIME * ABS(TL-TF)/DELT + 1.0
IF(TF.GT.TL) DELT=-DELT
2110 TSOL = TF
MTIME * NTIME
DAMK = DX*DECAY*RD/(V*V)
ALPHA=SORT(1.0+4.0*DAMK)
PE=V/DX
BETA = DX/DY
LAMBDA = 1.0/(12.566731*P«SORT(DX*DY) )
DO 26SO NT*1,NTIME
2120 LPRT = 1 .
LP » 1
NCFLG = 1
2140 NROW1 = 1
NROW2 = MAXRQW
2160 IFINROW2.GT.MAXRW) NROW2=MAXRW
DO 2580 NROW=NROW1,NROW2
GO TO (2180,2220.2200).NCFLG
2180 NCOL1 = 1
NCOL2 = MAXCOL
22OO !F(NCOL2 .GT.MAXC'L) NCOL2=MAXCL
NCOL = MAXCOL
!F(NCOL2.EO.MAXCL) NCOL=NCOL2-NCOL1+1
2220 1X1 = NCOL1
1X2 = NCOL2
DO 2300 L»1.MAXCOL
CON(L) = 0.0
2300 CONTINUE
DO 2^40 N=!/NS
0(1) = ST - YS(N)
IFIST.GE.0.9E32) D(1)=O.O
0(2) = YS(N)
' DO) = D( 1 )
COEF = 1.0
IF(D(1 ) .LT, 1 .OE-03.0R.D(2).LT.1.OE-03) COEF = 2.0
DO 2440 1=1X1.1X2
X = XF(1)'+ FLOAT(1-1 )«DEL( 1 )
IF(I.EO.NP(1)) X=XL(1)
XXS = X - XS(N)"
PEX = PE-XXS
Y = XF(2) + FLOATINROW-1)«DEL(2)
IF(NROW.EQ.NP(2)) Y=XL(2)
YYS » Y - YS(N)
PEY = PE*YYS
L = I-IX1 + 1 .
IFICON(L).LT.O.O) GO TO 2430
IF(ABS(XXS).LT.1.0.AND.A8S
-------
ZIMAGE' » (2.0-D(LM)+ZM) + 2.0»FLOAT(IM)«D(LM+1)
1 + FLOAT(2*IM-2)«D(LM)
PEZ = PE-ZIMAGE
CALL SOL2D(C,PEX,PEZ,TSOL,N,NR(N))
IF(C.LT.1.0E-06) GO TO 2312
CXYT « CXYT + 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 SOL2D(C,PEX.PEZ,TSOL.N,NR(N)')
IF(C.LT.1.OE-OS) GO TO 2320
CXYT * CXYT + COEF-C
2314 CONTINUE '
WRITE(NO,2315) MAXIMG.X.Y
2315 FORMATOX,<•»-»«« WARNING -- SOLUTION DID NOT',
1 ' CONVERGE USING',/,9X,12,' IMAGE WELLS AT X =',
2 F10.4.' Z •'.F10.4)
2320 CONTINUE
2325 IF(KFLOW.EQ. 1) CXYT*CXYT/2 .0
CON(L) = CON(L) + CXYT
GO TO 2340
2330 CON(L) • -9.9999
2340 ROW = Y
COL(L) = X .
2430 CONTINUE-
2440 CONTINUE
: PRINT CONCENTRATION DISTRIBUTION
GO TO (2460.2560), LPRT
2460' WRITEtNO, 147) NPAGE. (TITLE(I ) . I * 1,30)
. NPAGE = NPAGE + 1
IF(JFLOW.EQ.2) GO TO 2500
WRITE(NO,2465) TSOL,IT,IM1.IM2.IMS,
KLBLtLP.L ) ,L=1 ,6) . IL
2465 FORMAT(13X.'CONCENTRATION DISTRIBUTION AT '.F10.2.
11X.A2,' ('.3A2.') ',,'/. 13X.3X.6A2.//.
2'•'./,' ' X(',A2,')')
GO TO 2520
250O WRITE(N0.2505) IM1.IM2.IM3.(LBL(LP.L ) ,L=1,6 ) ,
1IL
2505 FORMAT(13X.'CONCENTRATION DISTRIBUTION AT STEADY .STATE',
1' C.3A2,') './/13X.3X.6A2,//.
2' *',/.' • X( ' ,A2. ' ).' )
2520 CONTINUE
WRITE(NO,2525). (COLlL).L=1,NCOL)
2525 FORMAT(' •' ,4X,7F10.2 )
WRITE(N0.2545) KHAR.IL
2545 FORMAT(1X,A1,•(',A2,') •',/,9X,'-')
2560 WRITEtNO,2565) ROW,(CON(L),L=1,NCOL)
2565 FORMAT(2X,F8.2.7F1O.4)
LPRT = 2
2580 CONTINUE
IF(NROW2.EO.MAXRW) GO TO 2600
NROW1 3 NROW1 •<• MAXROW
NROW2 - NROW2 + MAXROW
LPRT = 1
LP * 2
NCFLG - 2
GO TO 2160
2600 IF(NCOL2.EO.MAXCL) GO TO 2640
NCOL1 = NCOL1 + MAXCOL
NCOL2 = NCOL2 + MAXCOL
LPRT = 1
LP = 2
NCFLG = 3
PL2D491
PL2D492
PL2D493
PL2D494
PL2D495
PL2D496
PL2D497
PL2D498
PL2D499
PL2D500
PL2D501
PL2D502
PL2D503
PL2D504
PL2D505
PL2D506
PL2D507
PL2D508
PL2D509
PL2D510
PL2D51 1
PL2D512
PL2D513
PL2D514
PL2D515
PL2D516
PL2D517
PL2D518
PL2D51'9
PL2D520
PL2D521
PL2D522
PL2D523
PL2D524.
PL2D525
PL2D52S
PL2D527
PL2D52S
PL2D529
PL2053O
PL2D531
PL2D532
PL2D533
PL2D534
PL2D535
PL2D536
PL2D537
PL2D538
PL2D539
PL2D540
PL2D541
PL2D542
PL2D543
PL2D544
PL2D545
PL2D546
PL2D547
PL2D548
PL2D549
PL2D550
PL2D551
PL2D552
PL2D553
PL2D554
PL2D555
PL2D556
PL2D557
PL2D558
PL2D559
PL2D560
D-8
-------
2640
2660
** • * * I
400
401
405
410
415
r»
420
425
430
GO TO 2140
CONTINUE
TSOL = TSOL * DELT
IF(NT.EO.MTIME) TSOL=TL
CONTINUE
•• SECTION III -- PROBLEM REDEFINITION AND CONTROL OF EXECUTION
CONTINUE
IF(IEDIT.E0.2)GO TO 401
WRITE(NO,1001)KHAR.KHAR,KHAR.KHAR
IEDIT = 2
KNTL = 0
WRITE(N0.405)
FORMAT(//,3X,'ENTER NEXT COMMAND',/.' ?')
READ(NI.415) NEXT
FORMAT(A2)
DO 420 1=1.20
IF(NEXT.EO.ICd) ) GO TO 430
CONTINUE
WRITE(NO,425)
FORMAT(3X.'ERROR IN LAST COMMAND -- REENTER',/.' ?')
GO TO 410
GO TO (91.60.70,80,86.86,50,450.124,133,133,3060,1,
12OOO.700,146,1000,602,1320,720),I
C NEW SET OF X AND Y OBSERVATIONS
J5O KNTL = 1
GO TO • 124
C
C
C NEW .SOURCE/RATE SCHEDULE
3060 WRITE(NO.3065)
3065 FORMAT(3X;'ADD!A).DELETEID).MODIFY(M) A SOURCE OR RETURN!R)'
1: TO EDIT ?' )
3070 READINI.3075) ISK
3075 FORMAT! A', )
DO 3080 1=1.4
!F( ISK.EO.. IS( I ) ) GO TO 3090
3080 CONTINUE
WRITE(N0.3085)
3085 FORMAT(3X,'ERROR IN SELECTION -- REENTER ?')
GO TO 3070
3090 GO TO (400.3100', 3450, 349O) . I
C
C MODIFY SOURCE
C
3100 WRITE(N0.3105) NS
3105 FORMAT(3X.12,' SOURCES IN CURRENT SCHEDULE'./.'
13X.'ENTER SOURCE TO MODIFY'./.' ?')
READINI.1465,ERR=3ioo) FDUM
JS=FDUM
IF(JS.GT.O.AND.JS.LE.NS) GO TO 3220
WRITE(N0.3215) JS
3215 FORMATOX, 'SOURCE' . 14, ' NOT IN SCHEDULE')
GO TO 3060
3220 GO TO (3230,3260).JFLOW
3230 WRITE(N0.3235) JS.XS(JS),IL,KHAR,YS(JS),IL,IT.
1IM1,IM2,IM3.IL,IT.IL
: X =' .F8.2.A3.2X. ', ' ,A1 . '
'MASS RATE' , 14X, 'TIME (' .
A2. '/' ,A2, ' )/',A2,
3235 FORMAT(3X. 'SOURCE ',12,'
1F8.2.A3.//.3X, 'RATE' , 7X ,
2A2, ' )' ,/,4X, 'NO' ,3X, ' ( ' ,3A2. ' )(CU
35X, 'START' ,7X, 'END' ,/)
NRT = NR(JS)
DO 3250 J= 1 ,NRT
M = J + 1
WRITE (NO, 3245) J . 0( JS . M ) . T( JS . M- 1 ) . T( JS .M )
3245 FORMAT (4X, 12 , 5X . F 1 4 . 2 . 7X . F8 . 3 . 3X . F8 . 2 )
3250 CONTINUE
PL2D561
PL2D562
PL2D563
PL2D564
PL2D555
PL2D566
PL2D567
PL2D568
PL2D569
PL2D570
PL2D571
PL2D572
PL2D573
PL2D574
PL2D575
PL2D576
PL2D577
PL2D578
PL2D579
PL2D580
PL2D581
PL2D582
PL2D583
PL2D584
PL2D585
PL2D586
. PL2D587
' PL2D5S8
PL2D589
DL2D59O
PL2D591
PL2D592
' PL2D593
P.L2D5SJ
PL2D595
PL2D596
PL2D597
PL2D598
PL2D599
•PL2D600
PL2D60'
PL2D6C2
PL2D603
PL2D604
PL2D605
PL2D606
PL2D607
'PL206O8
PL2D609
PL2D610
PL2D611
PL2D612
PL2D613
PL2D614
PL2D615
PL2D616
PL2D617
PL2D618
PL2D619
PL2D620
PL2D621
PL2D622
PL2D623
PL2D624
PL2D625
PL2D626
PL2D627
PL2D628
PL2D629
PL2D630
D-9
-------
GO TO 3270
3260 WRITE(NO,32S5) JS,XS(JS).IL,KHAR,YS(JS),IL,0(JS,1) .
1IM1,IM2,IM3,IL.IT,IL
3265 FORMATOX, 'SOURCE ',12,': X =',F8.2,A3.2X, ' ,',A 1 , ' = '
1F8.2.A3./.3X.'STEADY-STATE MASS RATE «'.F16.4.
2' (',3A3,')(CU ',A2.'/'.A2,')/'.A2,/)
3270 WRITE(NO,3275)
3275 FORMATOX, 'CHANGE COORDINATES (Y/N)?')
REAO(NI,3075) JC
IF(JC.NE.IY) GO TO 3290
IF(KFLOW.E0.2) GO TO 3277
3276 WRITE(N0.1410) JS.IL
READ(NI,1425,ERR = 3276X XS(JS ) ,YS(JS)
GO TO 3290
3277 WRITEfNO,1415) JS.IL
READ(NI.1425,ERR=3277) XS(JS).YS(JS)
3280 IF(YS(JS).GE.O.O.AND.YS(JS).LE.ST) GO TO 3290
3281 WRITEfNO.1435) ST.IL
READ(NI,37,ERR»3281) YS(JS)
GO TO 3280
3290 GO TO (3300,3430),JFLOW
C
C TRANSIENT SOURCES
3300 WRITE(NO,3305) JS
3305 FORMATOX, 'MODIFY RATE SCHEDULE FOR SOURCE',13,' (Y/N)
READ(NI,3075) JY
' IF(JY.NE.IY) GO TO 3060
3310 WRITE(N0.3315)
3315 FORMAT!3X.'ENTER RATE TO BE CHANGED'./,
13X.'(ENTER 0 TO CHANGE ALL RATES)',/.' ?')
READlNI,1465,ERR=3310) FDUM
JR=FDUM
IF(JR.LE.O) GO TO 3350
IF(JR.LE.NR(JS)) GO TO 3330
WRITE(NO,3325) JR
3325 FORMATOX, 'RATE
GC TO 3300
3330 WR!TE(N0.3335) JS.JR.T(JS.JR)
3335 FORMATOX. 'SOURCE
.12,' NOT IN CURRENT SCHEDULE')
TdJS. JR*1 ) . IT
SATE '.12.' STARTS AT'.F8.2.
1' AND ENDS AT' , F8 . 2 . A3 . / . 3X, ' ENTER NEW MASS RATE'./.'
M = JR + 1
READ(NI,3345,ERR=3330) O(JS.M)
3345 FORMAT(FIO.O)
GO TO 3300
3350 NRT = NR(JS)
DO 3360 J=1.NRT
M = J + 1
0(JS,M) = 0.0
T(JS,M) = 0.0
336O CONTINUE
3370 WRITEfNO,1455) JS.MAXRT
3380 READ(NI.1465.ERR=3370) FDUM
NR(JS)=FDUM
IF(NR(JS).GT.O.AND.NR(JS).LE.MAXRT) GO TO 3390
WRITE(NO,1475) MAXRT
GO TO 3380
3390 CONTINUE
NRT = NR(JS)
DO 3420 J=1,NRT
M « J + 1
3394 WRITE(NO,1485) JS.J.T(JS,M-1),IT
READ(NI,1495,ERR=3394) 0(JS,M),T(JS,M)
340O IF(T(JS.M),GT.T(JS.M-1)) GO TO 3410
3404 WRITE(NO.1505)
READ(NI,37.ERR=3404) T(JS.M)
GO TO 34OO
3410 CONTINUE
3420 CONTINUE
GO TO 3060
PL2D631
PL2D632
PL2D633
PL2D634
PL2D635
PL2D636
PL2D637
PL2D638
PL2D639
PL2D640
PL2D641
PL2D642
PL2D643
PL2D644
PL2D645
PL2D646
PL2D647
PL2D648
PL2D649
PL2D650
PL2D651
PL2D652
PL2D653
PL2D654
PL2DS55
PL2D656
PL2D657
PL2D658
PL2DS5S
=>L2DSSO
PL2DSS1
PL2D6S2
PL2D6S3
°L2D6S-
PL2DS65
PL2D666
PL2D667
PL2DE63
PL2D669
=L2D67C
PL2D671
•PL2D672
PL2D673
PL2D674
PL2D675
PL2D676
PL2D677
PL2DS78
PL2D679
PL2D68O
PL2D681
PL2D682
PL2D683
PL2D684
PL2D685
PL2D686
PL2DS87
PL2D688
PL2D689
PL2D690
PL2D691
PL2D692
PL2D693
PL2D694
PL2D695
PL2D696
PL2D697
PL2D698
PL2D699
PL2D700
D-10
-------
c
C STEADY-STATE SOURCES
3430 WRITE(NO,3435) JS
3435 FORMAT(3X,'CHANGE STEADY-STATE RATE FOR SOURCE ',12,' (Y/N) ?'
READ(NI,3075) JC
IF(JC.NE.IY) GO TO 3060
3444 WRITE(NO.3445) JS
3445 FORMATOX,'ENTER NEW STEADY-STATE MASS RATE FOR SOURCE ',I2,/,
1- ?')
READ(NI,3345,ERR=3444) 0(JS,1)
GO TO 3060 . • .
ADD A NEW SOURCE
1
3450 NS = NS
OS = NS
IF(KFLOW.E0.2) GO TO 3455
3454 WRITE(NO.1410) JS.IL
READ(NI.1425,ERR=3454) XS(JS),YS(OS)
GO TO 3470
3455 WRITEfNO,1415) dS.IL
READ(NI,1425,ERR=3455) XS(dS ) .YS(dS )
3460 IFfYS(JS).GE.0.0.AND.YS(dS).LE.ST) GO TO 3470 '
3464 WRITEfNO,1435) ST.IL
READ(NI,37,ERR=3464) YS(dS)
GO TO 3460
3470 GO TO (3370.3480),JFLOW
C ' •
C STEADY-STATE SOURCES
348O WRITEINO,3485) JS
3485 FORMATOX, 'ENTER STEADY-STATE MASS RATE FOR SOURCE '.12.
M .' ?'•)
READINI,3345.ERR=3480) OfJS.1)
• NR( JS)- ° 0
'GO TO 3060
C ' • '
C DELETE A SOURCE
C
3490 IF(NS.GT.I) GO TO 35OC
WRITEINO.3495)
3495 FORMATOX,'ONLY ONE SOURCE IN SCHEDULE -- CAN NOT DELETE'./)
GO TO 3060
35OO WRITE(NO.3SO5) IL.IL.IL
3505 FORMATOX. 'SOURCE' ,6X. 'X ( ' , A2 . ' ) ' . 3X . ' Y . ( ' , A2 . ' ) ' . 3X , •
1 ' Z C, A2 ,')',/)
DO 3520 1=1.NS
WRITEINO,3515) I.XS(I ) .YS(I )
3515 FORMAT(5X,12,3X,F8.2.3X.F8.2)
3520 CONTINUE
3530 WRITEfNO.3535)
3535 FORMATOX. 'ENTER SOURCE TO DELETE'./.
13X.'(ENTER 0 .TO CANCEL)'./.' ?')
READ(NI,1465,ERR=3530) FDUM
dS=FDUM
IF(dS.LE.O) GO TO 3060
IF(JS.LE.NS) GO TO 3550
WRITEfNO,3545) dS
3545 FORMATOX,'SOURCE '.12,' NOT IN CURRENT SCHEDULE')
GO TO 3530
3550 WRITEfNO,3555) JS
3555 FORMATOX,'DELETE SOURCE ',12,' (Y/N)?')
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
00 3570 J=JS.NSD
XS(J) = XS(J+1 )
PL2D701
PL2D702
PL2D703
PL2D704
PL2D705
PL2D7O6
PL2D707
PL2D708
PL2D709
PL2D710
PL2D711
PL2D712
PL2D713
PL2D714
PL2D715
PL2D716
PL2D717
PL2D718
PL2D719
PL2D720
PL2D721
PL2D722
PL2D723
PL2D724
PL2D725
PL2D726
PL2D727
PL2D728
PL2D729
PL2D730
PL2D731
PL2D732
PL2D733
PL2D734
PL2D735
PL2D736
PL2D737
PL2D733
PL2D739
=L2D7iQ
PL2D7J1
PL2D.742
PL2D743
PL2D744
PL2D745
PL2D746
PL2D747
PL2D748
PL2D749
PL2D750
PL2D751
PL2D752
PL2D753
PL2D754
PL2D755
PL2D756
PL2D757
PL2D758
PL2D759
PL2D760
PL2D761
PL2D762
PL2D763
PL2D764
PL2D765
PL20766
PL2D767
PL2D768
PL2D769
PL2D770
D-ll
-------
YS(J) - YS(J-M)
NR(d) = NR( J-M )
NRT = NR(d)
DO 3570 K=1,NRT
M = K + 1
0(d,M) = 0(0+1.M)
T(d,M) = T(0+1,M)
3570 CONTINUE
3575 NRT = NR(NS)
00 3580 K=1,NRT
1
0.0
'0.0
M « K +
0(NS,M)
T(NS.M)
3580 CONTINUE
NR(NS) = 0
NS = NSD
GO TO 3060
C
C STEADY-STATE SOURCES
3590 IF(JS.EO.NS) GO TO 3605
DO 3600 U=OS,NSD
0(0,1) = 0(0+1.1)
XS(0) = XS(0+1)
YS(J) = YS(0+1 )
3600 CONTINUE
3605 0(NS.1) = 0.0
NS = NSD
GO TO 3060
r;
C CONFIRM WHETHER SATURATED THICKNESS IS A VARIABLE
602 IF(KFLOW.E0.2)GO TO 30
WRITE(N0.605)
605 FORMATOX.'SATURATED THICKNESS IS NOT A VARIABLE IN'./.
1.3X,'X-Y COORIDINATE SYSTEM (VERTICALLY AVERAGED SOLUTION)')
GO TO 400
C ' .
C ..MENU OF EDIT COMMANDS FOR PLUMES VERSION 2.02
1000 WRITE(NO,1001)KHAR.KHAR.KHAR.KHAR
1001 FORMAT!1H1,/,3X.'MENU OF EDIT COMMANDS'.//.
1' RETARDATION COEFFICIENT RD OBSERVATION POINTS OB'
'2' POROSITY' PO x COORDINATES xc-
3' SEEPAGE VELOCITY VX' .6X,A 1 . ' COORDINATES'.7X,A 1 ' C
4' X DISPERSION COEFFICIENT DX MENU OF COMMANDS MU
52X.A1,' DISPERSION COEFFICIENT D' ,A 1,6X, 'LIST INPUT DATA
. 6/,' DECAY CONSTANT
7' SOURCE RATE SCHEDULE
8' NEW PROBLEM
9' CHANGE SOLUTION/SOURCES
GO TO 400
700 STOP
END
DE ' RUN CALCULATIONS RN
RT DONE DN',
NP SATURATED THICKNESS ST'.
CS OBSERVATION TIMES TC')
/ _
/ ;
' ./.
/'. ••
LI ' .
I'./.
/ ,
'/'.
PL2D771
PL2D772
PL2D773
PL2D774
PL2D775
PL2D776
PL2D777
PL2D778
PL2D779
PL2D780
PL2D781
PL2D782
PL2D783-
PL2D784
PL2D785
PL2D786
PL2D787
PL2D788
PL2D789
PL2D790
PL2D791
PL2D792
PL2D793
PL2D794
PL2D795
PL2D796
PL2D797
PL2D798
PL2D799
PL2D80C
PL2D801
- PL2D802
PL 20 803
PL2D804
°L2D805
PL2D806
PL2D807
PL2D308
PL2D809
PL2D810
PL2D81 1
PL2D812
PL2D813
PL2D814
PL2D815
PL2D816
PL20817
PL2D818
PL2D819
PL2D820
PL2D821
PL2D822
D-12
-------
APPENDIX E
Listing of Utility Function Subroutines
-------
FUNCTION BIO(Z) BIO 001
C JAN WAGNER BIO 002
C SCHOOL OF CHEMICAL ENGINEERING BIO 003
C OKLAHOMA STATE UNIVERSITY . BIO 004
C STILLWATER, OK 74078 BIO 005
C TELEPHONE: (405) 624-5280 BIO 006
C BIO 007
C REVISED: 6 JANUARY 1983 BIO 008
C ' BIO 009
C EVALUATION OF MODIFIED BESSEL FUNCTION OF THE FIRST KIND BIO 010
C OF ORDER ZERO BIO 011'
C POLYNOMIAL APPROXIMATIONS ARE USED FOR IO(Z) ' BIO 012
C SEE SECTION 9.8 OF ABRAMOWITZ AND STEGUN (1966) BIO 013
DIMENSION A(9) ' BIO 014
COMMON/IO/NI.NO BIO 015
DATA A/0.9189385.4.32105045,6.09540829.6.45308739,4.6926023, BIO 016
13.88357842,3.63608323,4.10583047,5.5407023537 . BIO 017
C • BIO 018
IF(Z.LE.O.O) GO TO 200 BIO 019
f = Z/3.75 BIO 020
IF(Z.GT.3.75) GO TO 100 ' BIO 021
T2 = T*T BIO 022
T4 = T2-T2 ' • BIO 023
' T6 = T2-T4 BIO 02-
T8 = T2-T6 BIO 025
T10' » T2-T8 ' BIO 026
T12 =• T2-T10 BIO 027
BIO =' 1.0 + 3.5156229-T2 + 3.0899424-T4 * 1.2067492*T6 BIO 028
1 + 0.2659732-TS + 0.036O7S8-T1O+ C.0045813-T12 ' ' 310 O29
RETURN ' - EIO 030
100 CONTINUE ' - BIO 031
SUM = 0.0 ' BIO 032
DO 150 .1=1.9 ' BIO 033
SIGN- =' (-1 .0)-M 1*1 ) ' BI'O 034
IFd.EO.2) SIGN=1.0' ' BIC C35
ARG = -1.0"A(I) - 0.5«ALOG(Z) - FLOAT ( I - 1 ) - ALOG( T-) BIO 036
SUM = SUM * SIGN-EXP(ARG) ' BIO 037
150 CONTINUE . BIC 038
BIOLOG = ALOG(SUM) * Z EIO 039
BIO = EXP(BIOLOG) - BIO 040
RETURN . -BIO 041
200 CONTINUE • ' • BIC 042
IF(Z.LT.O.O) GO TO 300 .' BIO 043
BIO ='1.0 . ' BIO 044
RETURN " BIO 045
300 WRITEIN0.305) Z . BIC 046
305 FORMAT(6X.'ARGUMENT OF BESSEL FUNCTION IO(Z) IS NEGATIVE',/, BIO 047
16X.'Z = '.E12.S,' -- PROGRAM TERMINATED') EIO 048
STOP BIO O49
END BIO 050
E-l
-------
FUNCTION BIOLOG(Z)
JAN WAGNER
SCHOOL OF CHEMICAL ENGINEERING
OKLAHOMA STATE UNIVERSITY
STILLWATER, OK 74078
TELEPHONE: (405) 624-5280
REVISED:. 6 JANUARY 1983
EVALUATION OF THE NATURAL LOG OF A MODIFIED BESSEL
FUNCTION OF THE FIRST KIND OF ORDER ZERO
POLYNOMIAL APPROXIMATIONS ARE USED FOR IO(Z) '
SEE SECTION 9.8 OF ABRAMOWITZ AND STEGUN (1966)
'DIMENSION A(9)
COMMON/IO/NI.NO
DATA A/0.9189385,4.32105045.6.09540829.6.45308739,4.6926023,
13.88357842.3.63608323., 4. 10583047 , 5 . 540702353/
IF(Z.LE.O.O) GO TO 200
T = Z/3.75
IF(Z.GT.3.75) GO TO 100
T2 =• T-T
T4 = T2-T2
T6 = T2'T4
T8 = T2-T6
T10 = T2-T8 ' .
T12 " T2-T10
BIO = 1.0 + 3.5156229-T2 * 3.0899424-T4 * 1.2067492-T6
1 ' + 0.2659732-T8- * 0 . 0360768-T 10 + 0 . 00458 1 3-T 1 2
BIOLOG = ALOG(BIO)
RETURN
100 CONTINUE
SUM =0.0
DO 150 1=1.9
SIGN=(-1.0)"(I*1)
. IF(I.E0.2) SIGN=1.0
ARC = -1.0-AII) - O.S-ALOG(Z)
SUM = SUM + SIGN-EXP(.ARG)
150 CONTINUE
BIOLOG = ALOG(SUM) + 2
RETURN
200 CONTINUE • .
IF(Z.LT.O.O) GO' TO 300
' BIO « 1 .0
RETURN
'300 WRITE(N0.305) Z
305 FORMAT(6X,'ARGUMENT OF BESSEL FUNCTION IO(Z) IS NEGATIVE'./.
16X,'Z - '.E12.6,' -- PROGRAM TERMINATED')
STOP
END
- FLOAT(1-1 )-ALOG(T)
BIOL001
BIOL002
BIOL003
BIOL004
BIOL005
BIOL006
BIOL007
BIOL008
BIOL009
BIOL010
BIOL011
BIOL012
BIOL013
BIOL014
BIOL015
BIOL016
BIOL017
BIOL018
BIOL019
BIOL020
BIOL021
BIOLC22
BIOL023
BIOL024
BIOL025
BIOL026
BIOL027
BIOL028
SIOL029
BICL03C
BIOL031
BIOL032
EIOL033
3IOLC3-
B-IOL035
BIOL036
BIOL037
BIOL03S
BIOL039
BIOLOdO
BIOL041
BIOL042
BIOL043
BIOL044
BIOL045
BIOL046
BIOL047
BIOL048
BIOL049
BIOL050
E-2
-------
FUNCTION BKO(Z)
JAN WAGNER
SCHOOL OF CHEMICAL ENGINEERING
OKLAHOMA STATE UNIVERSITY
STILLWATER, OK 74078
TELEPHONE: (405) 624-5280
REVISED: 6 JANUARY 1983
EVALUATION OF MODIFIED BESSEL FUNTION OF SECOND KIND
OF ORDER 'ZERO
POLYNOMIAL APPROXIMATIONS ARE USED FOR KO(Z)
SEE' SECTION 9.8 OF ABRAMOWITZ AND STEGUN (1966)
COMMON/IO/NI,NO ' '
IF(Z.LE.O.O) GO TO 200
T = Z/2.0
T2
T4
T6
T-T
T2-T2
T2"T4
100
IF(Z.GT.2.0) GO TO
T8 = T2-T6
T10 = T2-T8
T12 = T2-T10
BKO = -1.0-ALOGfTI-BIO(Z) - 0.57721566
1 + 0.42278420-T2 + 0.23069756-T4 +
0.03488590-T6
> + 0.00262698-T8 + 0.00010750'T1O+ 0.0740E-04-T12
RETURN
CONTINUE . - ''
SUM = ( -'.25331414 - O.O7S32358/T - 0.0218955S/T2
1 - 0.01062446/(T-T2) * 0.005S7872/T4
2 - 0.00251540/(T-T4) + 0.00053208/T6)
BKOLOG = ALOG(SUM) - 2 - 0.5-ALOGI2)
BKO = EXP(BKOLOG)
RETURN "
200 CONTINUE
WRITEfNC.205)•Z
20E FORMAT ( 6X. 'A'RGUMENT OF 2ESSEL FUNCTION KO(Z) IS LESS THAN'.
-. ' OR EQUAL TC ZERO '
STOP
END
,/.6X.'Z = '.E12.6.
'-- PROGRAM TERMINATED' I
BKO 001
BKO 002
BKO 003
BKO 004
BKO 005
BKO 006
BKO 007
BKO 008
BKO 009
BKO 010
BKO 011
BKO 012
BKO 013
'BKO 014
BKO 015
BKO 016
BKO 017
BKO 018
BKO 019
BKO 020
BKO 021
BKO 022
BKO 023
BKO 024
BKO 025
BKO 026
BKO 027
BKO 028
BKO C29'
SKC 03C
BKO C31
BKO 032
5KO 033
BKO 034
BKO 035
BKO 036
BKO 037
BKG 03S
3KC CSS
5KO 04C
BKO 041
E-3
-------
FUNCTION BKOLOG(Z)
JAN WAGNER
SCHOOL OF CHEMICAL ENGINEERING
OKLAHOMA STATE UNIVERSITY
STILLWATER. OK 74078
TELEPHONE: (405) 624-5280
REVISED: 6 JANUARY 1983
NATURAL LOG OF MODIFIED BESSEL FUNTION OF SECOND KIND
OF 'ORDER ZERO
POLYNOMIAL APPROXIMATIONS ARE USED FOR KO(Z)
SEE SECTION 9.8 OF ABRAMOWITZ AND STEGUN (1966)
COMMON/IO/NI.NO
IF(Z.LE.O.O) GO TO 200
T « Z/2.0
T2 » T«T
T4 » T2*T2
T6 <= T.2-T4
IF(Z.GT.2.0) GO TO 100
T8 = T2-T6
T2"T8.
T2-T10
-1 .0-ALOG(T)-BIO(Z) - 0.577215S6
0.42278420"T2 * 0 . 23O69756*T4 +
0.00262698-T8 + 0 .00010750'T 10*
ALOG(BKO)
T12
BKO
0.03488590-T6
0.0740E-04-T12
BKOLOG '
RETURN
100 CONTINUE
SUM = (1
1
2
BKOLOG
RETURN
200 CONTINUE
WRITEIN0.205) Z
205 FORMAT'(6X . • ARGUMENT OF BESSEL FUNCTION KC( 2 ) IS LESS THAN •
25331414 - 0.07832358/T * 0.02189568/T2
- 0.01062446/(T-T2) + 0.00587872/T4
- 0.00251540/IT-T4) - 0.00053208/T6)
ALOG(SUM) - Z - 0.5-ALOG(Z),
1' OR EQUAL TO ZERO'
STOP
END
/.SX.'Z = '.E12.6.' -- PROGRAM TERMINATED'
BKOLOO1
BKOL002
BKOLOC3
BKOL004
BKOL005
BKOLOO6
BKOL007
BKOL008
BKOL009
BKOL010
BKOL011
BKOL012
BKOL013
BKOL014
BKOL015
BKOL016
BKOL017
BKOL018
BKOL019
BKOL020
BKOL021
BKOL022
BKOL023
BKOL024
BKOL025
BKOL026
BKOL027
BKOL028
BKOL029
BKOL030
BKOL031
BKOU032
BXOLG33
BKOLC34
BKOL-035
BKOLC36
.BKOL037
3KOL03S
3KOu03S
3KOL04C
BKOLO^1
E-4
-------
FUNCTION ERFC(Z)
JAN WAGNER
SCHOOL OF CHEMICAL ENGINEERING
OKLAHOMA STATE UNIVERSITY
STILLWATER, OK 74078
TELEPHONE: (405) 624-5280
REVISED: 6 JANUARY 1983
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( 2) •= 1 - ERF(Z)
ERF (-Z) = -ERF (Z)
REAL-8 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/I(1.ODOO' + 7.052307840-02'X + 4.22S20123D-02-(X--2
1 + 9.2705272D-03"(X-"3) + 1 . 520143D-04"-( X--4 )
2 + 2. 765720-04-(X--5) * 4 . 30638D-O5-(X"6 )) •- 16 ).
GO TO .100
FOR x>3 AN "ASYMPTOTIC EXPANSION OF THE COMPLIMENTARY ERROR
FUNCTION IS USED.
50. COEFLG = X-X + OLOG(X) + 0.57236494000
FX - 2.0DOO-X-X
SUM = 1 .ODOO-
TERMO = 1 .ODOO.
DO 60 1=2,50 '
01 = I
TERMI = -TERMOM2.0DOO-OI - 3.0D001/FX
IF(DABS(TERMI).GT.DABS(TERMO)) GO TO 70
SUM = SUM * TERMI
TEST = TERMI/SUM .
IF(ABS(TEST).LT.1.OE-16) GO TO 70
TERMO » TERMI'
60 CONTINUE
WRITE(NO,65)
65 FORMAT(6X.'-•" WARNING -- ASYMPTOTIC EXPANSION FOR ERFC DID NOT'.
1' CONVERGE WITH 50 TERMS IN THE SUMMATION')
70 SUM = DLOG(SUM) - COEFLG
IF(SUM.LT.-72.0000) SUM=-72.ODOO
DERFC * DEXP(SUM)
100 CONTINUE
FOR Z<0. ERFC(-Z) = 2-ERFC(Z)
ERFC = DERFC
IF(Z.LT.O.O) ERFC=2.0DOO-OERFC
RETURN
END
ERFC001
ERFC002
ERFC003
ERFC004
ERFC005
ERFCC06
ERFC007
ERFC008
ERFC009
ERFC010
ERFC011
ERFC012
ERFC013
ERFC014
ERFCC15
ERFC016
ERFC017
ERFCO18
ERFC019
ERFC020
ERFCC21
ERFC022
ERFCC23
ERFCC2-
I ERFCC25
ERFC026
ERFCC27
ERFCC28
ERFCC29
ERFC030
ERFCC31
ERFC032
ERFC033
ERFC03J
' ERFC035
ERFC036
ERFC037
ERFC038
ERFC039
ERFC040
ERFC041
ERFC042
ERFC043
ERFC044
ERFC045
ERFC046
ERFC047
ERFCO48
ERFC049
ERFCC'50
ERFC051
ERFC052
ERFC053
ERFC054
ERFC055
'ERFC056
ERFC057
E-5
-------
FUNCTION E1LOG(Z)
JAN WAGNER
SCHOOL OF CHEMICAL ENGINEERING
OKLAHOMA STATE UNIVERSITY
STILLWATER, OK 74078
TELEPHONE: (405) 624-5280
REVISED: 6 JANUARY 1983
EVALUATION OF THE NATURAL LOG OF THE EXPONENTIAL INTEGRAL
POLYNOMIAL APPROXIMATIONS ARE USED FOR E1(Z)
SEE SECTION 5.1 OF. ABRAMOWITZ AND STEGUN (1966)
COMMON/IO/NI,NO
IF(Z.LE.O.O) GO TO 200
Z2 = Z"Z
Z3 = Z-Z2
IF(Z.GT.1.0) GO TO 100
ARGUMENTS LESS THAN UNITY
E1 = - 0.57721566 + 0.99999193*2
1 * O.05519968-Z3 - 0.00976004-Z2-Z2
2 - ALOG(Z)
E1LOG = ALOG(E1) .
RETURN
100 CONTINUE
ARGUMENTS GREATER THAN UNITY .
E1LOG = ALOG(Z2*Z2 + 8.5733287401"Z3 + 18.0590169730-Z2
1 * 8.6347608925*2 * 0.26777.37343 ) '
2 - ALOG(Z2"22 *• 9.5733223-154-23 + 25.S329561486-22
3 - + 21.0996530827-Z * 3.9584969228 )
4 - ALOG(Z) - Z
RETURN
200 WRITE(NO,2O5) 2
205 FORMAT!6X.'ARGUMENT OF EXPONENTIAL INTEGRAL IS LESS THAN'
1" OR EQUAL TO ZERO',/.6X.'Z = '.E12.6.3X.•
2' TERMINATED' )
STOP
END
- 0.24991055*22
0.00107857-Z2-Z3
-- PROGRAM'
E1LGOO1
E1LG002
E1LG003
E1LG004
E1LG005
E1LG006
E1LG007
E1LG008
E1LG009
E1LG010
E1LG011
E1LGO'12
E1LG013
E1LG014
E1LG015
E1LG016
E1LG017
E1LG018
E1LG019
E1LG020
E1LG021
E1LG022
E1LG023
E1LG02-:
E1LG025
E1LG026
E1LG027
E1LG028
E1LG029
E1LG030
E1LG031
E1LG032
E1LG033
E1LG034
•E1LG035
E1LG036
E1LG037
E1LG038
E-6
-------
FUNCTION FUNCTN(Z) FUNC001
INTEGRAND OF HANTUSH WELL FUNCTION FUNC002
REAL-8 DB.Z.ARG.FUNCTN FUNC003
.COMMON BF FUNC004
DB=BF FUNC005
ARG = DLOG(Z) + Z + DB*DB/(4.ODOO-Z) FUNC006
FUNCTN ' DEXP(-ARG) FUNC007
RETURN FUNCOO8
END FUNC009
E-7
-------
FUNCTION GAUSS(A,B,FUNCTN)
NUMERICAL INTEGRATION BY 24 POINT GAUSS-LEGENDRE QUADRATURE
ZEROS AND WEIGHTING FACTORS ARE FROM TABLE 25.4, P916. OF
ABRAMOWITZ AND STENGUN(1966)
REAL'S A.B.C.D,FUNCTN,GAUSS,SUM, W, Z
DIMENSION Z(12).W(12)
DATA Z/0.064056892862065,0.19.1118867473616.0.31S042679696163.
1 0.433793507626045,0.545421471388839,0.648093651936975,
2 0.740124191578554,0.820001985973902.0.886415527004401.
3 0.938274552002732.0.974728555971309,0.995187219997021/
DATA W/0.127938195346752.0.125837456346828.0.121670472927803,
1 0.115505668053725,0.107444270115965,0.097618652104113,
2 0.086190161531953,0.073346481411080.0.059298584915436,
3 0.044277438817419,0;028531388628933,0.0123412297999S7/
.SET UP INITIAL PARAMETERS
C = (B-A)/2.0DOO
D = (B+A)/2.0DOO
...ACCUMULATE THE
SUM =0.0
00
SUM IN THE 24-POINT FORMULA
J = 1. 12
IF(Z(d).EO.0.0) SUM = SUM +
IF(Z-( J) .NE .0.0) SUM = SUM *
1 + FUNCTN(-Z(J)-C
.CONTINUE
W(J)-FUNCTN(D)
W(J)"(FUNCTN(Z( J)"
+ D))
D)
.MAKE INTERVAL CORRECTION AND 'RETURN
GAUSS = C-SUM •
RETURN
END
GAUS001
GAUS002
GAUS003
GAUS004
GAUS005
GAUS006
GAUS007
GAUS008
GAUS009
GAUS010
GAUS011
GAUS012
GAUS013
QAUS014
GAUS015
GAUS016
GAUS017
GAUS018
GAUS019
GAUS020
GAUS021
GAUS022
GAUS023
GAUS024
GAUSC25
GAUS026
GAUS027
GAUS028
GAUS029
GAUS030
GAUSO31
GAUS032
GAUS033
GAUS03-
GAUS035
GAUS03S
GAUS037
GAUS038
E-8
-------
SUBROUTINE SOL2D(C , PEX , PEY , TSOL , N , NR )
NUMERICAL EVALUATION OF ANALYTICAL SOLUTION
REVISED: 18 APRIL 1984
REAL LAMBDA
COMMON/RATE/0( 10. 12) ,T( 10. 12)
COMMON/PHYPRO/ ALPHA , BETA , DX , LAMBDA , PE . RD , V
PEXY = SORT(PEX«"
MT = NR + 1
IF(MT.GT. 1 ) GO TO
BETA«( PEY"2) )
10
STEADY-STATE SOLUTION
S = 0(N,1)
B = 0.5-PEXY*ALPHA
C » 2.0*LAMBDA«S«EXP(PEX/2.0
GO TO 50
BKOLOG(B))
TRANSIENT SOLUTION
10 C = 0.0
I.F(T(N,MT) .LT.TSOL) MT=MT+1
00 40 K=2,MT
IF(T(N,K-1).GT.TSOL) GO TO 50
S » Q(N.K) - 0(N.K-1)
PINO = V«V»(TSOL-T(N,K-1))/(DX«RD)
TAU = PINJ/(PEXY«PEXY)
U = .0. 25/TAU
.6=0.5-PEXY-ALPHA
IF (B.GT.20.0) GO TO 20
WF = W(U.B)
SUMLOG ='PEX/2.0 * ALOG(WF)
TERM = EXP(SUMLOG')
• GO TO 30 - . •
FOR LARGE VALUES OF B USE THE THIRD ORDER APPROXIMATION
OF WILSON AND MILLER (1979) JOUR. HYDRAULICS DIV.. ASCE,
VOL 105, NO HY12..P 1565.
NOTE: THE TERM EXP(PEX/2) . IS INCLUDED IN THE NUMERICAL
APPROXIMATION FOR W(U,B)'TO AVOID COMPUTATIONAL DIFFICULTIES
IN TAKING THE .PRODUCT EXP(PEX/2)-W(U,B )
20
30
40 CONTINUE
50 RETURN
END
TERM = WELPRD(U,B.PEX)
IFfTERM.LE.1.OE-32) TERM=O.O
C = C + LAMBDA-S-TERM
SOL2001
SOL2002
SOL2003
SOL2004
SOL2005
SOL2006
SOL2007
SOL2008
SOL2009
SOL2010
SOL2011
SOL2012
SOL2013
SOL2014
SOL2015
SOL2016
SOL2017
SOL2018
SOL2019
SOL2020
SOL2021
SOL2022
SOL2023
SOL2024
SOL2025
SOL2026
SOL2027
.SOL2028
SOL2C29
SCL2C30
SOL2C31
SOL2C32
SOL2C33
SOL2034
SOL2035
SOL2036
SOL2037
SOL2038
SOL2039
SOL2C40
SOL2041
SOL2042
SOL2043
SOL2044
SOL2O45
SOL2046
SOL2047
SOL2048
SOL2049
SOL2050
SOL2051
SOL2052
E-9
-------
FUNCTION W(U,B)
JAN WAGNER
SCHOOL OF CHEMICAL ENGINEERING
OKLAHOMA STATE UNIVERSITY
STILLWATER, OK 74078
TELEPHONE: (405) 624-5280
REVISED: 6 JANUARY 1983
EVALUATION OF THE WELL. FUNCTION FOR LEAKY' ARTESIAN AQUIFERS
THIS VERSION HANDLES ARGUMENTS OVER THE ENTIRE RANGE
REAL'S A1,A2,FUNCTN,GAUSS,DZ
EXTERNAL FUNCTN
COMMON BF
BF=B
IF(B.GT.O.-I) GO TO 200
IFIU.GT.1.0) GO TO 100
FOR B < 0.1 USE APPROXIMATIONS PRESENTED BY
HANTUSH, M.S. AND C.E. JACOB (1955)
TRANSACTIONS AMERICAN GEOPHYSICAL UNION,
VOL 36, NO. 1 PP. 95 - 100.
IFfU.LE.1.OE-10) GO TO 50
EQUATION 12. FOR U < '• . 0
E1LOGICON)I
+ E1U * (U«B-6/16.0)-( 1 .0 - U/9.0)
CON = B-S/U.O-U)
TERM1 = 2.0-BKOIB)
TERM2 = EXP(BIOLOGIB)
E1U = EXP(E1LQGIO)I
SUM = 0.57721566 + ALOG(U)
SUMLOG = ALOG(SUM)
TERM3 = EXP(SUMLOG - CON)
W = TERM1 - TERM2 •"• TERM3
RETURN
50 W = 2.0"BKO(B)
RETURN
RECIPROCAL RELATION, EQUATION 17, FOR U > 1.0
100 TERM1 * EXP(BIOLOGfB) + EILOG(U))
SUM = (B-B/4,0)-(1.0/U - 1.07(36.0-U-U))
' 1 + ( (B-B/4.0)T*2)-( 1 ,0/(4.0*U) - 1 .0/U.O-U-U) )
SUMLOG = ALOG(SUM)
TERM2 = EXP(SUMLOG - U)
W = TERM1 - TERM2
RETURN
200 CONTINUE
FOR 0.1 < B < 20.0 USE NUMERICAL INTEGRATION
FOR U < B/2, W(U,B) = 2KO(B)-INT(0--U) FUNCTION
FOR U > B/2 W(U.B) = INT(0--B*-2/4U) FUNCTION
A1 = 0.0
A2 = U
B2 = B/2.0DOO
IF(U.GE.B2) A2=B2*B2/U
DZ = GAUSS(A1,A2,FUNCTN)
Z = DZ
W = 2.0-BKOIB) - Z
IF(U.GE.B2) W=Z
RETURN
WELL001
WELL002
WELL003
WELL004
WELL005
WELL006
WELL007
WELL008
WELL009
WELL010
WELL011
WELL012
WELL013
WELL014
WELL015
WELL016
WELL017
WELL018
WELL019
WELL020
WELL021
WELL022
WELL023
WELL02-;
WELL025
WELL026
WELL027
WELL028
WELL029
WELLC30
WELL031
WELL032
WELL033
WELL03-
WEL'i.035
WELLC36
WELL037
WELL038
WELLC3S
WELLO-0
WELL041
WELL042
WELL043
WELL044
WELL045
WELL046
WELL047
WELL048
WELL049
WELL050
WELL051
WELL052
WELL053
WELL054
WELL055
WELL056
WELL057
WELL058
WELL059
WELL060
WELL061
WELL062
WELL063
WELL064
WELL065
WELL066
WELL067
WELL068
WELL069
WELL070
E-10
-------
END WELL071
E-ll
-------
FUNCTION WELPRD(U.B.PEX)
JAN WAGNER
SCHOOL OF CHEMICAL ENGINEERING
OKLAHOMA STATE UNIVERSITY
STILLWATER, OK 74078
TELEPHONE: (405) 624-5280
REVISED: 6 JANUARY 1983
THIS FUNCTION SUBROUTINE EVALUATES EXP(PEX/2)-W(U,B)
USING THE THIRD-ORDER APPROXIMATION FOR W(U,B) PRESENTED
BY WILSON AND MILLER (1979) JOUR. HYDRAULICS DIV., ASCE,
VOL 105,NO HY12, P 1565.
REAL-8 DI. F2, TERMI,
COMMON/IO/NI.NO
TERMO, SUM, Z
PAR « (B-2.0»U)/((4.0*U)**0.5)
.IF(ABS(PAR).GT.3.0) GO TO 50
TERM1 = (1.0- 1.0/(8.0-B))«ERFC(-PAR)
TERM2 = (PAR/(7.0898154"B))/EXP(PAR«"2)
C W » ((1.5707963/B)"-0.5)-EXP(-B)-(TERM1+TERM2) •
C WELPRD = EXP(PEX/2)-W(U.B)
C = ( ( 1 .5707963/B)"0.5)-EXP(PEX/2 - B ) -(TERM 1+TERM2)
SUMLOG * *LOG(TERM1 + TERM2)
WELOG = 0.5-(0.45158271 - ALOG(B)) * (PEX/2.0 - B) +' SUMLOG
IF(WELOG.LT.-72.0) GO TO 20
WELPRD « EXP(WELOG)
RETURN
20 CONTINUE
WELPRD = 1.OE-32
R'ETURN • •
50 CONTINUE
IFfPAR. LT .0.0) GO TO- 100
C FOR B>20 AND PAR>3.0 W(U,B)=2KO(B) ' '
WELOG = PEX/2.0 + BKOLOG(B) +-0.69314718
WELPRD = EXP(WELOG)
RETURN
1OO CONTINUE
C FOR PAR<-3.0 AN ASYMPTOTIC EXPANSION FOR ERFC(-PAR) -IS UTILIZED
C SEE SECTION .7.1 OF ABRAMOWITZ AND STEGUN (1966)
COEFLG = PEX/2.0 - B - PAR-PAR - ALOG(-PAR) - 0.5-ALOGl2.0"B)
1 + ALOGl1.0-0.1250/B)
IF(COEFLG.LT.-72.0) GO TO 200
Z = -PAR
FZ = 2.0DOO-Z*Z
SUM = 1.ODOO
TERMO = i.ODOO
DO 120 1=2.50
D: = i
TERMI = -TERMO"(2.0DOO"DI-3.0DOO)/FZ
IF(DABS(TERMI).GT,DABS(TERMO)) GO TO 150
SUM » SUM + TERMI
TEST = TERMI/SUM
IF(ABStTEST).LT.1.OE-16) GO TO 150
TERMO = TERMI
120 CONTINUE
WRITE(NO,500)
500 FORMAT(6X,'-•* WARNING -- ASYMPTOTIC APPROXIMATION FOR ERFC IN'
1 6X.' FUNCTION WELPRD DID NOT CONVERGE WITH
2 6X.' 50 TERMS IN THE SUMMATION')
150 SUMLOG = DLOG(SUM)
WELOG = COEFLG + SUMLOG
WELPRD = EXP(WELOG)
RETURN
200 CONTINUE
C FOR LARGE NEGATIVE VALUES OF PAR, ERFC(-PAR) -> 2
WELPRD = o.o
RETURN
END
'WELP001
WELP002
WELP003
WELP004
WELP005
WELP006
WELP007
WELPOO8
WELP009.
WELP010
WELP011
WELP012
WELP01'3
WELP014
WELP015
WELP016
WELP017
WELP018
WELP019
WELP020
WELP021
WELP022
WELP023
WELP024
WELP025
WELP026
WELP027
WELP028
WELP029
WELP030
WELP031
WELP032
WELP033
WELP034
WELPO35
. WELP036
WELP037
WELP038
WELP039
WELPO40
WELP041
WELP042
. WELP043
WELP044
WELP045
WELP046
WELP047
WELP048
WELP049
WELP050
WELP051
WELP052
WELP053
WELP054
WELP055
WELP056
WELP057
WELP058
./, WELP059
,/ WELP060
WELP061
WELP062
WELP063
WELP064
WELPO65
WELP066
WELP067
WELP068
WELP069
WELP070
E-12
------- |