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 (Rv. 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* ^
      OwlAUlJ_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


                  Rdx
                                               (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) - iXP(-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
      NI5
      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) KFLOW2
   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.ERR44) 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,ERR74)  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)
     NSFDUM
     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,ERR1406) 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   J1.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 L1.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*PSORT(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  L1.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.0FLOAT(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.0D(LM)-ZM) + 2.0FLOAT(IM)D(LM+1)
     1                    + FLOAT(2"IM)D(LM)
                  PEZ =  PEZIMAGE
                  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,ERR3281) 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.5ALOG(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  TT
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*LAMBDASEXP(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 = VV(TSOL-T(N,K-1))/(DXRD)
       TAU = PINJ/(PEXYPEXY)
       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 * (UB-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.0U)/((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

-------