MARCH 1989
KINEMATIC MODELING OF MULTIPHASE SOLUTE
TRANSPORT IN THE VADOSE ZONE
by
Randall J. Charbeneau
James W. Weaver
Virginiq J. Smith
Environmental and Water Resources Engineering
Department of Civil Engineering
The University of Texas at Austin
Austin, Texas 78712
CR-813080
Project Officer
John E. Matthews
Extramural Activities and Assistance Division
Robert S. Kerr Environmental Research Laboratory
Ada, Oklahoma 74820
ROBERT S. KERR ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U. S. ENVIRONMENTAL PROTECTION AGENCY
ADA, OKLAHOMA 74820
-------
DISCLAIMER
The information in this document has been funded wholly or in part by the United
States Environmental Protection Agency under CR-813080 to The University of Texas at
Austin. It has been subject to the Agency's peer and administrative review, and it has
been approved for publication as an EPA document
i i
-------
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 implement actions which lead to a compatible balance
between human activities and the ability of natural systems to support and nurture life.
/
v
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 programs to: (a) determine the fate,
transport and transformation rates of pollutants in the soil, the unsaturated 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.
The behavior of organic hazardous waste constituents that enter the soil
environment may be predicted using vadose zone transport and fate models. Such models
normally consider only the presence of water, air and solid phases of the natural media. An
oil phase often exists, however, which may greatly influence the fate and transport
characteristics of the hazardous waste constituents. This report presents two simplified
models which account for fate and transport based on multiphase partitioning between the
air, wastes, oil and solid phases of the waste-soil matrix. This information should prove
useful to those interested in predicting the potential fate and transport of hazardous
constituents in complex vadose zone environments containing an oil phase.
Clinton W. Hall, Director
Robert S. Kerr Environmental
Research Laboratory
111
-------
ABSTRACT
The goal of this research was the development of a computationally efficient simulation
model for multiphase flow of organic hazardous waste constituents in the shallow soil
environment Such a model is appropriate for investigation of fate and transport of organic
chemicals introduced to the soil through spills on the ground surface, leakage from surface
impoundments or underground storage tanks, or land treatment of hazardous wastes.
During the initial phases of a site investigation there usually does not exist sufficient data to
warrant the application of comprehensive though computationally expensive numerical
models, and simplified physically based models which can address the transport of an
organic constituent experiencing volatilization, multiphase partitioning, biodegradation, and
migration may be preferred. Two models based on the kinematic theory of multiphase flow
are developed and presented herein, along with a number of illustrative examples. The
Kinematic Oily Pollutant Transport (KOPT) model assumes steady infiltration of water
based on the expected annual infiltration rate, while the Kinematic Rainfall and Oily
Pollutant Transport (KROPT) model includes transient hydrologic phenomena (evaporation
and infiltration) along with a model for stochastic generation of rainfall. The examples
presented suggest that the KOPT model may be preferred for most applications.
This report was submitted in fulfillment of CR-813080 by The University of Texas at
Austin, under the sponsorship of the U.S. Environmental Protection Agency. This report
covers a period from April, 1986 to September, 1988, and work was completed as of
September, 1988.
IV
-------
CONTENTS
DISCLAIMER ii
FOREWORD iii
ABSTRACT iv
FIGURES vii
TABLES ix
1. INTRODUCTION 1
2. CONCLUSIONS AND RECOMMENDATIONS 3
3. RESEARCH BACKGROUND AND SCOPE 5
4. MULTIPHASE FLOW IN POROUS MEDIA 8
5. THE KINEMATIC MODEL 17
MAJOR ASSUMPTIONS OF THE MODEL 17
DERIVATION AND SOLUTION OF THE MODEL EQUATIONS 19
THE SOLUTION FOR THE DISSOLVED CONSTITUENT 24
INITIAL AND BOUNDARY CONDITIONS 31
NUMERICAL TECHNIQUES 31
6. MODEL IMPLEMENTATION AND RESULTS 37
KINEMATIC OILY POLLUTANT TRANSPORT MODEL (KOPT) 37
KOPT Model Results 44
Example 1: Oil Spill into a Fine Sand 44
Example 2: Gasoline Leak into Silt Loam 47
Example 3: Land Treatment in a Sandy Clay Loam 53
KOPT Model Implementation 58
-------
KINEMATIC RAINFALL AND OILY POLLUTANT TRANSPORT
(KROPT) 61
The Formation of Oil Banks 61
Quasi Steady Oil Banks 62
Unsteady Oil Banks 69
KROPT Model Results 77
Example 4: Oil Spill in Coarse Sand 77
ExampleS: Gasoline Leak into Sandy Clay Loam 83
Example 6: Land Treatment 91
KROPT Model Implementation 97
COMPUTATIONAL PERFORMANCE 99
7. DISCUSSION AND CONCLUSIONS 102
PHYSICAL PHENOMENA OCCURRING IN MULTIPHASE FLOW 103
NUMERICS OF THE METHODS SELECTED FOR THE SOLUTION 104
REFERENCES 105
APPENDIX A KOPT USERS GUIDE 109
APPENDIX B KOPT PROGRAM LISTING 114
VI
-------
FIGURES
1. Typical Two Phase Relative Permeability Curves 11
2. Water, Oil and Air Distribution within a Strongly Water Wet Medium 14
3. Oil Relative Permeability 15
4. Profile for Generalized Solution .- 22
t *'
•
5. Control Volume for Volatilization Model 29
6. Flow Chart of Runge-Kutta-Fehlberg Algorithm 35
7. Oil Spill in a Fine Sand 40
8. Oil Profiles 43
9. Example 1 - Base Characteristic Plane for Oil Spill in Fine Sand, 180 Duration 45
10. Example 1 - Conservative and Decaying Profiles at 60 Days and 180 Days 46
11. Sensitivity of Results to Residual Oil Saturation 48
12. Example 2 - Base Characteristic Plane 49
13. Example 2 - Constituent Mass Balance 50
14. Example 2 - Concentration Profiles 51
15. Example 2 - Oil Profiles 52
16. Example 3 - Base Characteristic Plane 54
17. Example 3 - Concentration Profiles 55
18. Example 3 - Concentration Profiles (continued) 56
19. Example 3 - Volatilization Comparison 57
20. Example 3 - Concentration Histories at 30 and 150 cm Depths 59
21. Quasi-Steady Oil and Water Migration 63
22. Quasi-Steady Oil Profile 64
23. Effective Oil Conductivity and Oil Bank Speeds 67
24. Unsteady Oil and Water Migration 70
25. Hypothetical Unsteady Oil Profiles 71
VII
-------
26. Oil Bank Function 73
27. Profiles for Unsteady Flow Example 75
28. Example 4 - Base Characteristic Plane 78
29. Example 4 - Profiles at 11,46,56, and 72 Hours 80
30. Example 4 - Base Characteristic Plane (continued) 81
31. Example 4 - Profiles at 120 and 140 Hours 82
32. Example 4 - Histories at 50,100, and 200 cm Depths 84
33. Water Relative Permeability in a Sandy day Loam Soil 85
34. Example 5 - Base Characteristic Plane 86
35. Example 5 - Profiles at 20,40, and 60 Days 88
36. Example 5 - Constituent Distribution 89
37. Example 5 - Concentration Profiles, (a) Toluene and (b) Benzene 90
38. Example 6 - Base Characteristic Plane 94
39. Example6-WaterBalance 95
40. Example 6 - Constituent Mass 96
Vlli
-------
TABLES
1. Initial and Boundary Conditions 32
2. Input Data for Examples 1,2 and 3 .41
3. KOPT Programming Synopsis 60
4. Input Data for Steady State Oil Bank Discussion 65
5. Examples 4,5 and 6 Input Data 79
6. KROPT Programming Synopsis 98
7. Solution Performance for Examples 1-6 100
IX
-------
SECTION 1
INTRODUCTION
The focus of this research has been on the fate and transport of organic hazardous
waste constituents in the shallow soil environment Organic hazardous waste constituents
may be introduced into the soil by spills on the ground surface, leakage from surface
impoundments or underground storage tanks, or land treatment of hazardous wastes.
Oftentimes in these circumstances, an "oil" phase is present in addition to the water, air,
and soil of the natural media, and the fate and transport characteristics are determined by the
mobility of the water and oil, as well as by multiphase partitioning of the constituent
between the air, water, oil and soil matrix. The emphasis of the research has been on
development of simplified models for fate and transport of organic contaminants which
may be appropriate for the initial screening or investigation of a site when the available data
base does not warrant the application of more general models.
The primary objectives of this cooperative agreement have been to support
RSKERL in their modeling studies of subsurface contamination from organic wastes,
including hazardous waste land treatment operations, and to investigate the application of
kinematic models for multiphase processes including the flow of water and oil, as well as
multiphase-partitioning, volatilization, and biodegradation.
The Kinematic Oily Pollutant Transport (KOPT) and Kinematic Rainfall and Oily
Pollutant Transport (KROPT) programs are implementations of a kinematic multiphase
transport model that was developed as a part of this research. The KOPT program is a
simplified implementation of the kinematic multiphase transport model, intended to be used
as a screening tool for hydrocarbon spills or near-surface releases. It addresses the
questions of how far an oil release might go into the soil and how soon it might get there.
In the KOPT implementation, steady water infiltration is assumed to occur, and the
specified water saturation is taken as representative of climatic conditions. The KROPT
program is the full implementation of the kinematic multiphase transport model and
includes modeling of transient hydrologic phenomena (evaporation and infiltration), along
with a model for the stochastic generation of rainfall. This model handles multiple
-------
loadings or releases of oily wastes at or near the ground surface, multiple rainfall events,
potential oil migration, sorption, volatilization, and biodegradation.
-------
SECTION 2
CONCLUSIONS AND RECOMMENDATIONS
.: The results of this research found that kinematic theory was successfully applied to
multiphase flow of water and oil. Kinematic models have an advantage over traditional
numeric models in that the kinematic models are well-suited for studying the hyperbolic or
wave behavior of the multiphase flow equations. Inclusion of sharp fronts is natural to the
kinematic models, while finite difference or finite element models encounter difficulties
when derivatives in the governing equations tend to become undefined. The kinematic
model simplifies the governing equations revealing some of the fundamental behavior of
the solution and showing how the oil and water fronts interact Solution of the kinematic
model indicates that the incoming water must displace all of the mobile oil for mass to be
conserved. Another important result is that oil banks are formed when the incoming water
displaces the oil that is present. The main limitation to the kinematic model is that
displacement of oil by water may be unstable depending on the fluid properties of the oil
and the speed of the incoming water.
The two computer codes written to implement two versions of the kinematic model,
KOPT and KROPT, run in a computationally efficient manner, and may be used to
estimate the vertical migration of oil and a dissolved constituent within the soil profile,
including the effects of volatilization and biodegradation. Running these codes showed
that the amount of residual oil saturation strongly affects the total depth over which the the
oil will be predicted to spread.
In light of the success of kinematic models for multiphase flow problems, further
development of the KOPT and KROPT programs. Specifically, the following
recommendations are suggested:
• Generalize KOPT and KROPT models for use in simulation of facilitated
transport. Such an application could look at the use of solvents (such as heptane) for
flushing a soil profile to remove low mobility organic hazardous wastes, and then model
the much more rapid volatilization and degradation of the solvent
-------
• Development and implementation of user-friendly versions of the KOPT and
KROPT programs.
• Expand KROPT program to handle unstable displacements of oil by water,
thereby eliminating the main limitation of the model and enabling the use of randomly
generated rainfalls that may exceed the kinematic rate. This program enhancement would
also handle cases where oil is the displacing fluid, instead of water, such as when the oil
loading rate exceeds the kinematic capacity of an initially nonuniform water profile.
• Expand KROPT to handle layered systems, where the effects of the fluids
crossing into media of differing intrinsic permeabilities and pore size distributions can have
an important effect on their migration.
• Modify KROPT to have the option to use other oil relative permeability models,
instead of the Brooks-Corey approach, such as proposed by Parker et al. (1987) and/or
Stone (1970 and 1973).
• Modify KROPT to handle depth and time variability in oil residual saturations.
-------
SECTION 3
RESEARCH BACKGROUND AND SCOPE
The objective of this research was to apply kinematic modeling theory to fate and
transport of hydrocarbons in the vadose zone for oil spills and land treatment sites. The
purpose of the work was to provide a way to estimate the gross movement of pollutants in
a multiphase system resulting from spills, leaks and at land treatment sites. In the case of
land treatment systems the long term behavior of the site is of interest, so computationally
efficient simulations are needed. This was the reason for using a simplified approach.
When placed in porous media, immiscible pollutants retain their unique properties.
Although these pollutants are largely immiscible in water, they are capable of causing
groundwater contamination by their dissolved constituents. Immiscibility leads to distinct
bodies of pollutant in close contact with both air and water, where present. Above the
water table, the addition of an immiscible pollutant causes the original two-phase water/air
system to become a three-phase water/pollutant/air system. When flowing, the transport
and fate of immiscible fluids are governed by a set of relationships which are, basically,
expanded forms of the single phase governing equations. The former provides the
theoretical framework for solving immiscible (or multiphase) flow problems in general.
The specific application of the theory depends on the initial and boundary conditions
imposed during pollutant migration as well as other features unique to pollution incidents
(Weaver, 1984).
MULTIPHASE FLOW MODELS
Currently available analytic and semi-analytic multiphase flow models share the
limitation that they rely on an assumed shape of the oil profile that is not based upon
solution of the governing physically based equations (Dracos, 1978, Mull, 1970 and 1978,
and Raisbeck and Mohtadi, 1974). After an oil release ends, the oil is gradually replaced
by air as it drains. A drainage profile, similar to those observed for water is expected,
since the same forces that act on the water also act on the oil. Other phenomena outside the
scope of these models are biodegradation which also causes variable oil saturations, and
-------
simultaneous unsteady flow of the water. The latter is typical of actual dine-varying
conditions in the field.
Only a few numeric solutions of the multiphase flow equations have been presented
for the specific problem of shallow aquifer contamination (Abriola and Finder, 1985a and
1985b, Corapcioglu and Baehr, 1987, Faust, 1985, Kuppusamy et aL, 1987, Hochmuth
and Sunada, 1985, and Osbourne and Sykes, 1986). A limitation of the numeric approach
t •'
is that the models are computationally very intensive. To exploit fully their capacity for
modeling geologic variability, a large amount of data is required Because data is usually
sparse, the full capacity of a numeric model may never be reached when applied to specific
site investigations.
KINEMATIC THEORY APPLIED TO MULTIPHASE FLOW
The analytic and semi-analytic flow models illustrate that if a constant amount of
water is present in the pore space, only the mass conservation equation for the oil needs to
be solved. A solution, that has not been used for this purpose, is the kinematic model for
soil moisture transport (Sisson et al., 1981, Smith, 1983, and Charbeneau, 1984). In
addition to a solution for the duration of water application, this model gives a solution for
the unsteady flow occurring in soils undergoing drainage (replacement of soil moisture by
air). The solution is obtained by solving an approximate governing equation for the water
by the method of characteristics (MOC). With the proper constitutive relationships for two-
phase flow, the water solution is analytic. Adapting this solution for the oil phase gives a
more realistic model than the analytic and semi-analytic models mentioned above, since the
oil drainage profile is determined by solution of the governing equation. An analytic
solution is not possible for the oil because the relative permeability function is too
complicated and depends on the amount of water in the pore space. However, a semi-
analytic solution for the oil phase can be achieved, allowing relatively low computation
times to be achieved.
The following assumptions are used for the model formation. The model equations
are solved for one-dimensional flow. This is one of the necessary assumptions to extend
kinematic theory to multiphase flow. Lateral migration, however, reduces the amount of
pollutants and water that move downward. The one-dimensional formulation is
conservative as it will determine the maximum amount of substances that will reach a given
depth. A two-component pollutant is assumed. The oil acts as an inert carrier for the
-------
dissolved pollutant. Various transformations, other than biodegradation, that the oil can
undergo are not modeled. Advective transport is used to model the transport of the
constituent by oil and water. Equilibrium, linear partitioning relationships are used for
sorption, dissolution into the water and for volatilization into the air phase. Volatilization is
modeled by applying the wet zone/dry zone model of Thibodeaux and Wang (1982) as
modified by Short (1986). To extend kinematic modeling to the oil phase three more
assumptions are used. First the flow of the air does not impede the flow of the liquids.
Second the liquid phase transport is dominated by gravity. Lastly, the medium is assumed
to be strongly wetted so that the fluid saturation distributions maintain the size ordering
shown in Figure 2. These assumptions are discussed in Sections 4 and 5 below.
Under the conditions to be discussed, the kinematic model can be applied 1) to the
unsteady flow of both the oil and the water and 2) to the fate and transport of the dissolved
constituent. If only the advective motion of the constituent is modeled, the resulting
equations are compatible with the numerical approach taken for the oil equations.
Two separate implementations were created. The first (KOPT) was of a kinematic
model for the oil migration only, with a constant amount of water present in the pore space.
This model is analogous to the simplified analytic and semi-analytic flow models noted
above, but with the capability to simulate degrading oil and the transport of the dissolved
constituent In related research this model was coupled with a model of oil spread along
the water table, which in turn supplies input for a model of solute transport in the saturated
zone. Here the overall model was intended as a screening tool for oil spills. A Users
Guide for the KOPT model is presented in Appendix A while the listing of the code is
presented in Appendix B.
The second implementation (KROPT) was a model to investigate the interactions
between incoming rainfall and oil moving within the profile. This model is suited for
studying the effect of a series of rainfall events, and irrigations and cultivations on a land
treatment site. There are practical limitations to application of the KROPT model and it is
considered a research code at this time.
-------
SECTION 4
MULTIPHASE FLOW PRINCIPLES FOR POROUS MEDIA
The equations that govern multiphase flow in porous media are coupled, nonlinear,
partial differential equations (PDEs). There is a mass conservation equation for each fluid
and for each dissolved constituent in the system. There are also auxiliary relationships
which incorporate various physical phenomena into the PDEs. By nature, these latter
relations are empirical, since the phenomena are highly complex and thus are not suited to
exact mathematical expression. Taken together these equations form a mathematical model
of the flow and associated solute transport (e.g., Peaceman, 1977, pp. 24-27).
A typical mass conservation equation may be written in one dimension as
fr-"-
where m is the mass per bulk volume of porous media, Jm is the mass flux, and Pm is the
source strength. Throughout the rest of this work, the liquids will be assumed
incompressible and the medium non-deformable. The resulting continuity (or phase
conservation) equation is
where qi is the volume flux of fluid i, Si is percent of the pore space filled by fluid i, and T|
is the porosity of the soil. • Sj is the called the degree of saturation of fluid i, or more
commonly, simply the saturation. In a multiphase flow system,
(3)
-------
where i represents each fluid present. In order to model the liquids with the phase
conservation equations 2, the concentration of the dissolved constituent is assumed to have
no effect upon the fluid's migration. Fluids for which partitioning phenomena can
significantly alter the saturation during the period of simulation are not suited to this
approach. Examples are the volatilization of benzene or dissolution of perchloroethene
(PCE) into water.
The flux term of equation 2 can be related to fundamental properties of the system
by the usage of Darcy type equations. Specifically,
(4)
where z is directed positive downward, k is the intrinsic permeability of the media, kn is
the relative permeability of the media to fluid i, m is the dynamic viscosity, pi the density,
g is the acceleration due to gravity, Pi is the fluid pressure, and S stands for the set of fluid
saturations. The quantity Kei(S) is called the effective hydraulic conductivity of the media
to fluid i and is defined by
Pi g „. . ,_.
(•>)
W
where Kg is the fully saturated hydraulic conductivity to i. This definition underscores that
the intrinsic permeability, density and viscosity play an important role in determining the
flux rates in multiphase flow systems.
Equation 4 is an extended form of Darcy's Law. It incorporates the physical
phenomena occurring in a porous medium containing immiscible fluids into a form that can
be used in the continuity equations. Knowledge of the geometry of each individual pore
and of the fluid interactions within the matrix is not needed. This information is contained
in an integrated sense in the relative permeability and capillary pressure functions. The
dependence of these functions on the other fluids results in the coupling of the conservation
equations.
In multiphase systems the fluid pressures are related by the capillary pressure, Pc.
This fundamental relationship is defined for each pair of fluids present by
-------
Pc = Pnw-Pw (6)
where the pressure, Pw, is the pressure in the wetting phase - the fluid most strongly
attracted to the solid. Pnw is the pressure in the nonwetting fluid. For simple geometries
(i.e., single circular tubes), there are analytic expressions for the capillary pressure (see
e.g., Bear, 1972, pg. 446). Soils, however, are composed of a variety of irregular,
tortuous pores. The relationship expressed by equation 6 still holds, but must be expanded
to represent the variation in Pc which can occur over a representative elementary volume.
In contrast to uniform geometries, different sized pores are filled by different fluids, and
their contributions to the overall pressure difference varies.
Numerous characteristics of both the capillary pressure and relative permeability
functions can be related to the underlying physical phenomena (see Weaver, 1984, pp. 57-
60, 62-76). Interpretation of directly measured capillary curves reveals information
concerning these phenomena. Notably, capillary pressure (and relative permeability) are
subject to hysteresis. That is, the capillary pressure depends on the displacement direction,
being higher during drainage (displacement of a wetting fluid by a nonwetting fluid) than
during imbibition (the opposite displacement).
The relative permeability function, kn(S) is the fraction of the fully saturated
hydraulic conductivity existing for the set of saturations S. Numerically it is a number
from 0 to 1. It quantifies the notion that the presence of other fluids reduces the hydraulic
conductivity of the medium to a given fluid. Physically, this occurs due to reduction in the
numbers of pores available to the fluid and due to the increased tortuousity of those
remaining. These effects are illustrated by the typical two phase relative permeability
curves in Figure 1. The permeabilities drop off rapidly in a characteristic manner when the
saturation is reduced. At the trapped or residual saturations Sir, the relative permeability is
zero, implying that the fluid is not present as a continuous phase. For wetting fluids this is
as fluid held in the smallest pores, while residual nonwetting fluids are found in larger
pores as isolated droplets.
Due to the considerable difficulties of measurement, three-phase relative
permeability is rarely measured (see Stone, 1970). Alternatively, various models of
relative permeability are used. One way to develop these is to begin with conceptual
10
-------
Saturation
Figure 1: Typical Two-Phase Relative Permeability Curves
models of the porous medium, presumed distributions of fluids within the medium, and a
solution of laminar flow through the medium (see Bear, 1972, pg. 463). The Burdine
(Burdine, 1953, Wyllie and Gardner, 1958) equations form one such model. For each
phase in the system there is a relative permeability equation:
"VI
/
/
dS_
P2
dS.
Pc2
-------
tj s.
f dS
Pc2
kro = ['Vr.r^r -,— (8)
^
Pc2
Originally, the three-phase integration of the Burdine equations contained no provision for
a residual oil saturation (Wyllie and Gardner, 1958). Equations 7 and 8 have been
modified to fit the requirements of the definition of relative permeability: that kro equals
one only when the pore space is completely filled with oil and that k^ equals zero when the
oil is at its residual saturation. The modification is to the squared term of equation 8, which
represents the ratio of the tortuosities at full saturation to that at any saturation, S0 (Brooks
and Corey, 1964). At S0 = S0r the tortuosity is infinite, since the oil is found as
disconnected blobs in the pore space. This condition gives zero relative permeability.
When the oil fills the pore space (meaning that there is no residual water saturation), the
tortuosity is equal to the saturated tortuosity.
The capillary pressure as a function of saturation is needed to evaluate the integrals
in equations 7 and 8. Brooks and Corey (1964) proposed a power law model for drainage
capillary pressure,
where Se is the effective water saturation, P^ is the minimum capillary pressure required to
force a non-wetting fluid into a fully saturated medium (the bubbling pressure), and X is
indicative of the pore size distribution. Using equation 9 in equations 7 and 8 results in the
following model of drainage relative permeability for water and oil:
12
-------
2+3X ^ . . , . . ,
where e = . Despite its apparent analytic nature, certain parameters must always be
A,
measured to fit the model to a specific situation. For equations 10 and 11 these are the
residual water and oil saturations, and the Brooks and Corey exponent. The values of SWT.
X, and Pce are obtainable through a procedure outlined by Brooks and Corey (1964).
They are essentially the result of a capillary pressure curve measurement An extensive
tabulation of these parameters for many different soil types has been developed, making
this an attractive model to use (Brakensiek et ah, 1981, and Rawls et al., 1983). The value
of the residual oil saturation is more difficult to obtain as less work has been done on it
Parker and Lenhard (1987) have recently published methods for estimating residual oil and
air saturations. In the present work, however, a constant value of SOT is assumed
Equation 10 is the same as would result from a two-phase integration of the
Burdine equations. This is because in a strongly wetted system the strongest wetting phase
(say water) is most strongly attracted to the solid. The most nonwetting fluid (air) is
repelled. Water, then, is found in the smallest pores and grain contacts; air in the largest
pores. Figure 2 illustrates this fluid distribution, and how oil acts as an intermediate phase.
Visual studies confirm this distribution for some systems (see e.g., Schwille, 1988).
Figure 3 shows oil relative permeabilities from equation 11. As required, the
relative permeability at the residual saturation, Sor> is zero. This can be verified by
substituting Sor into equation 11. At the other endpoint, 1 - Swr - Sor, the hydraulic
conductivity is somewhat less than one. This is also a required end point condition,
because by definition, kro = 1 only when the pore space is completely filled with oil (see
equation 5). Figure 3 resembles measured curves presented by Leverett and Lewis (1941).
Between the endpoints, the relative permeability at a given oil saturation increases
with increasing water saturation. This is due to the fluid distribution illustrated in Figure 2.
If the water saturation increases from that shown in Figure 2a to that in 2b, the oil is forced
into larger pores. The larger pores have a lower surface area to volume ratio, which results
in less resistance to flow and greater relative permeability. In the model equation this
behavior is simulated by the limits of the integral in equations 7 and 8 being based on the
fluid distribution.
A comparison with Stone's (1970) relative permeability model, which is based on
pore blocking, showed general agreement with these values. At low water and oil
13
-------
frequency
size
(a)
•frequency
size
(b)
Figure 2: Water, Oil and Air Distribution within a Strongly Water Wet
Medium; Water Saturation in (b) is greater than Water Saturation in (a)
14
-------
1.0
0.9--
0.8--
0.7--
0.6--
Kro(So,Sw) 0.5--
Sw
Sw
Sw
Sw
Sw
0.17
0.30
0.40
0.50
0.60
- Swr
Figure 3: Oil Relative Permeability
15
-------
saturations, however, increasing the water saturation reduces the oil relative permeability
due to pore blocking. Throughout this region, the agreement of the models was poor,
apparently because the qualitative nature of the models is different at low saturations. Both
models indicate very low values of kro at these saturations. The comparison indicated that,
while the Burdine equations model certain phenomena involved in oil phase permeability
(relative tortousity and resistance of pores occupied), they do not necessarily include all the
phemonema. In their review, Saraf and McCaffery (1981) recommended using either the
Burdine equations or Stone's (1970) model. Delshad and Pope (1987) reviewed a variety
of kro models with regard to available data. The latter review showed that different models
fit different sets of measured relative permeability data best This indicates that three-phase
relative permeability modeling is an area that needs further development The three-phase
Burdine equations were selected for this work, since their behavior is very regular.
In most cases of interest the oil applied at the ground surface will contain dissolved
constituents. These may partition between the water, oil, air, and soil phases. A mass
conservation equation for these may be written by using equation 1 and appropriate
expressions for the mass of the constituent per unit volume and the mass flux. The total
mass of constituent per unit bulk volume, me, is
me = n (S0Co + Sw Cw + SaCa) + pbCs (12)
where T\ is porosity, Q>, Cw, and Q are the constituent concentrations in oil, water, and
air, pb is the soil's bulk density, and Cs is the sorbed concentration expressed as mass of
constituent/mass of soil. The mass flux is
Jc = qoQ> + qwCw + qaCa - D0 - DW - Da (13)
dz dz dz
where D0, Dw, and Da are the dispersion coefficients for the constituent in the oil, water,
and air phases. Using partitioning relationships of the form
Cj=fj(Cw) (14)
where the subscript j refers to oil, air or soil, allows the constituent mass conservation
equation to be written in terms of one unknown concentration, say, Cw.
16
-------
SECTION 5
THE KINEMATIC MODEL
Three physically based assumptions give rise to the kinematic model. They are of
central importance to the model, and will be discussed first in this section. Following the
discussion of the assumptions, the derivation of the model equations and their solution will
be presented.
MAJOR ASSUMPTIONS OF THE MODEL
The first assumption of the model is to follow Richards (1931) and neglect
resistance to air flow in the mass conservation equation for water. This is in accordance
with common soil science practice and many of the other multiphase infiltration models.
Richards' assumption means that the flow of the air does not affect the flow of the liquid,
and this assumption is usable as long as pressure does not build up in the air phase
(Youngs and Peck, 1964). The presence of air is accounted for via the usage of a three
phase relative permeability function, and by assuming that there is a residual air saturation.
Experience shows that the maximum effective hydraulic conductivity to water is 40% to
60% of Kws (Bouwer, 1966). For the modeling, 50% Kws is used to determine the
trapped air saturation. Typically this is a fairly small value (< 0.05), since the relative
permeability drops off rapidly.
Second, a gravity drainage or kinematic model is used. Several comments on the
kinematic assumption are warranted. The theory of kinematic waves was presented by
Lighthill and Whitham (1955) for flood waves in rivers. They noted that whenever an
approximate functional relation exists at each point between the flux and concentration
(saturation) then the wave motion follows directly from the equation of continuity. For
most cases of interest the functional relation is nonlinear, and the kinematic waves are either
self-spreading, self-sharpening, or regions of constant state. For multiphase flow the self-
17
-------
spreading wave corresponds to a region of internal drainage and redistribution. Here, the
saturation gradients decrease with time as does the role played by pressure gradients. The
self-sharpening wave corresponds to a wetting front While capillary pressure gradients
determine the actual shape of a wetting front, the kinematic model is able to move the
sharp-front representation of the wetting front downward at the correct speed so that mass
is conserved. For these reasons, kinematic models are able to represent the essential
features of nonlinear wave phenomena.
In addition to the comments concerning the ability of kinematic models to represent
the important features of soil water and hydrocarbon waves, one must distinguish between
the capillary pressure in individual pores and the macroscopic capillary pressures appearing
in the flux equation. Two ways of representing the underlying phenomena on a
macroscopic scale are with the relative permeability and capillary pressure functions.
Research done in studying the effects of reducing interfacial tensions (Bardon and
Longeron, 1980) shows that in the absence of interfacial tension (i.e., pore level capillary
pressures), the relative permeability is equal to the saturation, otherwise kr is less than S.
Using a kr(S) function similar to that appearing in Figure 1 assures inclusion of capillary
phenomena in the model. The same point is illustrated indirectly by equations 7 and 8,
which model kr as a function of the macroscopic capillary pressure function.
For a hydrocarbon release, initially, strong capillary gradients drive the flow.
These are included through use of a Green and Ampt (1911) type of model for both the
water and hydrocarbon phase with a corresponding wetting front suction pressure.
Following the initial wetting period there is a transition period, which leads later to gravity
dominated flow (Philip 1969, pg. 280). Under the kinematic assumption only the term
involving the gravity force is retained explicitly in the liquid flux equations (equation 4),
while the effects associated with capillary pressures are approximately handled implicitly
within the model
The third assumption for the model is that in a strongly wetted system, the relative
permeability to the wetting phase is independent of the other fluids present. The reason for
this is that the most wetting phase is so strongly attracted to the medium that it occupies the
smallest pores no matter what other fluids are present (see Figure 2). Payers and Matthews
(1984) reviewed experimental results extensively and concluded that most researchers
would agree with this assumption in strongly wetted systems. Note that most of the
commonly used models of relative permeability (Burdine, Stone, Parker et al, etc.) are
18
-------
predicated upon the assumption of a distribution of fluids within the medium that would
exist in a strongly wetted system. Consequently, relative permeability models for the
strongest wetting phase are independent of the other saturations.
DERIVATION AND SOLUTION OF THE MODEL EQUATIONS
t
The derivation of the model proceeds from the governing equations and
assumptions as follows. The water flux, qw, is given as the gravity component of equation
4
qw = Kew(Sw) (15)
Next the continuity equation for water can be written in terms of its unknown saturation,
Sw, by expanding the z derivative. The resulting first order hyperbolic equation is
+ w = 0
at dsw az
where its source term is zero because it is assumed for the present argument that water is
conserved within the profile (later the question of water loss through evaporation will be
addressed). In equation 16 all quantities are expressed in terms of kinematic variables
Qength and time only), hence the name kinematic model
After application of the kinematic assumption, the oil equation can be written in
terms of the unknown saturation, SQ. The resulting first order hyperbolic equation is
oKeo(So,Sw) Sw _ .
dt dSo dz dSw dz
where the first order decay term has been used for the source term to account for the
possible biodegradation of the oil since oil might not be conserved; X0 is the decay
constant The third term on the left hand side results from the dependence of the oil relative
permeability on Sw. As noted below, for both water and oil the drainage relative
permeability functions are used (equations 10 and 11).
19
-------
The simplified model equations (16 and 17) can be written in matrix form as
£*A.jS-B (.8)
dt 9z
with
[Swj
ISoJ
A = -
dS
w
0
B =
The system of equations is hyperbolic if the eigenvalues of the matrix A are real and
distinct (Whitham, 1974 and Lax, 1973). The eigenvalues of A are easily determined as
they are identical to the diagonal entries (Kreyszig, 1972, page 277) divided by f\. These
are always real, since they are the derivatives of the real functions kw and km (equations
10 and 11). The values are also distinct for cases of practical interest
The consequences of strong wetting on kinematic flow are apparent in matrix A.
For water in a strongly water wet system, the kinematic flux (relative permeability) is
independent of the other fluids. This results in the partial derivative of the water flux with
respect to oil saturation S0, being equal to zero. The mass conservation equation for the
water is independent of the other equations, and can be solved first This is true only as
long as the flow remains kinematic, and the "0" element remains in matrix A.
Since the system of equations is hyperbolic, the method of characteristics can be
used to determine the solution. Where the water saturation distribution varies
continuously, the classical method of characteristics solution of equation 16 is
dz i = 1 dKeW(Sw)
dt 'Sw -n dSw
(19)
20
-------
Equation 19 is valid only where the derivatives appearing in equation 16 exist at each point
of the solution domain (see Rozdestvenskii and Janenko, 1980). When the derivative,
^^g , is a known function, equation 19 can be integrated analytically, leading directly
Qo\v
to an expression of z as a function of time. Then for a given saturation and time, the
associated depth is known immediately. This is the case with relative permeability equation
10 (Charbeneau, 1984). /'
>'
After obtaining the solution of the water equation, equation 17 for the oil phase may
also be solved by the method of characteristics. Application of the method of
characteristics gives the result*
dS0 1 9Keo(So,Sw) 9SW
~dt = " -- ^ -- ^T~
m TJ 3SW 9z
along
dz _ 1_ 9Keo(S0.Sw)
" ( }
Equation 21 is the projection of the characteristic path onto the base characteristic plane (the
z-t plane). In cases where Sw is constant, equation of 20 has the well-known analytic
solution for exponential decay
S0(t) = SO(T) exp[ -X0 (t - t)] (22)
where t is the time origin of the characteristic. The solution, equation 20, also exists only
where the saturation derivatives appearing in equation 17 exist in the solution.
The method of characteristics has been used to reduce the coupled set of first order
quasi-linear partial differential equations 18 to the set of ordinary differential equations
specified by equations 19, 20 and 21. The resulting ordinary differential equations are
easier to solve numerically. All of the remaining parts of the solution will also be written as
ordinary differential equations or, when possible, analytic solutions.
21
-------
For the kinematic model the saturation derivatives may fail to exist, because smooth
solutions of quasilinear hyperbolic equations do not necessarily exist throughout the whole
solution. A discontinuity may develop after a finite period of time (Jeffrey and Taniuti,
1964). In this case integral equations are used to find alternative solutions that are called
generalized or weak solutions. One place where this occurs is at the leading edge of a
liquid infiltrating into the soil, also called a wetting front. In the absence of capillary
pressure gradients, the gravity drainage model causes the derivative, 3S/3z, to tend to
infinity for the infiltrating fluid. The derivation of the generalized solution will be
illustrated by the water equation.
An integral form of the continuity equation 2 applied to a control volume around the
front gives
*>z *-i
J ^ at z + J az
(23)
where the integrals are evaluated from points z\ to Z2 shown on the water profile of Figure
4. The integrals can be evaluated analytically if the true front is replaced by the sharp front
shown in Figure 4. Note that the idealized profile contains the same volume of water as
does the true profile. For the idealized profile,
zf
J Sw dz + JSW dz
zf
+ Kew(SW2) - KeW(Swi) = 0 (24)
The exchange of differentiation and integration is permissible since z\ and Z2 are fixed
locations. The pressure gradient term of the flux vanishes because the integration is carried
out between points where 8Sw/9z is zero. Completing the integration gives
Kew(Swl) - Kew(Sw2) (25)
dt TKSwl-Sw2)
for the water. The similar oil equation is
Keo(S0i.Swi) -
,_,,.
(Sol - so2) (26)
22
-------
Equations 25 and 26 are the solution of equation 18 in the weak or integral sense. They are
the well-known jump conditions (Whitham, 1974, page 138). In obtaining these results,
the true fronts were approximated by sharp fronts. After the initial period, when the
capillary gradients affect the infiltration rate, their primary effect is to spread the fronts. As
a result some amount of the liquids will be deeper than this model predicts. At all times,
however, the mean displacement speeds of the true fronts are the same as those used in the
kinematic model and the model is conservative (Charbeneau,1984).
Saturation
Death
Figure 4: Profile for Generalized Solution
Taken together, the method of charactisitics solution and the jump equation 19 and
25 comprise the approximate solution for the water. Equation 19 is used where the water
saturation is continuous, while equation 25 is used to determine the location of the front.
Note that in the limit as Swi - SW2 approaches zero, equation 25 becomes equation 19.
Charbeneau (1984) continued by using equations 19 and 25, along with a two phase
version of the Burdine-Brooks-Corey model (equation 10) to develop an analytic solution
for a single water infiltration event He also provided integral solutions for multiple rainfall
events. The infiltration model of Green and Ampt (1911) with the time to ponding
determined by the model of Mein and Larson (1973) are used to supplement the kinematic
23
-------
infiltration model (see the development in Chow, et al, 1988). In addition, the evaporation
model of Charbeneau (1988) is used to modify the kinematic formulation for the water.
For the oil, the method of characteristics solution and jump equation (20, 26) and the
relative permeability model (eq. 11) describe kinematic unsteady oil infiltration in a strongly
water wet medium, where Richards' assumption holds. Due to the relative complexity of
the oil relative permeability function, no analytic solution exists for unsteady flow of oil.
THE SOLUTION FOR THE DISSOLVED CONSTITUENT
In this work the kinematic model is also applied for the dissolved constituent. This
is justified in part from the fact that advection, multiphase partitioning, volatilization, and
degradation are the major controls on transport and fate for many organic contaminants.
Dispersion serves to spread concentration fronts and will result in an earlier arrival time at
the water table. However, the kinematic model does move the concentration fronts at the
correct speed for mass conservation. In addition, it is much simpler to implement, it easily
includes other processes which are important, and it is much less data and computationally
intensive than full dynamic models. For these reasons the kinematic model appears to be
especially attractive for application during the initial phases of site investigations.
The dissolved constituent model solution proceeds as follows. For the
nondiffusive migration of a constituent initially dissolved in the oil, equation 12 and the
advective portion of equation 13 are substituted into the general mass conservation equation
1,
9 3
- (n B(So,Sw) Cw) + — (q0C0 + qwCw) = - AC TI B(S0,SW) Cw (27)
at dz
where B = B(S0,SW) = (Sw+ S0 ko) + Pb ksAl and AC is the constituent first-order decay
constant For consistency with the assumption that the flux of air is neglected; partitioning
into the air phase is not used in the mass balance equation for the constituent. A separate
model is used for volatilization, and treats volatilization as a loss from the top of the
constituent profile. This will be discussed below with the volatilization model. Expanding
the derivatives, while applying the linear partitioning relationships
24
-------
Co
cw
(28)
kw = 1
gives /
(29)
After the saturations of water and the oil are known, Cw is the only unknown in equation
29.
Applying the continuity equations for oil and water to the third term of equation 29
gives the result that
3C 3C
TlB—* + Zojki —* + Cw (-XoSotlko) = -AcTiBCw (30)
at az
which becomes, upon rearrangement
~ qw + qpko dCw _ , w „ w „
» ~" ~ /V£;V^\y "T" TJ/C C \
3t T| B(S0,SW) oz o^o.^wJ
Applying the method of characteristics results in
~df = GW(-^-C +B(S0,SW)
along
dz _ qw + q0k0
T| B(S0,SW)
The variables in equation 32 may be separated and integrated over a time interval from ti to
t2 to give
25
-------
An analytic solution can be obtained where the water saturation is constant. The oil
saturations in the third term of equation 34 vary in time according to the exponential decay
relation (equation 22). For a constant water saturation and with equation 22, equation 34
becomes
Cw(t2) = Cw(tl) exp{-Xc(t2 - ti)} (35)
This equation shows that in the absence of oil decay the constituent concentration changes
in accordance with its own first order decay. When the oil decays, B(S0(ti),Sw) >
B(S0(t2)>Sw) as there is less of the constituent dissolved in the oil, merely because there is
less oil present at time \2- All the concentrations must increase as there is less total solvent
present.
When the water saturation varies, the integral in equation 34 must be evaluated
numerically over the interval from ti to t2. This integration must be carried out along the
path that the constituent characteristic takes during the time interval. Thus the constituent
concentration is given by
I? 1
Cw(t2) = Cw(ti)exp{-Xc(t2- ti)} expj J BffiK dt f (36)
In °* w J
Evaporation of water from the profile removes solvent in an analogous manner to
biodegradation of the oil. Where evaporation is occurring locally from the profile, the
source strength in the continuity equation 2 is equal to minus the local evaporation rate, e*.
To account for the evaporation a third multiplicative correction term is added to equation
36. This term is
26
-------
exp<
t2
f
J
TiB(S0,Sw)
(37)
Qiarbeneau's (1988) evaporation model tracks the cumulative evaporation, E, from the soil
profile. The profile evaporation rate, e, is equal to -jr. In turn, the local evaporation rate is
related to e through e = J e* dz, where the integral is over the depth of the uppermost soil
water wave. For the kinematic model, the water saturation profile is known analytically,
and one can separate the saturation change at a given depth due to drainage from that
change due to evaporation. With this separation, the local evaporation rate is determined
from
3S
e* = TI — - I (evaporation) (38)
dt z
Quantitatively, the concentration change associated with evaporation is not found to be
significant, and this representation is discussed in more detail in Charbeneau, et al (1988).
For the constituent decay the analytic solution, with corrections for oil decay and
water evaporation was selected in favor of solving equation 32 numerically. Indeed the
inflation factors in equations 36 and 37 are evaluated numerically in the code (via the
trapezoidal rule). This approach was selected, since the oil's half life is expected to be long
relative to the solution step size, leading to a small correction in concentration per step. The
evaporation correction was found to be very small and difficult to implement (see example
5).
Equations 35 or 36 give the concentration at time t2, if the concentration at time ti is
known. The concentration at the beginning of simulation is determined as follows. The
constituent concentration in the released oil changes instantly upon placement in the soil
because of the local equilibrium assumed for the constituent For situations where the oil
flux is specified, the constituent concentration in water in the soil can be determined by
balancing mass fluxes across the surface of the soil:
Q>(surface) Qo = Cw(initial) (qw + qo *o) (39)
27
-------
This relation supplies the initial concentration, Cw(initial), for each constituent characteristic.
Equation 39 shows that only if the water flux is zero will the initial concentration in waste
oil be the same as that in the soil. Otherwise, Q)(soil) < Q>(surface)- When a specified
volume of oil is appb'ed to a zone of thickness, dpz, the initial concentration is calculated via
equation 12:
Q>(surface) Vo = Cw(initial) B(So,Sw) dpz (40)
where V0 is the volume of applied oil.
In the same way as integral equations apply at discontinuities in the liquid
saturations, an integral conservation law applies at discontinuities (which are assumed to
exist due to nondiffusive transport) in the constituent concentration. It is easy to show that
for a linear governing equation, the solutions are the same for the integral and the
differential forms. Equation 31 is semi-linear since the partial derivative with respect to z is
multiplied by a function of z, but the result still holds. Equation 32 is the solution for both
the integral and differential conservation equations. Another characteristic of linear and
semi-linear equations is that the speed given by 32 is independent of the concentration.
When the dissolved constituent is volatile, emissions from the soil may contribute
to atmospheric contamination. The tendency for the constituent to volatilize is related to its
Henry's law coefficient, kh- Volatilization of the constituent is modeled according to a
simplified method proposed by Thibodeaux and Wang (1982). AxMry" zone is assumed
to develop in the soil that lacks any of the volatile constituent. This zone begins at the
ground surface and becomes deeper as volatilization proceeds. Short (1986) modified the
model by adding resistance due to diffusion through a boundary layer at the surface. The
volatile flux is determined by diffusion through the soil and through the atmospheric
boundary layer.
Placing a control volume on the boundary between the wet and the dry zones, as
shown in Figure 5, gives the following integral mass conservation relationship
w
Im T3J
•dz = 0 (41)
28
-------
wet zone
Figure 5: Control Volume for Volatilization Model
where d and w refer to the dry and wet zones, m is the mass of constituent and J is the
mass flux. Applying Liebnitz rule,
w
,
I
, j
• m dz + md -g-
dZw r
- mw -- + Jw
n
= 0
(42)
Allowing the control volume to move with the speed of the discontinuity, such that the
mass inside the control volume is constant gives
dz _'
dt m<] -
- J
w
(43)
The mass in the dry zone is approximated as being zero and the flux there is the volatile
flux, -Jv. Since the top of the zone contaminated by the constituent can move downward
due to the flow of the oil and the water, Jw is equal to the constituent flux given by the
advective portion of equation 13. The speed of the advancing dry zone is given by
dz _ Jv + (q0k0 + qw)Cwd
dt B(So,Sw)Cwd
(44)
29
-------
where Cwd is the constituent concentration just ahead of the dry zone and the oil and water
saturations are those on the boundary. Equation 44 is equal to the speed of the nonvolatile
constituent (equation 32) plus the speed due to the volatilization. The additivity is due to
the linearity of the governing equation. If the speed determined by equation 44 is greater
than the speed of the leading edge of the dissolved constituent, then at some time the entire
amount of the; constituent is volatilized.
*'•'
This approach uses only the diffusive flux to model the volatilization of the
constituent The volatilized constituent is assumed to rise to the surface in response to the
concentration gradient. Implicitly, the Vapor is taken to be of neutral buoyancy. Some
vapors, however, may be denser than air and so sink into the soil. An indication of the
tendency for the vapor to sink can be determined by the following. A pure chemical's
vapor density relative to dry air, RVD, at 20 degrees C may be calculated from the
molecular weight, MW, and the vapor pressure in ton, p0, of the constituent by
29.0
RVD = _ (45)
according to Schwille (1988). (The factor of 29.0 changes to 28.5 if air at 100% relative
humidity is used instead of dry air.) The vapor pressure may be approximated from the
boiling point and the heat of vaporization using the Qausius-Claperyron equation. The
partial pressure resulting from the chemical in solution can then be determined by (Stumm
and Morgan, 1981, pg. 728)
Pa = Pog; (46)
where Csw is the water solubility of the constituent Substituting the last equation into
equation 45 gives the density of the vapor resulting from the constituent dissolved in
water. Lower vapor densities are obtained for constituents below their aqueous
solubilities, which is the anticipated condition for this modeling work.
30
-------
INITIAL AND BOUNDARY CONDITIONS
Table 1 shows the boundary conditions that are applied to the model. For oil spills
constant fluxes are specified for the oil, water and constituent fluxes. Initially constant
concentration and saturations are specified. A generalization of the initial condition for
water is discussed below in the KROPT model results discussion. For land treatment the
conditions are generally similar, except that the initial concentration and saturation are
specified over the depth of the plow zone, dpz.
NUMERICAL TECHNIQUES
The result of the assumptions used in the above discussion is that the method of
characteristics can be used for water, oil and the dissolved constituent The original
nonlinear system of coupled PDEs has been reduced to a system of nonlinear ordinary
differential equations (ODEs). The weak solutions for the discontinuities and the
volatilization model are also in the form of ordinary differential equations.
The solution of the model equations is obtained through the use of an ordinary
differential equation solver as suggested by Greenberg (1978, pg. 569). With such a
technique the solution of a system of n equations of the form
df = fn(zi, Z2, zs,..., zn) (47)
may be obtained, if appropriate boundary conditions are specified. The functions fn may
be coupled and/or nonlinear functions of the zj. For the models that will be presented a
Runge-Kutta-Fehlberg (RFK) method was selected. The details of this type of method was
originally derived by Fehlberg (1969). A sample implementation can be found in Forsythe,
Malcolm and Moler (1977).
There are several advantages to using a differential equation solver instead of finite
difference or finite element techniques for this problem. Use of the differential equation
solver leads to simplified programming in the sense that the programmer supplies the
functions on the right hand side of equation 47 to the solver in a subroutine. No additional
treatment is needed for their solution, except as detailed below. In discretization methods,
31
-------
Table 1: Initial and Boundary Conditions
Boundary Conditions: Oil Spills
Initial Conditions: Oil Spill
q « constant
q - constant
duration of oil
duration of oil
otherwise
otherwise
spill
spill
So-
C »
Sw-
S >
a
0
0
constant
Sor
ar
0
0
0
0
<
<
<
<
z
2
z
2
<
<
<
<
-
•
-
•
["constant, or
^declining
lamount
duration of rainfall
otherwise
Boundary Conditions: Land Treatment
Initial Conditions: Land TreatBent
t > 0
t > 0
C(0)
0 < z < d
pz
0 < z < d
['constant, or
q m ^declining
lamount
duration of rainfall
S - 0
o
otherwise
sv - sy(o)
S > S
a ar
0 < z <
0 < z <
32
-------
however, knowledge of specialized procedures is needed, particularly for nonlinear
equations such as these.
Second, the RKF methods test for truncation error in the solution. The solver
maintains the truncation error below a specified tolerance, by reducing the time step. This
feature operates automatically during the program execution, resulting in a variable time
step routine. The method used here, RKF1(2), uses a first order scheme to get the
solution, and an embedded second order scheme to check the truncation error. This
method was selected since it requires the fewest number of function evaluations of the
Fehlberg methods. The trade off is that smaller time steps may be required to maintain the
truncation error tolerance, since the method is of low accuracy.
Third, under certain conditions only a few of the equations need to be solved. For
example, the characteristics for nondegrading oils with constant water saturation are
straight lines, and only the equations for the discontinuities are solved. Much of the
computational work is eliminated, since the characteristic equations have an analytic
solution (i.e. a straight line).
A final advantage is that the location of the fronts is determined directly as a
function of time. In a discretization method, node points must be placed on points of
discontinuity. In these problems the discontinuities move with time, making it difficult to
specify node point locations in advance. Further, the equations for the discontinuities are
used precisely because the saturation derivatives do not exist at sharp fronts. Despite the
simplification by assuming sharp fronts, in the true fronts the derivatives tend to become
undefined at these locations. This presents difficulties for finite difference and finite
element methods.
A departure from typical usage of RKF solvers was made for these codes. In
addition to the truncation error control of the method, a fairly complicated and specialized
system of ad hoc controls is needed to assure the accuracy of the solution. For example the
following features of the solution must occur at times the solver picks for the solution:
• the beginning and end of rainfall and oil events,
• the end of maximum oil and water saturation regions,
• the origin of any characteristic,
• the time when water or oil fronts pass each other,
33
-------
• the time when "%oil banks" are created,
• irrigation and loading times for land treatment,
• and others.
The checks are implemented in a controller routine that is called by the RKF solver. A
beneficial side effect of using the controller is that the number of steps rejected for
truncation error violations is reduced, since the solver is guided to critical times in the
solution by the controller. ;
The overall structure of the program remains simple: a system of ordinary
differential equations solved by the RKF solver, subject to a set of controls. A flow chart
illustrating the program structure is presented in Figure 6. The method requires at most
three function evaluations per time step (blocks A,B and C of Figure 6). If the solution for
a time step is accepted, the third evaluation is equal to the first evaluation for the next time
step. The first evaluation is only needed for starting the solver and when steps are rejected.
In practice approximately 2.2 to 2.4 evaluations per step were needed on the average. If
the step is acceptable in regard to both the controller and the truncation error checker, then
the time step is increased and the solution proceeds. Otherwise the time step is cut and the
step is retaken.
Several additional numerical techniques are used in the RKF solver for special
purposes. In many cases the relationships in the model are nonlinear, necessitating the
solution of single nonlinear algebraic equations. An example is the inversion of the oil
relative permeability equation to determine the saturation associated with a given flux. The
bisection method is used for these problems, since it is guaranteed to converge if the search
interval contains the solution. The main disadvantage of the bisection method is that it is
the slowest method of solving nonlinear equations (e.g. Forsythe, Malcolm and Moler,
1977, pg. 157). The secant method was tried for some of the nonlinear equations but
discarded because in some cases the solutions oscillated and/or failed to converge in some
cases.
34
-------
A. Evaluate Model Differential Equations at
V
B. Evaluate Model Differential Equations at
to* O.SAtj Zj - zo « 0.5f(zo)
C. Evaluate Model Differential Equations at
<„ * «'• '2 - V 256 f(V * IH f(2l>
D. Check for Truncation Error
TE- |5IT(f<2o> -«2
E. Check For Critical Conditions
in Controller
F. If
TE Within Tolerance
and
No Controller Checks
Violated
C. Yes Advance Time Step
At
H. No Reject Tine Step
ot «• 0.5* At
Figure 6: Flow Chart of Runge-Kutta-Fehlberg Algorithm
35
-------
Two variants of the bisection method are used for functions known at discrete
points. The first of these is used to find between which two entries in an array another
value lies for the purpose of interpolation. This is used extensively in computing the
saturations at various depths in profiles, from values on characteristics. An efficient binary
search technique was implemented for this purpose (the method was based on a routine
appearing in Forsythe, Malcolm and Moler, 1977, pg. 79). The second is used in cases
where many bisection solutions would be made of the same function. To save
computational effort for these functions, a set of 50 function values are saved at the
beginning of the program. When the inverse of the function is needed, a binary search is
done on the saved values, and inverse quadratic interpolation (Abramowitz and Stegun,
1965, pg. 882) is used to determine the interpolated value.
Functions are integrated with a 7 point Gauss Quadrature routine (Abramowitz and
Stegun, 1965, pp. 887,916). This was chosen to provide reasonably high accuracy, while
utilizing relatively few function evaluations. Quadrature is used to obtain the mass of
substances in the soil, by integrating the saturation (and concentration) profile. Note that
mass balances reported in the discussion of results are obtained numerically, and contain
some error.
Both programs produce the following three types of output. First, after echo
printing the input data, a record of the locations of the fronts and their saturations,
concentrations and mass balances at each point in time of the solution is produced in an
output file. Second, in a separate file, concentration and saturation vs. depth profiles for
selected times are written. Third, at selected depths, concentration, saturation, and flux
histories are written in another file. These provide the user with a record of the flux of
pollutants passing a certain depth at each time. Specific details of these and other features
of the implementations are presented below.
36
-------
SECTION 6
MODEL IMPLEMENTATION AND RESULTS
Two implementations of the kinematic model will be presented in this section.
First, a simplified implementation is presented that assumes a uniform constant water
saturation is present throughout the profile. This model is called KOPT (for Kinematic
Oily Pollutant Transport), and is intended to be used as a screening tool for oil spills. It
addresses the questions of how far spilled oil might go into the soil and how soon it might
get there. In this case, no rainfalls are assumed to occur, the specified water saturation is
taken as representative of climatic conditions. This model is usable with other models of
subsurface oil transport. Specifically, it can be used to supply input data to the model of oil
spread along the water table developed by Roberts-Schultz (1988). This discussion will
also be used to introduce the elements of the method of characteristics solution, in
preparation for the more complicated variations that occur in the full model.
The full implementation, presented later in this section, includes the effects of
rainfall on the spilled oil, and is called KROPT (for Kinematic Rainfall and Oily Pollutant
Transport model). It includes the modifications of the kinematic model for water
(infiltration and evaporation) discussed above, along with a model for the stochastic
generation of rainfalL The character of the results will show that the incoming water has an
extreme influence on the flow of the oil. In some cases, the water flux is too high for the
kinematic model, and the method breaks down. A discussion of how to generalize the
model for these cases is presented in Weaver (1988).
KINEMATIC OILY POLLUTANT TRANSPORT MODEL (KOPT)
In the KOPT model, rainfall is represented by using an average infiltration rate.
The corresponding water saturation is determined through inversion of the effective
hydraulic conductivity relationship:
qw = Kew(Sw) -> Sw(avg) = Kew-1(qw(avg)) (48)
37
-------
which for the water relative permeability (equation 10) is
Sw(avg) = Sw + (1 - Sw) j\ £ - (49)
Sw(avg) is used as the water saturation for the entire profile. When qw is a constant, the
phase conservation equation for water shows that dSw/dt = 0, and the equation is
effectively solved.
The oily pollutant is assumed to enter the soil and fill some portion of the pore
space that is not occupied by water or trapped air. If the oil leaks at a rate which is less
than the maximum rate that the soil can accept, then equation 26 is used to determine the
(kinematic) speed of the invading oil. The oil's effective hydraulic conductivity, which is
equal to its flux rate in the kinematic model, is inverted to obtain the oil saturation during
the application, So(max)- In doing this, the average water saturation, Sw(avg)» is used in all
the oil equations.
When the application rate of the oil exceeds the kinematic capacity of the soil, the
Green-Ampt (Green and Ampt, 1911) model can be used for the oil flux. For cases with
applied head of Hs,
qo(So(max)) = Keo(So(max)) l + (50)
where zf is the depth to the oil front Initially the pressure gradient Hs/zf is infinite, because
zf is zero, but is reduced as the oil front invades the soil. The head Hs may be due to
surface ponding and/or capillary pressure. This situation adds a constant head boundary
condition to those shown in Table 1. Even though there is a head applied at the surface,
the amount of water in the soil is unchanged in this situation, as the water saturation is
considered representative of climatic conditions.
In cases where the supply of the oil is finite, the kinematic model ties in with the
end of the oil event to model the drainage. Drainage modeling is not possible in a NVpure"
Green-Ampt model and is a special feature of the present model.
38
-------
The solution is obtained in the following manner. The flux from equation 50 or the
kinematic flux is used in the jump equation, 26, to give the oil front speed,
dz
I! So(max)
The method of characteristics solution (equation 20 and 21) is used to give the characteristic
directions: either through inversion of the derivative of the relative permeability curve when
AO is zero, or by solution of the differential equations when AO is not zero. The distribution
of the constituent is obtained from the method of characteristics solution of the constituent
equations (32 and 33), and the volatilization model (equation 44). The constituent moves at
rates that are determined by the partitioning phenomena, and the speeds of the oil and the
water. There can be an apparent cleansing of the oil. This can be in front of the
constituent's maximum depth, if the constituent moves slower than the oil. If the oil moves
slower than the constituent, then the constituent can be leached entirely out of the oil. In
reality dispersive forces and non-linear partitioning smear the assumed sharp constituent
fronts, and so make the cleansing less complete.
A brief example of oil migration will be presented to introduce the types of results
obtainable from the model. More examples are presented in following sections. Figure 7
shows the base characteristic plane and the projection of a few displacement paths in z-t
space. Shown here is the result of a sixty day simulation from a ten day oil release at 2
cm3/cm2/day into a dry sand. Other details of this case (example 1) are found in Table 2.
During the oil release, the oil fills 37.6% of the pore space and moves at its maximum
speed. This is indicated by the straight line labeled AC on Figure 7.
Once the supply of oil ends, drainage begins and there is in the profile a region of
constant oil saturation (triangle ABC), associated with the oil release, that will now be
replaced by a region of variable saturation, associated with the method of characteristics
solution for oil drainage. The displacement path (called a characteristic) of the residual
saturation, Sor, remains at the surface, because its speed is zero by definition. Bounding
the region of maximum saturation, ABC, is the characteristic corresponding to S0 = So(max)
= 0.3458. In between are characteristics for all intermediate values of S0. They are
represented in the figure by characteristics for each 0.01 increment of oil saturation. The
greatest oil saturations are found the deepest in the profile, because the derivative of the kro
39
-------
Depth
(cm)
20 30 40
Time (days)
— oil Front
— Oil Characteristics
Figure 7: Oil Spill in a Fine Sand
40
-------
Table 2: Input Data for Examples 1, 2 and 3
Soil (a)
K
xvs
n
pb
*!L
cm/day
3
gr/cm
2
Example 1
Fine Sand
209.7
3.7
0.314
1.72
Example 2
Silt Loam
39.4
0.207
0.462
1.70
1
oc
Water
Svr
qVC
u
Oil
u
P
S
half life
lag time
loading
time
volume
depth
cm/day
gr/cm/day
gr/cm/day
gr/cm
days
days
cm/day
dags 2
cm /cm
cm
0.1669
0
872
1285
0.73
0.10
2.0
10
0.039
0.015
872
471
0.65
0.05
360
45
1.0
25
Constituent
C Unit)
half life
lag time
mg/1
days
days
oc
Climate
Temperature
Evaporation
Relative
Humidity
cm/day
200
134 (b)
0.15 (b)
60 (b)
20
0.1
0.5
Example 3
Sandy Clay Loam
118.5
0.3446
0.328
1.70
1
0.2287
0.104
872
1000
0.65
0.10
180
30
1.7
30
200
360
45
50
0.10
10
15
0.15
0.5
a) Example 1 soil data from Brooks and Corey, 1964; Example 2
soil data from Brakensiek et al., 1981
b) from Ryan et al., 1987
& 3
41
-------
function is such that the highest oil saturations have the highest speeds. Thus the pattern of
characteristics which determines the shape of the oil profile results from the method of
characteristics solution of the oil equation and the relative permeability function. This
pattern is called a centered simple wave.
Figure 8 shows oil profiles for 10, 1 1, 30 and 60 days. In these and all subsequent
profiles, water and the total, liquid saturations are plotted, and the oil saturation is read as
t
.the difference between the two. At 10 days the profile has not yet begun to drain and the
oil saturation is equal to 0.3458 at every point above the oil front. By 1 1 days drainage is
occurring from the surface down to a depth of 80 cm. The region from 80 to 202 cm still
has the original saturation of 0.3458. Once the oil saturations that are less than So(max)
reach the oil front, then its speed begins to be reduced in accordance with equation 26. As
time goes on the speed of the front slows as it intersects slower and slower characteristics.
Figure 8 also shows the oil profile for 30 and 60 days where the entire profile is
undergoing drainage, and the oil front saturations have dropped to 0.2201 and 0.1809,
respectively. In the kinematic model it is only after an infinite amount of time that the oil
saturation is at residual throughout the entire profile at which time, the motion of the front
stops. Practically, however, the speed of the front becomes negligible at much shorter
times. When this occurs the depth of the front, and thus the depth of the entire oil
contaminated zone, ZQCZ. is given by the simple calculation,
zocz = (52)
where q0 is the application rate of oil, and Ato is the duration of the oil spill. Equation 52
gives the ultimate depth of a nondegrading oil. It is a useful quantity, if the residual
saturation is known.
As noted previously, the relative permeability functions are subject to hysteresis.
Solution of the kinematic model shows that, where the classical method of characteristics
solution holds (equations 20 and 21) drainage is occurring in the oil profile (see Figure 8).
The drainage relative permeability model (eq. 11) applies to this situation. Imbibition (or
oil wetting) occurs at the oil front. The true shape of the oil profile at the oil front is
determined by imbibition relative permeability and capillary pressure functions. An
imbibition relative permeability curve is necessary for modeling this shape. In this model,
however, the shape of imbibing profile is not modeled, as the front is assumed sharp. The
42
-------
u
-50
-100
-150-
-200-
-250-
Depth ?nn
in cm "30°"
-350-
. -400-
-450-
-500-
-550-
water
•
^ V Xfc 1 1 1
•\ v
' i
1 i
1 i
"i "1
1 1
1
\
\
total
liquid
profiles
air
-4 1 1 1
0.0 0.2
0.4 0.6
Saturation
------------- 10 days
---- 11 days
- 30 days
------- 60 days
0.8 1.0
Figure 8: Oil Profiles
43
-------
consideration that must be addressed is the speed of the front, which is determined by the
hydraulic conductivities on either side of the front Since the hysteresis loops connect the
drainage and imbibition capillary pressures at the point where the draining profile becomes
an imbibing profile the relative permeability is the same. The imbibing kr is not necessarily
the same as the drainage kr, since a scanning loop connects the drainage and imbibition
curves. The drainage kr is in error over the range of saturations where the capillary
pressure follows the scanning loop. The error is expected to be within the same range as
those caused by the other assumptions of the model, since the magnitude of the difference
between imbibition and drainage kf is often < 10%.
KOPT Model Results
Example 1: OH Spill into a Fine Sand-
Figure 9 shows the maximal extent of the oil for the same oil loading as in Figure 7,
but extended for 180 days. For clarity none of the characteristics are shown, but the oil
saturations vary as discussed for Figure 7. The maximum numerically computed volume
balance error for the oil is 0.15% and occurs at 180 days. Also shown in Figure 9 is a
simulation where the oil has a half life of 180 days in the soil. Note that for both cases
most of the total depth achieved is done so in the first 20 days of the simulation. However,
during the first few days there is virtually no difference between the two cases. At those
short times (relative to the half life), the decay does not significantly reduce the effective
hydraulic conductivity of the oil. Later as the decay becomes significant, the front for the
decaying oil slows. The KOPT model also allows a lag time to be incorporated to account
for acclimation of soil micro-organisms (see Table 2). In contrast, the non-decaying oil
continues to move downward throughout this time, reaching a depth of 504.7 cm at 180
days, 101.9 cm deeper than the degrading oil Equation 52 predicts the final depth for this
case to be 796 cm, which will not be achieved by the degrading oil.
Profiles at 60 and 180 days are drawn in Figure 10. At 60 days and for the
decaying oil at depths above 3.5 cm, all the oil has been reduced to saturations below
residual and is thus immobilized. The oil's decay is creating a region of immobilized oil in
the profile. This begins at the ground surface and works downward with time. The dotted
line on Figure 9 shows the locus of points where decay has reduced the oil saturations
below residual. Thus all points above and to the right of this curve are at oil saturations at
44
-------
-600
50 100
Time (days)
150
Oil Front (no decay)
Oil Front (180 day half life)
Locus of Residual Oil Saturations
Figure 9: Example 1 • Base Characteristic Plane for Oil Spill in
Fine Sand, 180 Duration
45
-------
V-
-50-
-100-
-150-
-200-
-250-
Depth _,OQ.
in cm 30°
-350-
-400-
-450-
-500-
-550-
-600-
0.
. , M
u
ti
;•
u
t\
«• *
j
t
-
i
t
i
t
>
t
i i i i
0 0.2 0.4 0.6 0.8 1.
0
Saturation
60 days conservative
180 days conservative
60 days degrading
— - 180 days degrading
Figure 10: Example 1 - Conservative and Decaying Profiles at
60 days and 180 days
46
-------
or below residual. At approximately 160 days the oil saturation at the front is reduced to
residual and all motion of the oil stops.
The sensitivity of the solution for the nondecaying oil to the magnitude of the
residual oil saturation is shown by Figure 11. All parameters of example 1 were
unchanged, except for the residual oil saturation. The greatest depths were attained for the
smallest residuals. This result is due directly to the increased oil relative permeability at a
given oil saturation, with decreasing residual oil saturation (equation 11).
Example 2: Gasoline Leak into Silt Loam—
The next example is of a gasoline leak in a silt loam soil (Table 2). The water flux
is 0.055 cm/day (20 cm/year), which results in a water saturation in this soil of 0.6108. A
constituent initially present in the gasoline at 200 mg/1 is partitionable into the water and
slightly sorbed. The partition coefficients are for benzene (Ryan et al., 1987). When
placed in the soil, the concentration in the gasoline instantaneously changes to 199.70 mg/1,
according to the flux relationship (equation 39). The water concentrations are 134 times
lower than these values, according to the partitioning relationship. If only water samples
were analyzed, much of the mass of the constituent would be missed.
For a volatile constituent (kh = 0.15), the upper boundary is indicated as "volatile
front" in Figure 12. The mass balance for the constituent is shown in Figure 13, where
the sum of the mass remaining in the soil and the mass volatilized is seen to remain very
close to the amount placed into the soil (0.5 mg/cm). The maximum mass balance error for
the constituent is 0.3%. Volatile losses are calculated using the volatilization model
described above. Profiles for the volatile constituent are shown in Figure 14. These show
that the contaminated zone moves downward with time, and that concentrations increase
throughout the simulation. The concentrations are calculated by equation 35, which models
the effect of oil decay on the constituent concentration. The mass balance indicated in
Figure 13 includes the increase in concentration. The limiting constituent concentration in
water is 24.28 mg/1. None of the concentration profiles contain this value, since the oil is
present throughout the simulation. The decaying oil profiles are shown in Figure 15.
Gasoline is highly volatile, because of the large amounts of volatile constituents it
contains, so some loss of the liquid is expected. This model (and also the KROPT model)
47
-------
Oil Spill
Fine Sand
Depth -400
(cm)
50 100
Time (days)
150
0.01 Sor
0.05 Sor
• 0.10 Sor
0.15 Sor
0.20 Sor
Figure 11: Sensitivity of Results to Residual Oil Saturation
48
-------
-80
-160--
-180--
-280
volatile
front
-4-
100 200 300 400
Time in Days
Oil Front
500
. constituent Limits
Figure 12: Example 2 - Base Characteristic Plane
49
-------
0.55
0.50- - r-~- •"•»"«•
Mass
in mg/cm2
0.45-
0.40-
0.35-
0.30--
0.25-
0.20-j
i
0.15-'
200 400
Time in Days
- • Mass Remaining in Soil
— Mass Volatilized
— Total Mass
Figure 13: Example 2 - Constituent Mass Balance
50
-------
v/-
-20-
-40-
-60-
-80-
-100-
-120-
-160-
-180-
-200-
•220-
•240-
260-
280-
0.
I i
0 0.5 1.0 1.
.ISS^fe
i i
1
i
i i
T
*
I
*
J
I
5 2.0 2.5 3.0 3.5 4.
0
Aqueous
Concentration mg/1
45 Days
180 Days
365 Days
1 500 Days
545 Days
Figure J4: Example 2 - Concentration Profiles
51
-------
in cm
v-
-50-
-100-
-150-
-200-
-250-
III
water
1 1 1
!\
I
I
I
i
I
ii
ii
]
il
I
i
i
ii
:•
B
i*
it
ii
L — ;
air
1
0.0 0.2 0.4 0.6
Saturation
0.8
1.0
45 Days
180 Days
• 365 Days
: 500 Days
545 Days
Figure 15: Example 2 - OH Profiles
52
-------
are not suited for modeling the loss of the immiscible liquid. Only the volatile loss of a
dissolved constituent is considered. Some loss of the gasoline is expected, but the model
results are conservative for these applications insofar as no loss of gasoline gives the
maximum amount of gasoline in the soil.
Example 3: Land Treatment in a Sandy Clay Loam-
A hypothetical land treatment example is presented in Figures 16-20. There are
330 barrels of oil per acre (1.70 cirP/cm2/day), corresponding to the maximum typical
industrial loading reported by Ryan et at (1986), uniformly placed in the top 30 cm of the
soil. Since its maximum saturation is above the residual (0.1727 vs. 0.1000), the oil
moves 12.7 cm out of the bottom of the plow zone before it is immobilized. If the oil were
nondegradable, equation 52 indicates that it would have moved 22 cm out of the plow
zone. With a half life of 180 days, 30 days required for acclimation, and constant residual
saturation of 0.10 assumed, the oil is immobilized after 104 days.
The vertical extent of volatile constituent is indicated in Figure 16. The infiltrating
water leaches the constituent from the oil. After 455 days all of the constituent is deeper
than the now immobilized oil. Once this happens all constituent characteristics are straight
lines, since the motion is only due to the water flux, and there is no oil present (see
equation 33). Where the constituent is above the oil front, the characteristics' curvature (as
evidenced by the dotted line) is caused by partitioning into the degrading oil. In this case
and as would be typical for land treatment, the oil only moves a small distance. Its motion
has a minimal effect on the ultimate migration of the constituent
Concentration profiles for times from 30 to 730 days are presented in Figures 17
and 18 for the nonvolatile case. At 30 days the concentration is equal to the initial value
throughout the profile, since neither the oil nor the constituent have started to decay. At 60
days biodegradation of the constituent has reduced the concentration at the leading edge,
while concentrations above the oil front are increased because of oil biodegradation. The
straight line connecting these values is an artifact of the number of characteristics used for
the constituent At 30,60 and 180 days only one constituent characteristic is below the oil
front and linear interpolation is used to find intermediate concentrations. By 270 days
enough characteristics are below the oil front to give the profile curvature. The general
53
-------
Depth
in cm
v-
-20-
.-40-
-60-
-80-
-100-
—120-
-140-
-160-
-180-
-200-
-220-
-240-
_?£r>-
* — -.. I
V
••;* _ .._
"x.^'x """••
^ X:-^
N \ \
N \ \
X \ X
\ \
X \
\ \
V \
N \
\ \
^ \
\'
1
I I
"""\
*.
*,
t
%
\
\
\
\
x \
V. \
^ \
\ \
N \
\ \
\ \
x \
\ x-
X \
\ \
\»
»
\ '>
\
1 1
200 400 600
Time in Days
Oil Front
Constituent—Leading Edge
Volatile Front k^ *• 0.00
Volatile Front kh - 0.10
Volatile Front kh - 0.20
Figure 16: Example 3 • Base Characteristic Plane
54
-------
-50-
-100--
Depth
(cm)
-150-
-200--
-250- •
Concentration in Water in mg/1
30 Days
- 60 Days
. go Days
180 Days
270 Days
Figure 17: Example 3 - Concentration Profiles
55
-------
-50--
-100--
Depth
in cm
-150--
-200--
-250--
0 2 4
Concentration in Water in mg/1
360 Days
450 Days
• 540 Days
-• 630 Days
730 Days
Figure 18: Example 3 - Concentration Profiles (continued)
56
-------
0.035
Mass
mg/cm2
0.010-- /
0.005-/.
0.000
200 400 600
Time in Days
- Mass in Soil kh D 0.00
............. Mass in Soil kh = 0.10
---- Mass in Soil kh « 0.20
------- Mass Volatilized kh - 0.10
---------- Mass Volatilized kh - 0.20
Figure 19: Example 3 - Volatilization Comparison
57
-------
trend of the 60 and 90 day profiles continues until the constituent is entirely ahead of the oil
front (at 455 days). All concentrations then decrease due to constituent biodegradation.
The same case was run with kh = 0.1 and kh = 0.2 to illustrate the effects of
volatilization. The solid line in Figure 19 shows the nonvolatile case. The mass begins to
be reduced after the lag of 45 days has ended. At 405 days, one half of the constituent's
mass is gone. This corresponds to the half life of 360 days. Volatilization reduces the
mass in the soil further. Note that for both volatile cases shown, there is more total mass
(mass volatilized + mass in the soil) than for the nonvolatile case. This is because
volatilized constituent has left the soil and cannot be degraded. It is only the remaining
fraction of the mass that is subject to degradation.
Concentration histories are shown in Figure 20. At the bottom of the plow zone
(30 cm), the concentration is at the initial value (equation 40) until the oil degradation lag
ends. The kink in the concentration history curve occurs when the constituent
biodegradation begins (after 45 days). Thereafter the concentration increase caused by oil
biodegradation is reduced by constituent biodegradation. In this case the half life of the oil
is the shortest, so that the concentrations always increase in the presence of decaying oil.
This occurs until the water leaches the constituent out of the plow zone at 501 days. After
400 days the constituent reached 150 cm in depth (approximately the depth to the lower
treatment zone). Lower concentrations are encountered, since oil biodegradation no longer
inflates the concentrations.
KOPT Model Implementation
The specific details of the implementation of the KOPT code are presented in
Weaver (1988). The programming synopsis is shown in Table 3. For each step in time,
the right hand sides of the oil equations are evaluated first As the RKF routine performs
each intermediate step, the oil solution must be known before the constituent equations are
solved. There is storage space reserved for the oil front and oil characteristics. The oil
characteristics are only solved when the oil is biodegradable. Otherwise the solution of the
oil characteristics are straight lines which need not be evaluated numerically. A set of 51
characteristics is used for the oil. Ten of the set are used during the oil's application and
the remaining 41 are used for the drainage pattern.
58
-------
5.0
Cone
mg/1
3.5-
3.0--
2.5--
2.0--
1.5--
1.0--
0.5-
0.0
200 400
Time in Days
600
Concentration in 'water 30 cm Depth
Concentration in Water 150 cm Depth
Figure 20: Example 3 - Concentration Histories at 30 and 150 cm Depths
59
-------
Table 3: KOPT Programming Synopsis
Fronts Ch«r«ct«ritst>cs Notes
Water 0 0
Oil
Nondecaying 1 0
Decaying 1 SI
Constituent
nonvolatile 0 11
volatile • 1 11
\
\
Constant value used
10 characteristics used during oil release
41 characteristics used for sinple wave
Front is for the volatilization Model
Itea
Controller Checks
Notes
Water
Ron*
Oil
Constituent
Simulation
Beginning of Event
Cnd of Event
End of S
o(aax)
At le«st S steps through release
Starting tiae of each characteristic
See figure 4.1, point D
Rondecaying oil
Decaying oil
Tiae when characteristic passes oil front
Volatilisation of entire >ass of const. Remove all constituent equations
Place solution on even tiae step
Siaulation Ending tiae
Profile tiaes
60
-------
For cases with nondegrading oils, no oil characteristics need to be computed by the
RKF solver, since they are straight lines. This makes for a very efficient solution as only
one equation is needed to determine the oil solution (that for the front). If the oil is
biodegradable, all oil saturations decrease with time. Thus the speed slows down and the
characteristics bend, in accordance with equation 21. Characteristic speeds are still given
by equation 21, but are time dependent according to equation 22 and the kro function.
Points defining the curved characteristics must be stored temporarily, since equation
26must be used to determine the speed of the front Note that many characteristics are
used to determine the solution initially, but as time goes on fewer and fewer are used in the
solution. At later times, however, the oil profile is flatter and thus the loss in accuracy due
to fewer characteristics is offset by the shape of the profile to some degree.
This problem does not exist for a nonvolatile constituent, since the linearity of the
solution assures that characteristics never intersect the leading edge of the discontinuity.
Thus the solution is always determined by the same number of characteristics. The
equations for these are solved last, since the oil and water saturations and fluxes must be
known. Eleven characteristics are used for the constituent. Since the constituent equation
is linear under the model assumptions, no equation is needed for a front; the deepest
characteristic corresponds to an invading discontinuity in concentration. One equation is
used for the volatilization model
KINEMATIC RAINFALL AND OILY POLLUTANT TRANSPORT -
(KROPT)
In this section the implementation of the kinematic model for flow of both the water
and oil is presented Instead of using constant water saturations and fluxes, the unsteady
migration of water is also modeled. This makes the model more realistic by considering the
effect of a series of rainfall events on the oil and the dissolved constituent
The Formation of Oil Banks
A difficulty of the method of characteristics solution approach for this problem is
that the qualitative nature of the solution must be determined in advance by the programmer
(i.e., not by the user). The classical method of characteristics solution does not indicate
61
-------
where discontinuities appear. If the proper discontinuities are not included, then the
solution fails to conserve mass. A crucial example occurs as water begins to infiltrate into
the soil. Ahead of the water front, the water saturation is at some pre-existing value,
while behind the water front it is determined by the incoming water. Oil moving at constant
flux across the water discontinuity would also experience a change in saturation, due to the
dependence of oil relative permeability on water saturation (Figure 3). At the location of
the water front, there is a discontinuity in the oil saturation. This means that both saturation
derivatives in the oil (equation 18) fail to exist at the water front, implying that the jump
condition (equation 26) must be applied to the oil.
The speed of the jump is fixed by the water front speed in the kinematic model, and
the upstream oil saturation is fixed by the boundary condition. The oil saturation ahead of
the water front, then, will depend on these two. At some point further just ahead of the
water front, the oil saturation is again a value determined by the boundary condition of the
problem, since it is as yet unaffected by the invading water. Thus there is a second
transition (which is also a jump) down to the saturation at the deeper depth. This
phenomenon also occurs during oil production activities (e.g. Helfferich, 1981), where the
region of increased oil saturation is called an oil bank. Oil banks will be discussed at length
in the following section. Control volume analyses of oil banks for steady and unsteady
flows are presented in Weaver (1988).
Quasi Steady Oil Banks--
The first step in studying the oil bank phenomena is to consider a simplified case
where there is no drainage. In such a situation the method of characteristics solution plays
no role, and the jump equations can be studied by themselves. As an example, oil is
applied at a constant rate for ten hours to a sand where water is at residual saturation
throughout the profile. Five hours later, additional water begins to infiltrate, also at a
constant rate. Figures 21 and 22 show the base characteristic plane and a profile for this
case. Table 4 contains the input data used.
Ahead of the water front (line BH on Figure 21), the relative permeability of the oil
is subject to a water saturation of SWT (0.1669). The oil flux (0.2 cm/hr) gives, through
inversion of the hydraulic conductivity for oil, the oil saturation, S03, as 0.3730. Down to
the depth of the water front, the water saturation is Sw(max) = 0.4536, determined by
62
-------
-10--
Depth -15--
-20--
-25--
10 12 14
Time
Water-
Leading Edge of Oil Bank
Oil Front
Constituent
Figure 21: Quasi-Steady Oil and Water Migration
63
-------
Profile at 6 Hours
u-
-2-
-4
-6-
zfs
Depth -8-
-10-
-12-
-14-
1 C_
— AOT
0.
wa
•
0 (
— i 1
ter 5
Sc2
oil
So3
t i
).2 0.4
D
C
— 1
1
ai
i
i
).6 C
i
r
1
1.8 1.
0
Saturation
Oil
Water
Figure 22: Quasi-Steady Oil Profile
64
-------
Table 4: Input Data for Steady State Oil Bank Discussion
SoU
on
Pb
fbc
Water
WT
loading
time
P
SOr
loading
time
Constituent
Co(init)
ko
koc
cm/hr
gm/cm/hr
cm/hr
hr
gm/cm/hr
gm/cm'
cm/hr
hr
mg/1
Fine Sand
8.74
3.7
0.314
1.50
0.01
0.1669
36.36
0.2
1
55.44
0.73
0.02
0.2
1
100
200
10
inversion of hydraulic conductivity for water. The boundary between the two regions is the
discontinuity, which moves at the speed of the water front:
_ Kcw(Swfmax))
(53)
Applying the jump equation to the flow of the oil, with the speed fixed, prescribes a
relationship between the oil saturations on either side of the water front as follows:
65
-------
dzwf _ dzrfr _ Keo(So2,S wiO-KeoCSpi ,Sw(max))
Tj(S02-Sol)
dt ~ dt ~ -- - » (54)
The oil saturation behind the water front, S0i, can be determined, since the oil flux at the
surface remains the same. Inverting the relative permeability function gives this value
(0.3149). Since the effect of the water is modeled as causing an increase in the
permeability of the medium to oil, then this new oil saturation will be less than the original
oil saturation, S03 (0.3730). This indicates that a lower oil saturation is required to give a
certain flux, when the water saturation is increased.
A second jump is created, that is, a jump from S02 down to the original oil
saturation, S03- The speed of this jump is given by
_ Keo(So2.Swr)-Keo(So3,Swr)
dt TKS02-S03)
The space between the two jumps (region BEIHF on Figure 21) is an oil bank. The speed
of the second jump must be faster than the speed of the first for this solution to make sense.
Figure 23 shows the speed of the water front (rear of the oil bank) compared with the front
the oil bank. Only one kr curve is illustrated here, as the water saturation ahead of the
water front is constant. The lines drawn on the curve are directly proportional to the speed
of the fronts (see equation 26). As long as S02 is also constant, the increase in oil
saturations from S0i to S03 assures that the speed of the front of the bank is greater than
speed of the rear, 7.03 cm/hr and 2.22 cm/hr, respectively.
Examination of the saturation profile in Figure 22 reveals that the volume of oil that
would have been in the region where water is present (zwf (S0i-S03) TJ) is exactly balanced
by the extra oil in the bank ((zfl,- zwf) (S02 -S0l) "H). In effect the infiltration of the water is
causing this amount of oil to be displaced in front of the water. Eventually the front of the
oil bank overtakes the original oil front and S03 disappears from the profile. Later the water
front overtakes the front shock and S02 in turn disappears.
The solution was obtained by using only the jump equations (for the water front,
the leading and trailing edges of the oil bank, and the oil front), without resorting to volume
66
-------
Keo(SO.Sw) 1.0--
0.0
0.0
Sw - 0.17 m Swr
Rear of Bank
Front of Bank
Figure 23: Effective Oil Conductivity and Oil Bank Speeds
67
-------
balance computations. A volume balance of the oil is then an independent check of the
solution. Obviously in a correct solution of a conservative problem, the volume of oil
applied must at all times equal the volume of oil in the soil. For all the cases run, only
insignificant errors in volume were found. Since this is a necessary condition for the
solution to be correct, the correct volume balances prove existence of the solution. Further,
it is the unique solution, as the two jump equations which are in effect being solved
simultaneously have two unknowns: the speed of the leading edge of the bank, and the oil
saturation in the bank S02- The other oil saturations (S0i and S0s) are determined by the
boundary condition of the problem.
The necessity of assuming the existence of the oil bank can be shown as follows.
If the discontinuity in oil saturation were neglected, the boundary condition would still
specify that S03 is greater than S0i- As the water front invades the region where the oil
exists, the oil saturation would drop from S03 to S0i- The total mass, however, is the sum
of the oil mass above and below the water front. Above the water front there is less mass
than before the water infiltration started, while below the water front there is the same
amount of mass as would exist if no water was invading. Thus there would be a net loss
of mass caused by the infiltrating water.
A dissolved constituent carried by the fluids is indicated in Figure 21. Here the
constituent was hydrophobic and slightly sorbed (ko = 200, ks = 10, pt = 1.5). The
deeper line labeled ^constituent" (AG) on the figure shows the leading edge of the
constituent Its speed is given by equation 33. The speed is independent of concentration
and depends only on the saturations, fluxes, and partitioning relationships. The oil bank
causes the speed of the constituent to increase at point C, because the oil flux and saturation
are elevated in the bank. When the bank passes (point F), the speed slows in accordance
with the oil and water fluxes and saturations behind the bank.
Initially, equilibrium concentrations are established in the soil according to equation
39. When the boundary condition changes, as at time t = 5 hr, so does the initial
concentration in the soil. A discontinuity in concentration is created that then begins to
migrate downward at the concentration velocity (line BD). Above this line the
concentration in water is 0.4975, while below it is 0.5000. The small difference in
concentration results from the constituent being hydrophobic; most of the constituent
resides in the oil. Note that this solution results in near perfect mass balances for the
constituent.
68
-------
Unsteady Oil Banks—
Of course, problems of realistic interest include drying of the profile and the
analogous draining of the oil. Once again the increased water saturations due to a rainfall
event cause the oil permeabilities to increase. Behind the water front the oil would move
faster than if no infiltrating water was present The abrupt change in oil flux occurring at
the location of the water front again must lead to the existence of an oil bank, if mass is to
be conserved.
Figure 24 shows a case where oil was applied to the soil for an hour. Later water
was applied, also for an hour. Now it is the drainage or unsteady portion of the oil event
that is affected by the water. There is again a discontinuity in the oil's domain, driven
solely by the water flux. The jump condition can be applied to both edges of the oil bank.
In general the conditions for the rear of the bank are
_ dzrb _ Keo(S02,Sw2)-Keo(Soi>Swi)
dt dt Tl(So2-Sol)
and the leading edge
dzfb _ ,w,w ™.
dt
In this example the water saturation ahead of the water front is at residual (SW2 = SW3 =
SWr)» where the numeric subscripts on the saturations refer to Figure 25. Since the oil
saturations may vary throughout the profile, the various saturations are those immediately
adjacent to the jumps (Figure 25). At the water front the oil saturation jumps from S0i to
So2» while at the leading edge of the bank it jumps from S03 down to So4.
A complication arises in that there are three unknowns: the oil saturations S0i, S02»
and the speed dzfb/dt, while there are two jump equations (S03 depends on S02 and So4 is
known from the previous oil drainage pattern). This presents a difficult problem, which
has been resolved via numerical experimentation. This work has shown that the only
solution that conserves mass is the solution where the water displaces all the mobile oil
69
-------
-10--
Depth _.r _
in cm ~
V-'•.-.»• -. \ • XV.-\\-» \ N- \ >.
%\\\\ \ W\\\\\ x \ V-^.
^\\XH\V\\V
-20--
-25--
Time in Hours
Oil Front
• Water Front
Oil Bank
Oil Characteristics
Figure 24: Unsteady OH and Water Migration
70
-------
Saturation
Depth
fowf _ dZffr _ Keo(S02,Sw2)
dt *
Figure 25: Hypothetical Unsteady Oil Profiles
from behind the front The oil is left at its residual saturation after the passage of the water
(S0i on Hgure 25 is always equal to SOT). Equation 26 becomes
(58)
Here again, existence of a solution requires conservation of the oil volume. In effect, an
integral continuity equation, which is applied to the entire profile, is being added to the
system of equations that are solved.
The base characteristic plane shown in Figure 24 and the accompanying profiles
represent a solution that conserves volume throughout the entire event. The maximum oil
volume balance error (computed numerically) was 0.083% at 6 of the 8 solution points
between 6.25 and 6.40 hours. Using a code where the water solution was purely analytic,
maximum errors of 0.02% were achieved. Thus existence of the solution is proved
71
-------
numerically. Several other proposed solutions which allowed oil drainage behind the
water front did not conserve volume, and therefore are not valid solutions. They fail to
conserve mass as soon as the water event begins. This is when dzrt/dt is a constant and
thus no unsteady water flow is involved. This means that the attempted solutions fail for
the simplest case: unsteady flow in the oil only. Since the system has been shown to
remain hyperbolic when the water saturation is constant, that these trial solutions fail is not
due to the system of equations becoming nonhyperbolic.
The oil bank is driven by the invading water, the speed of which solely determines
the bank saturation (in the solution that conserves mass). The water causes the
discontinuity in the oil saturation that occurs across the rear of the oil bank. Therefore, the
discontinuity is not then, the result of the oil's governing equation, and only exists because
of the invading water. The KOPT model shows that where there is no water front, mass is
conserved with only a discontinuity at the oil front. The rear of the oil bank is of a
different character than the shock that is created due to the nonlinearity of the governing
equations (i.e., the oil front).
The bank saturations are determined by finding the saturation, S02» at which
Tl C S02 - Keo(S02,Sw) = TI C Sor = f(S0) (59)
where C is the speed (or celerity) of the water front. The last equation is equation 58
rearranged with S0i = Sor. This equation has a trivial solution at S02 = S0r. since
Keo(S0r»Sw) = 0. There may also be a solution for the oil bank: S0=S02>Sor. Figure 26
shows values of the function f(S0) for gasoline in a sandy clay loam soil (e = 8.8) and
various water front speeds. Bisection is used to locate the solution S02 in an interval from
the relative maximum of these curves, Si, to the maximum oil saturation, So(max)- The
saturation Si is found first by solving
a(r1CS0-Keo(So.Sw)) = Q
3S0
which has a relative maximum at
72
-------
f(So)
JV
40-
30-
20-
10-
-
-10-
-20-
-30-
i 1 1 1
s'" "\
/ \
» »
.' X " • '.
' x" N '
/ x' N \ •
/ x' ^ — ~ N » «
/' x' x \ \
/'." X- ^..--"" "\ \ V
x..---" \ »
\\ -
\
1 1 1 1
10 cm/day
100 cm/day
150 cm/day
200 cm/day
250 cm/day
0.0 0.2 0.4 0.6 0.8 1.0
So
Figure 26: Oil Bank Function
73
-------
The last equation is solved by bisection over the interval SOT < S i < So(max). By locating
the relative maximum at Si first, S02. is assured to be in the search interval for the oil bank
saturation determination. Since the oil bank function is double valued, this two step
procedure is necessary to have a reliable numerical method for determining the oil bank
saturation.
The implication of the solution presented in Figure 24 will be discussed
subsequently. Most strikingly, all of the mobile oil behind the water front is displaced by
the water. This is directly analogous to the steady state case, where the saturation S0i is
determined by the surface oil flux. Here, though, the flux is zero and the corresponding oil
saturation is S0r- In both cases the boundary condition controls the magnitude of the oil
saturation immediately behind the water front
Note that the gross features of this unsteady event are similar to the quasi-steady
case. For both there are water and oil fronts and an oil bank driven by the speed of the
water front. In the unsteady case, the oil ahead of the bank, as well as that in the bank
undergoes drainage. This involves the method of characteristics part of the solution, along
with the equation for the fronts, in determining the oil and water distribution.
In contrast to the steady state case, the oil saturation is only constant throughout the
bank under restricted conditions. If the saturation S0i and the speed dzwf/dt are constant,
then the bank saturation S02 is also constant (see the 6.0 and 6.1 hour profiles on Figure
27). This makes the saturations constant throughout the bank (S02 = S0s) as at each point
in time of the solution the same saturation is found in the oil bank. As long as the
conditions that created the bank remain the same, no mechanism exists for causing a change
in saturation in the bank. This occurs in the present case when the water front moves at a
constant speed. The solution of the model shows that variable saturations behind the bank,
S0i, do not occur. So the constant speed water front always gives a constant oil bank
saturation.
Alternatively, variable water front speeds, dzwf/dt, create unsteady oil banks
characterized by variable bank saturations (S02 * S03). Since the speed of the water front
slows as the water undergoes unsteady drainage, the oil saturations in the bank at the
location of the water front continuously are reduced. Continuously decreasing oil
74
-------
6.0 Hours
OH IT.
-10--
Depth
in cm
-20--
-30-r
0.0 O.S 1.0
Saturation
6.1 Hours
-10- •
-20--
-30
0.0 0.5 1.0
Saturation
10.0 Hours
OH CT 1
-10--
Depth
in en
-20--
-30
-+-
0.0 0.5 1.0
Saturation
15 Hours
-10-- :
-20-
-+-
0.0 0.5 1.0
Saturation
Hater
Oil
Figure 27: Profiles for Unsteady Flow Example
15
-------
saturations means that the method of characteristics solution holds. Therefore, oil
characteristics must originate along the waterfront.
As time goes on and the water front slows, each characteristic carries lower oil
saturations. The simple wave pattern emerges, but is not centered as for an undisturbed
drainage pattern. Profiles drawn through the bank during this time show that the bank
itself is undergoing drainage, illustrated by the profiles at 10 and 15 hours in Figure 27. In
general the bank saturations decrease with time as the oil is spread out at ever lower
saturations. Although the immediate effect of the infiltration of water is to increase oil
saturations locally (in the bank), the overall effect is to spread the oil out over a region at its
residual saturation quicker than the undisturbed drainage would.
Mathematical analyses of shock fronts (which are the analogous counterpart to
wetting fronts in gas dynamics) indicate that when characteristics are followed forward in
time, they intersect shock fronts (Rozdestvenskii and Janenko, 1980, and Lax, 1973).
When they are followed backward in time, they do not intersect shocks. The solution
presented above, apparently is in violation of these rules. Oil characteristics originate at the
water front; tracing them backward intersects the discontinuity in both oil and water
saturations. As discussed above, the rear of the oil bank is not a shock as it results from
the change in oil flow domain, caused by the water.
An important consideration, which has not yet been addressed is the mobility ratio
(KSW/KSO) of the liquids. For highly viscous oils, the maximum conductivity,
Keo(So(max))> can easily be less than values that would be required by the invading water.
Such cases represent unstable displacements, where the oil cannot be completely displaced
by the incoming water. Some amount of mobile oil remains behind the water front. In
these cases, "fingering" occurs as the water displaces some of the oil (Homsy, 1987, and
Kueper and Frind, 1988). This indicates that flow of the liquids is no longer kinematic and
that pressure gradients are important. The assumption that the water solution is
independent breaks down as all the saturations must be determined dynamically.
In Figure 26, the curves for the 200 and 250 cm/day water front speeds have no
nontrivial solution in the range of possible oil saturations. This can be seen in the figure,
since the solution occurs at the same function value as the S0 = Sor end point For the 250
cm/day speed the function value at Syntax). 22.094, is above the endpoint value, 8.200.
At this speed, then there is no value of the oil saturation that satisfies the jump conditions,
76
-------
and the kinematic model breaks down. The simulation ceases and the message "Unstable
displacement—model terminated" is output. The factors which control this condition are the
density and viscosity of the oil, the speed of the water front, and the pore size distribution
of the soil. A proposed method to extend the kinematic model to cases where the oil banks
cannot form under kinematic conditions is discussed in Weaver (1988).
KROPT Model Results
Example 4 Oil Spill in Coarse Sand-
The first 80 hours of the base characteristic plane for an oil spill is shown in Figure
28. No characteristics are drawn in the figure, but they exist in all of the oil banks and the
original drainage region. The input parameters can be found in Table 5. This example is
similar to the unsteady oil bank example, but is meant to illustrate the effect of multiple
rainfalls on the oil. Each rainfall that passes the rainfall event that is driving the oil bank,
causes the water front speed to increase (points A and B on Figure 30). Each time this
occurs a new oil bank is created, due to the abrupt increase in water front speed, and
corresponding jump in oil saturation (lines AD and BE). Since the banks move faster than
the oil front, they surpass the oil front, causing its speed to increase (points C, D and E).
Four oil profiles are shown in Figure 29. By 11 hours after the beginning of the oil
spill, the first rainfall has created the first of the oil banks, with a constant saturation of
0.8195. The original undisturbed oil drainage profile is below the bank. At 46 hours this
bank has long since passed the original oil front and replaced it The second rainfall has
begun and reached a depth of 14 cm. By 56 hours the second rainfall has passed the first,
triggering a new oil bank. The water front speed slows continuously from point A on, so
the second oil bank contains ever decreasing oil saturations. By 72 hours the third rainfall
has passed the second, triggering the third oil bank.
The entire simulation of 180 hours is shown in Figure 30. The last rainfall (line
EF) is of such a great magnitude, that it drives a very narrow oil bank ahead of it (Figure
31). After 137.9 hours the rainfall is spread out over 457.4 cm at its residual saturation of
0.05. The maximum oil volume balance error for the oil +0.29%, while the maximum
error for water is +0.09%.
77
-------
Depth -120--
(cm)
20 30 40 50 60
Time in Hours
80
Oil Front
Water Fronts
Oil Banks
Figure 28: Example 4 - Base Characteristic Plane
78
-------
Table 5: Examples 4, 5, and 6 Input Data
Soil <•)
•ater
Oil
u
c«/d«y
«r/c.5
fr/ca/day
gr/ca/d»y
*r/e»J
bill life dayt
lag UM days
loading »/day
«««* days
voluaw ea^/M*
depth c.
Conatilveat
C Unit) «g/l
kilt life da/a
lag tla* daya
Coarse Sand Sandy Clay Sandy Clay
loae, Leaa
720.0
C.5X)
O.J49
J.«
n/a
0.0500
0.2000
672
1285
0.73
0.05
2.0
0.147
lli.J
0.3*44
0.328
1.70
1
0.2787
0.4000
872
471
0.45
0.10
(e)
4.0
5.0
200
240 (b)
0.10 (b)
«0 .
lit.} (a)
0.3444 (a)
0.3900 (a)
0.2041
0.4240
(72
(b)
1000 (e)
0.65 (b)
0.10 (d)
360 (d)
0.3597
IS.2
(b)
Toluene r>rene
396. 16*9 a»pl« 6
Rainfalls
Intensl ty
c./hr
7-i
1.5
3.5
2.0
4.0
1.0
Intensity
c»/day
(eg In
hr
10.0
45.0
60.0
100.0
130.0
150.0
t«gini
oay
nlng tvllng
hr
13.0
49.0
41.0
104.0
135.0
152.0
nlng Ending
day
16.0
5.0
30.0
30.0
5.0
a) toll column co*pari>on
'• U.O
1- 18.0
3- 18.0
18.0
18.0
18.0
9.0
15.}
20.0
40.S
52.0
57.0
90.0
119.0
151.0
183.0
207.0
b) Randomly generated rainfall* lor Austin. T«
S«« llfuft 4.34 lor rainfall t«qu«nc«.
9,
15,
20.
40.
52.
17.
90.
119.
151.
161.
202.1
Temperature
Evaporation
Relative
timidity
25
0.41
20
0.32
21.4
0.1243
0.5
(b)
(d)
a) sell data fro. 8raJkciuUk «| al.. 1981
b) fro. Kyan «t al.. 1987
c) »•» text
d) value >«tiMted (or example 6
•) >«tntlal «»aporailon raiei calculated (roc Auitlo, n clleatlc
data
79
-------
Profile at 11 hr.
Profile at 56 hr.
«•
-50-
-100-
Depth
in cm
-150-
-200-
IT-
1
u-
-50-
-100-
-150-
-200-
-?«;n-
A .
~^^.^j .^
0.0 0.5 1.0
Saturation
Profile at 46 hr.
0.0 0.5 1.0
Saturation
Profile at 72 hr.
-50-
-100-
Depth
in cm
-150-
-200-
-250-
0.
Y\
]!
0 0.5 1.
u-
-50-
-100-
-150-
-200-
-250-
0 0.
: |.:
• — i
0 0.5 1.
Saturation Saturation
Water Hater
.. rti 1 ........ rtl 1
Figure 29: Example 4 - Profiles at 11, 46, 56, & 72 Hours
80
-------
50 100 150
Time in Hours
Oil Front
Water Fronts
Oil Banks
Figure 30: Example 4 • Base Characteristic Plane (continued)
81
-------
Profile at 120 hr.
Profile at 140 hr,
-450--
-500
0.0 0.5 1.0
Saturation
Water
Oil
-500
0.0 0.5 1.0
Saturation
Water
Oil
Figure 31: Example 4 - Profiles at 120 and 140 Hours
82
-------
In cases like this one, where a relatively small volume of oil is released, a series of
rainfall events may cause the oil to be immobilized over the depth given by equation 52. In
the kinematic model with no rainfalls, the oil will only be immobilized after an infinite
amount of time. The rainfalls may greatly reduce the time required for immobilization;
137.9 hours in this case.
Saturation histories for 50,100 and 200 cm are presented in Figure 32. At 50 cm
the first two oil banks are visible as an abrupt increase in oil saturation. Once the rainfall
events pass the history depth, the oil is left at its residual saturation for the remainder of the
simulation. The rainfalls can be seen at the deeper depths also and corresponding events
are marked on the 3 histories.
Example 5: Gasoline Leak into Sandy Clay Loam—
The next example is of a gasoline leak into a sandy clay loam soil (Table 5 and
Figures 33 through 36). The soil in this example is of a wider pore size distribution than
the sand in example 4, as indicated by the Brooks and Cory exponents - 6.6 for the sand
and 8.8 for the sandy clay loam. The water relative permeability curve for this example is
given in Figure 33. In this type of media the relative permeability drops very rapidly once
the saturation drops. After a rainfall event in such soil, drainage will appear to cease at
relatively high water saturations, because the relative permeability is so low. This
phenomena is responsible for the poorly defined concept of field capacity, which is defined
as the moisture content after gravity drainage has ceased. (The kinematic model, of course,
shows that gravity drainage never truly ceases.) In example 5 the influence of soil type
appears as relatively slowly moving fronts, with high saturations.
Figure 34 shows a 60 day simulation for infinite (nondegrading) and 360 day
gasoline half lives. Three of the rainfalls are intense enough to create oil banks (those
beginning at points A, C, and D), while others (points B, E, and F) remain behind the
deepest water front The gasoline is more mobile (Kso > KSW) so that it is difficult for the
water to go deep enough to bypass the gasoline. The gasoline can move as fast as required
by the water, in this case. (The possibility of unstable displacement remains, since the
initial infiltration capacity of the soil is infinite).
83
-------
History at 50 cm Depth
1.0
Saturation 0.5--
0.0
B C
50
100
150
1.0
1.0
Saturation 0.5--
0.0
History at 100 cm Depth
50
100
150
History at 200 cm Depth
D £
50 100 150
Time in Hours
Hater
Figure 32: Example 4 - Histories at 50, 100, and 200 cm Depths
84
-------
relative
permeability
0.0 0.2 0.4
0.8 1.0
Sw
Figure 33: Water Relative Permeability in a Sandy Clay Loam Soil
85
-------
E F
-500
10 20 30 40 50 60
Time in Days
Oil Front
Water Front
Oil bank
Oil Front 360 Day Half Life
Oil Bank 360 Day Half Life
Figure 34: Example 5 - Base Characteristic Plane
86
-------
Oil profiles shown in Figure 35 also reflect the high mobility of the gasoline,
because the oil saturations in the banks are lower than the water saturation in the fronts
driving them. The oil saturation jumps are of such a low magnitude that for all practical
purposes they are undetectable. The jumps are needed, however, to maintain the accuracy
of the computation. Somewhat unexpectedly, the oil banks persist over long periods of
time. The relative permeability function leads to low effective hydraulic conductivities as
the oil saturations become reduced as the banks move away from the water front. The
slowness of the banks means that more than one bank may be present at any given time.
As a result of this example, the code has been modified for this possibility.
When the oil decays, the features of the solution move slower as time goes on. The
reduction in effective hydraulic conductivity is evident in both Figures 34 and 35. The
profiles show decreased extents of gasoline, with small reductions in gasoline saturation.
Another feature, visible in the 40 and 60 day oil profiles, is that the water front always
leaves the oil at its residual saturation at the water front. This oil then degrades, so its
saturation is reduced further. In contrast to example 1, the locus of points at residual
saturation is defined mainly by the water front.
The extent of a dissolved constituent, with the properties of toluene (Ryan et al.,
1987) is indicated in Figure 36. The leading edge of the constituent generally mimics the
oil front, but is retarded due primarily to sorption. Only a small amount of this material
volatilizes. Concentration profiles are found in Figure 37a, and show increasing
concentration with time, due to oil biodegradation. For comparison, Figure 37b shows
profiles for the same initial concentration in gasoline of a constituent with benzene's
partitioning properties (ko = 134, kh = 0.15, ks = 6.0). The partitioning relationship
(equation 39) gives higher water concentrations for the benzene since it is less hydrophobic
than toluene. As a result the toluene can migrate deeper that the benzene, since it is carried
with the faster moving oil.
The examples presented above illustrate the fluid's migration where the initial water
saturation is a constant The equations for oil migration reduce to equation 21, since
dSw/dz and XQ are equal to zero. An attempt was made to simulate a case where there was
an initially uniform water profile in the soil, that began to evaporate before the oil was
released. Since 3Sw/9t is not now equal to zero, the oil saturations vary along the
characteristics, curving even when XQ is zero. Once out in front of the variable water
profile the oil characteristics become straight lines.
87
-------
Profile at 20 Days Profile at 40 Days Profile at 60 Days
Depth
in cm
V
-100
-200
-300-
-400-
-500-
^Y'-
)-
r ;
. — S
i
w
• -100
• -200-
• -300-
-400-
-500-
wate
-*>
.r
oi
.....
i>v-
\
\
\
I
....
a:
V
1
•
ir
-100
-200-
-300-
-400-
-500-
\
\
•
•
0.0 0.5 1.0 0.0 0.5 10 oTo 0.5 TTo
Saturation Saturation Saturation
Water
Oil no decay
Oil 360 Day Half Life
Figure 35: Example 5 - Profiles at 20, 40, and 60 Days
-------
-500
10 20 30 40 50 60
Time in Days
Constituent Leading Edge
Volatile Front
Figure 36: Example 5 - Constituent Distribution
89
-------
(a)
(b)
Depth
in cm
-200--
-400--
Concentration in Water
mg/1
20 Days
40 Days
. 60 Days
u-
200-
400-
1
i
•j "
Concentration in Water
mg/1
20 Days
40 Days
60 Days
Figure 37: Example 5 • Concentration Profiles, (a) toluene and (b) benzene
This presents some extreme difficulties for the model. First, the water profile has
an infinite slope at the ground surface. This causes a problem for the oil characteristics that
originate at the surface during the oil release. They face an infinitely large decrease in
saturation, resulting from equation 20, in an infinitesimally small depth. The RKF1(2)
solver uses function evaluations at the surface, and two depths for the first time step (for
each characteristic). Error must be introduced into the solution to avoid the saturations
becoming large negative numbers as the solution starts. In a fine sand the error was
approximately 5%, while in a silty sand the error was about 50%. For this reason this case
was not pursued any further.
When water enters a soil containing variable water saturations, the water front
speed increases at first, since the water saturations ahead of the front increases. In this
case, the oil bank saturations increase as time goes on, instead of remaining the same or
decreasing as in the previous examples. The speed of each emerging oil characteristic is
90
-------
higher than that of the previous characteristic. The characteristics could cross, creating a
new discontinuity of some type. It was found, however, that the 3Keo/dSw term in
equation 20 acts in such a way to reduce saturations so that the characteristics do not cross.
Variable water saturations make the solution very complicated from the point of
view of managing the characteristics. Not only must the curving characteristics be
accounted for, but also the equations for the changing oil saturations. There may be several
sets of characteristics, associated with different oil banks. In the face of the complexity of
the solution, using finite difference techniques for the solution of the characteristic
equations appear to be an attractive alternative. Little efficiency would be lost, since the
calculation of the characteristics is no longer insignificant
Example 6 Land Treatment-
Before presenting the land treatment example, some features of the KROPT code
specifically designed for land treatment will be discussed. Because of the unresolved
limitations of the kinematic model, the land treatment features of the model were built
around an immobile oil phase. Oil applied at saturations slightly above residual will not
move very far (see example 3 and equation 52). The oil's low flux contributes little to the
constituent migration. For a land treatment case the oil, constituent and applied water are
incorporated over some depth (the plow zone) into the soil. Multiple waste loadings,
irrigation to a specified moisture content and cultivation are capabilities of the model.
The KROPT model requires the initial concentration in oil to be input This can be
calculated from the concentration in the waste according to the mass balance
CwasteVwaSte = CoCVo + Vw/ko) . (62)
This implies that
r1 — r , _r \ - CwasteVwaste ,,-o\
Co - ^surface) - VQ + Vw/ko . (63)
Q>(surface) is then used as input in equation 40 to give the initial concentration in the oil.
91
-------
The land treatment example is based on soil column experiments performed by
Ryan et al. (1987). Input data is found in Table 5. Since the experiments were not
designed with the KROPT model in mind, some of the input parameters must be estimated
(indicated in Table 5). The experiments were run for 6 months, at the end of which time the
columns were destructively sampled for a variety of waste constituents. Toluene and
pyrene were selected for comparison of the experiments with the KROPT model. The
columns were cultivated every two weeks and checked to assure proper moisture content
(80% field capacity —> Sw = 0.4240). Approximately 5 centimeters per month of water
was added to the columns to simulate rainfall. The waste was composed of 18.9% oil and
grease, 66.6% water, and 14.5% solids by weight. The effect of the solids was neglected
in the model. The column simulated had a waste loading that was 8% by weight of the
total of soil and waste.
The experimental results showed no migration of any constituent out of the bottom
of the plow zone; except for a few measurements judged anomalous by Ryan et al. (1987).
This is close to the result shown in Figure 38 for the soil columns. The toluene and pyrene
migrate a small distance (approximately 3 cm over 6 months), because of the high water
fluxes from individual rainfall events. The result is due to the constituents being
preferentially partitioned to the immobile oil and the low overall water flux. The jumps
visible in the leading edge are associated with the arrival of an intense rainfall. The
simulated mass balance of the water is shown in Figure 39. A series of rainfalls
approximated the times, fluxes and volumes of water added in the experiment. Less water
was added for irrigation in the simulation, than in the experiment
Figure 40 shows the loss of pyrene and toluene from the columns with time.
Virtually complete removal is achieved, while the experimental results indicate 85%
removal of pyrene and removal of toluene to below the detection limit Since pyrene is
modeled as a nonvolatile constituent, its loss is due solely to its biodegradation half life; the
average measured by Ryan et al. (1987) of 35 days was used in the simulation. Toluene,
in contrast is lost only through volatilization. The top of the volatile front is indicated in
Figure 38. At each cultivation the volatile front is brought back to the surface, causing the
rate volatilization to increase temporarily. In the model, irrigation and cultivation occur at
the same frequency (14 days), unlike the experiment
Simulation of the same loading of waste and toluene was performed for a series of
rainfall events generated by the methods of Charbeneau et al. (1988) for Austin, Texas
93
-------
-6--
-8--
Depth 1Q.
(cm) 1
-12--
-14--
-16--
-18--
-20-1
50 100
Time in Days
150
Pyrene Depth {soil columns)
Toluene Depth (soil columns)
Toluene Volatile Front (soil columns)
Toluene Depth (random rainfall)
Toluene Volatile Front (random rainfall
Figure 38: Example 6 - Base Characteristic Plane
94
-------
50
45--
40--
35--
30--
cm3/cm2 25--
20--
15-
10
5-
_ r
_ j
50 100
Time in Days
150
Cumulative Infiltration (soil columns)
Cumulative Evaporation (soil columns)
Cumulative Infiltration (random rainfal
Cumulative Evaporation (random rainfall
Figure 39: Example 6 - Water Balance
95
-------
0.26
0.24--
0.22-J
0.20--
0.18--
0.16--
Mass
Mg/cm2
0.00
100
Time in Days
150
Pyrene Mass (soil columns)
Toluene Mass (soil columns)
Toluene Mass (random rainfall)
Figure 40: Example 6 - Constituent Mass
96
-------
climatic conditions. Similar results were obtained for the depth of the toluene determined in
the soil column simulation. More water, however, was added to the profile by the random
rainfall generation, suggesting that the transport of hydrophobic constituents is relatively
insensitive to the rainfalls.
KROPT Model Implementation
The implementation of the full solution in the KROPT code is much more difficult
than for the simplified KOPT model. The following features must be included in KROPT,
• for water
- solution of the water equation for multiple rainfall events,
- evaporation of the current rainfall,
- Green-Ampt flow of the water during the rain,
• for oil:
- creation of the oil banks,
- motion of the oil banks,
• for the constituent:
- biodegradation of oil with varying water saturations,
- volatilization with variable evaporation rates,
• for land treatment cases:
- irrigation of the land treatment site for maintaining mininmrn
moisture content,
. cultivation of the site, and
- multiple waste loadings at specified intervals.
Table 6 contains a synopsis of the programming details used in the KROPT code,
for each of the above items. Included are the number of fronts and characteristics used in
the program and the treatment in the controller routine. Since more characteristics are
needed than for the KOPT model, more attention must be paid to the data structures. The
97
-------
Table 6: KROPT Programming Synopsis
tqmtion
front*
Characterit»ticj
Motet
W»t«r
10
Oil
Mondecaying
becayin?
Oil bank
Constituent
nonvolatile
volatile
ttOB
1
1
I
0
1
Control ler
0
SI
Si/bank
11(22)
11(22)
Checks
Vp to SI characteristic origins saved per bank
10 characteristics used during oil release
41 characteristics used for simple wave
Maximum characteristics used per bank
22 Characteristics (or Land Treatment
front is for the volatilisation model
Hetes
Water
Oil
Beginning of Event
End o( tvent
tad of S
v
-------
underlying philosophy in this implementation is to store only the data that must be used
again during a time step, and as in KOPT, to store data only for one time step.
The solution algorithm proceeds through the following order: water, oil and then
constituent equations are evaluated. Figure 6 shows the subroutines used in obtaining the
function evaluations for each partial step taken by the RKF solver (boxes A, B and C of
Figure 6). Next to each subroutine is a description of the operations performed in it The
controller structure is also indicated in the figure.
Integration of the profile is done in KROPT by breaking up the profile into intervals
and using the Gaussian quadrature routine to integrate each interval of the profile. The
breaks are picked so that each discontinuity (oil front, oil banks, water fronts, leading and
trailing edges of the constituent) fall on the end points of an interval. An advantage of the
quadrature routine is that the function is not evaluated at the endpoints of the intervals,
o
which is where the saturations and concentrations are double valued. A result of this
procedure is that the number of function evaluations for the profile varies throughout the
solution. The accuracy of the mass balance computation varies, then, throughout the
solution. Some cases were found to have errors that, in reality, resulted from the profile
integration instead of error in the solution.
COMPUTATIONAL PERFORMANCE
A summary of the code performance is presented in Table 7 for examples 1 through
6. In the first column is the example number and the duration of the simulation, whatever
time unit used. Information concerning the RKF1(2) solver steps taken is indicative of the
computational effort needed to obtain the solution. In columns 2 through 5 are tabulated
the total number of steps taken, the percentage of steps that were rejected, the average
number of function evaluations per step, and the average number of model equations used
for each evaluation. The CPU time in seconds is reported in column 6 and the time per
RKF step in column 7. The runs were made on a Microvax II. Each time listed is the
average of three runs made when there were no other users and includes all processing and
the I/O. Mass balance computations were performed for each accepted step of the solution,
adding, not insignificantly, to the computational burden.
99
-------
Table 7: Solution Performance for Examples 1-6
(I
Ou
IB
Ik
2«
2k
2e
3d
2*
31
1
4
SB
5b
Sc
3d
5*
if
«B
4k
4e
(Bl
• BBl»
(1*0)
(1*01
(5451
(5431
(545)
(545)
(545)
(7)01
(1101
(el
(el
(el
(c)
1401 (c)
(401 (c)
KOI le|
(40) (e|
(2101
(2101
12101
dDTBtiea
HMO) St.pi
l«1«ct«d step
212
251
2*4
331
313
)47
112
411
(70
455
541
43*
71*
0)5
4)4
*14
1)4
1757
42)*
4.44
12.25
11.4)
11.2)
».;«
».42
».oo
7.34
10.)*
7.1)
10.15
11.41
11.4*
11. »4
10.03
4.1*
1 .12
•1 BlBBlBtiBB IB
2.1)
2.1)
2.24
3.3*
2.)*
2.21
2.21
2.27
2.2)
2.24
2.
2.
2.
2.
2.
2.
2.
2.
25
40
42
41
4)
34
14
10
ctu ti
txluition
1
10
1
12
*
33
44
10
>
.00
.40
.00
.00
.72
.14
.24
.70
.5)
13^54
15.75
75.74
0*.ll
25.))
24.44
27.25
4.70
12.14
*.74
24.14
25.41
24.7)
44.02
42.04
10.45
11. l»
74.02
131.10
144.47
154.34
203.52
1)1.5*
273.27
5X.13
a* in i*e
Jt.p
0.0)14
0.0512
0.0331
0.0774
0.0445
0.0474
0.1152
0.102)
0.0*25
0.}2*1
0.13*1
0.2033
0.20*3
0.1147
0.2433
0.1574
0.1547
0.124*
•ot«»
(kl
ad.ac
d.ac
Bd.ee
ad.e
ad.ve
d!c
d.»e
d.»e
ad.ae
ad(ae
Bd.e
Bd.»e
d,ae
d.e
d,»e
c
»e
•e.rr
whBt««*r tia» aalt u«»d
(kl •bbt«»i«ti»o§: c, vita a«av*l*til* c*a«tita»t; d, d»c«riaf •!!: ae. with a* eaa>tit«*ati
ad, a«ad»c«riat oil; cf r«ad»aly a*a*r«t*d riiaftll •voati; »e, with rolttil* c*a-
•titu«at
Icl r«iult> aot pr»*Bt*d ia th* t»t
100
-------
Example 2 was selected for a comparison of various KOPT model simulations.
The first three (2a - 2c) have nondecaying oil and thus the fewest number of equations to
solve. For a nondegrading oil with no dissolved constituent exactly 1 equation is solved
per step. Fairly large steps were taken, which required a high percentage of the steps to be
rejected. The dissolved constituent adds exactly 11 equations to the simulation (example
2b). Smaller steps were taken in order to satisfy the truncation error criteria on all 12
equations. When the constituent is volatile (2c) another equation is added for the
volatilization model. The average number of model equations is lower than in case (2b),
since all of the constituent is volatilized before the end of the simulation. When the oil is
degradable (2d), the number of equations jumps, because the characteristic equations are
solved to find the characteristic locations. The average number is less than the total allowed
(51), because some characteristics are dropped from the solution set as time goes on. The
number of steps increases, again for the reason that more equations must not exceed the
truncation error criteria. The same trend continues when the dissolved constituent is added
(2e and 2f). The CPU times for these cases increase as expected as the complexity of the
simulation increases.
A similar comparison was made for Example 5, which gave in general similar
results. The typically higher number of equations results from the various water fronts,
and the oil bank. The increase from 5a to 5d is higher than for the KOPT code (2a to 2d),
because more oil characteristics are needed to simulate the decaying oil bank. One notable
fact is that when no mass balance computations were performed for example 5f, the
average CPU time dropped to 180.89 seconds. The increase in efficiency is obtained at the
cost of knowing how accurate the solution is.
In the land treatment cases, no oil equations are solved. The computation time for
the volatile constituent is much higher than that for the nonvolatile one, due to the high
volatile fluxes that occur each time the plow zone is cultivated. More time is needed for the
randomly generated rainfalls, because of the greater number of rainfall events and the
occasionally higher fluxes than in the controlled column study.
101
-------
SECTION 7
DISCUSSION AND CONCLUSIONS
This work has shown that kinematic theory can be applied to multiphase flow
problems. Kinematic models have an advantage over traditional numeric models in that the
kinematic models are well-suited for studying the hyperbolic or wave behavior of the
multiphase flow equations. Inclusion of sharp fronts is natural to the kinematic models,
while finite difference or finite element models encounter difficulties when derivatives in
the governing equations tend to become undefined. Solution of the kinematic model shows
how the oil and water fronts interact Thus by simplifying the governing equations, some
of the fundamental behavior of the solution is clearly revealed. The effect of the neglected
capillary pressure gradients are to increase the infiltration capacity of the soil (accounted for
by Charbeneau et al., 1988) and smear the fronts. Much of the important physics is
retained, however, in the kinematic model
Computationally, by reducing the original system of coupled partial differential
equations to a system of ordinary differential equations, the difficulty of finding a
numerical solution is reduced greatly. Resulting computer codes are well suited to solving
problems where there is a constant amount of water in the pore space. Such a situation
corresponds to a soil at the so-called field capacity and can represent typical conditions of
the profile. The best performance of the KOPT model is obtained when the oil is
nondegrading and the oil characteristics are straight. Curving characteristics, due to
biodegradation, increase the computational effort
The implementation of the more general model (KROPT) reveals the interactions
between the oil and water fronts. The model is limited to situations where the initial water
profile is uniform. Once again the best performance is obtained when the oil characteristics
are straight, since the characteristic equations are not solved explicitly. Performance is
degraded where rainfalls exceed the kinematic capacity of the soil. Since the infiltration
capacity is infinite initially, very large water fluxes are possible until runoff is produced.
During this period small time steps are required.
Two categories of conclusions will be discussed:
102
-------
• physical phenomena occurring in multiphase flow,
• numerics of the methods selected for the solution.
PHYSICAL PHENOMENA OCCURRING IN MULTIPHASE FLOW
This modeling work enhances understanding of multiphase flow through the results
obtained from the example problems presented in Section 6. Most important is the material
concerning the oil banks. If incoming rainfalls displace oil, then the banks are created.
Although this is a result of simplified modeling, the notion that incoming water displaces
oil into a bank, moving ahead of the water front, is a more general result In situations
where the flow of the fluids is not kinematic, banks will also exist In certain cases, the
banks may be very pronounced, although short lived. In others the banks may persist for
longer periods, but would be virtually undetectable in the laboratory and the field or by
numerical methods that are not based upon the kinematic model.
The kinematic model results indicate that the water displaces all of the mobile oil.
This result was contrary to what was expected, but was established through conservation
of mass. Some displacement of oil would be expected, along with the creation of an oil
bank, but the model result may or may not be confirmed by experiment
The value of S0r is critical for the results obtained from the model This parameter
could be easily manipulated to show that oil would spread over a smaller depth (see
equation 52).
Biodegradation of the oil has only a small effect upon the ultimate depth attained by
the oil. It does effect the mobility of the oil (reducing it). An interesting effect on the
constituent is caused by biodegradation of the oil. Concentrations increase as the oil
degrades because of the loss of solvent If the constituent is also biodegradable, the shorter
half life determines if the concentration will increase or not
As in field cases, where constituent concentrations in water are measured, a large
mass of constituent may reside in the oil phase and not be detected. In example 2
concentrations in the oil are 50 times greater than in the water. If only the water phase
concentrations are measured, most of mass in the soil is missed. This is particularly a
problem if the oil is immobilized, since the oil cannot flow into an observation well
103
-------
(phenomena akin to end effects may prevent oil flow into wells even at higher saturations).
The presence of the oily material is not detected by the observation wells; the core material
itself must be analyzed
NUMERICS OF THE METHODS SELECTED FOR THE SOLUTION
Using a differential equation solver is well suited to solving the kinematic model,
since the model equations can be written as a system of ordinary differential equations.
The speed of the solver is degraded as more and more equations are added. For these
situations, it might be advantageous to solve the method of characteristics portion of the
solution by finite difference techniques instead of tracking the characteristics. The
overhead associated with using the RKF solver may be comparable to that for the
discretization methods. Also, the necessity to track sets of oil bank characteristics makes
the program exceedingly complex. Extension of the model as discussed below may not be
practical, because of this.
Severe difficulties are encountered in solving the oil equations when the initial water
profile is not constant This occurs because the slope of the water profile is infinite at the
surface. The term in the method of characteristics solution that gives the change in oil
saturation (equation 20) has an infinite negative slope as the solution commences.
Modification of the model so as to ignore this allows a solution to be obtained, but that
solution does not conserves mass.
104
-------
REFERENCES
Abramowitz, M. and I. A. Stegun, 1965, Handbook of Mathematical Functions. Dover,
New York.
Abriola, L. M. and G. F. Finder, 1985a, A multiphase approach to the modeling of porous
media contamination by organic compounds, 1 Equation development, Water Resources
Research. 21. 11-18.
Abriola, L. M. and G. F. Finder, 19855, A multiphase approach to the modeling of porous
media contamination by organic compounds, 2 Numerical simulation, Water Resources
Research. 21,19-26.
Bardon, C. and D. G. Longeron, 1980, Influence of very low interfacial tensions on
relative permeability, Society of Petroleum Engineers Journal. 269.
Bear, J., 1972, Dynamics of Fluids in Porous Media. American Elsevier, New York.
Bouwer, H., 1966, Rapid field measurement of air entry value and hydraulic conductivity
of soil as significant parameters in flow system analysis, Water Resources Research. 2,
729-738.
Brakensiek, D. L., R. L. Engleman, and W. J. Rawls, 1981, Variation within texture
classes of soil water parameters, Transactions of the American Society of Agricultural
Engineers. 335-339.
Brooks, RJL, and A.T. Corey, 1964, Hydraulic Properties of Porous Media, Colorado
State University Hydrology Paper No. 3, Fort Collins, Colorado.
Burdine, N. T., 1953, Relative permeability calculations from pore size distribution data,
Transactions American Institute of Mining Engineers. 198,71-78.
Charbeneau, R.J., 1984, Kinematic models for soil moisture and solute transport, Water
Resources Research. 20(6), 699-706.
Charbeneau, R. J., 1988, Liquid Moisture Redistribution: Hydrologic Simulation and
Spatial Variability, Proceedings of the NATO Workshop on Modeling Flow Processes in
the Unsaturated Zone, in press.
Charbeneau, R. J. and R. G. Asgian, 1988, A model for the simulation of the transient
soil water content profile, submitted.
Chow, V. T., D. R. Maidment, L. W. Mays, 1988, Applied Hydrology. McGraw Hill,
New York.
Corapcioglu, M. Y. and A. L. Baehr, 1987, A compositional multiphase model for
ground water contamination by petroleum products: 1. Theoretical considerations, Water
Resources Research. 23,191-200.
105
-------
Delshad, M. and G. A. Pope, 1987, Comparison of the three-phase oil relative
permeability models, submitted to Journal of Transport in Porous Media.
Dracos, T., 1978, Theoretical considerations and practical implications on the infiltration of
hydrocarbons in aquifers, Proceedings International Symposium on Groundwater Pollution
by Oil Hydrocarbons, Prague, June 5-9,1978, International Association of Hydrologists.
Faust, C. R., 1985, Transport of immiscible fluids within and below the unsaturated zone -
a numerical model, Water Resources Research. 21,587-596.
Payers, F. J. and J. D. Matthews, 1984, Evaluation of normalized Stone's methods for
estimating three-phase relative permeabilities, Society of Petroleum Engineers Journal. 24,
224-232.
Fehlberg, E., 1969, Low-order Classical Runge-Kutta Formulas With Stepsize Control
and Their Application to Some Heat Transfer Problems, NASA TR R-315.
Forsythe G. E., M. A. Malcolm, and C. B. Moler, 1977, Computer Methods for
Mathematical Computations. Prentice-Hall, Englewood Cliffs, NJ.
Green, W.H., and G.A. Ampt, 1911, Studies on soil physics, Part I, The flow of air and
water through soils, J. Agnc. Sci.. 4(1), 1-24.
Greenberg, M. D., 1978, Foundations of Applied Mathematics. Prentice-Hall, Englewood
Cliffs, New Jersey.
Helfferich, F. G., 198\i Theory of multicomponent, multiphase displacement in porous
media. Society of Petroleum Engineers Journal. 21.51-62.
Hochmuth, D. P. and D. K. Sunada, 1985, Ground-water model of two-phase immiscible
flow in coarse materials. Groundwater. 23.617-626
Holzer, T. L., 1976, Application of groundwater flow theory to a subsurface oil spill,
Groundwater. 14, 138-145.
Homsy, G. M., 1987, Viscous fingering in porous media, Annual Reviews of Fluid
Mechanics. 19.271-311.
Jeffrey, A. and T. Taniuti, 1964, Non-linear Wave Propagation. Academic, Orlando, Fla.
Kreyszig, E., 1972, Advanced Engineering Mathematics. 3rd ed., Wiley, New York.
Kueper, B. H. and E. O. Frind, 1988, An overview of immiscible fingering in porous
media, Journal of Contaminant Hydrology. 2,95-110.
Kuppusamy, T., J. Sheng, J. C. Parker, and R. J. Lenhard, 1987, Finite-element analysis
of multiphase immiscible flow through soils, Water Resources Research. 23,625-631.
Lax, P. D., 1973, Hyperbolic Systems of Conservation Laws and the Mathematical Theory
of Shock Waves, Society for Industrial and Applied Mathematics. Philadelphia.
Leverett and Lewis, 1941, Steady flow of gas-oil-water mixtures through unconsolidated
sands, Transactions. American Institute of Mining Engineer^. Vol. 142, T.P 1206 pp
107-152.
106
-------
Lighthill, M. J. and G. B. Whitham, 1955, On kinematic waves, 1, Flood movement in
long rivers, Proc. R. Soc. London. Ser. A. 229, 281-316.
Mein, R.G., and C.L. Larson, 1973, Modeling infiltration during a steady rain, Water
Resources Research. 9(2), 384-394.
Mull, R., 1970, The migration of oil-products in the subsoil with regard to ground-water
pollution by oil. Proceedings of the Fifth International Conference on Advances in Water
Pollution Research, San Francisco and Hawaii, S. A. Jenkins ed., Pergammon, Oxford,
vol. 2, HA 7(A)/l-8.
Mull, R., 1978, Calculations and experimental investigations of the migration of oil
products in natural soils, Proceedings International Symposium on Grpundwater Pollution
by Oil Hydrocarbons, Prague, June 5-9,1978, International Association of Hydrologists,
167-181.
Osbourne, M. and J. Sykes, 1986, Numerical modelling of immiscible organic transport at
the Hyde Park landfill, Water Resources Research. 22,25-33.
Parker, J. C., R. J. Lenhard, and T. Kuppusamy, 1987, A parametric model for
constitutive properties governing multiphase flow in porous media, Water Resources
Research. 23, 618- 624.
Parker, J. C. and R. J. Lenhard, 1987, A model for hysteretic constitutive relations
governing multiphase flow: 1. Saturation-pressure relations, Water Resources Research.
23, 2187-2196.
Peaceman, D. W., 1977, Fundamentals of Numerical Reservoir Simulation. Elsevier
Scientific, Amsterdam.
Philip, J.R., 1969, Theory of Infiltration, in Advances in Hydroscience. ed. by V.T.
Chow, vol. 5, pp. 215-296.
Raisbeck, J. M. and M. F. Mohtadi, 1974, The environmental impacts of oil spills on land
in the arctic regions, Water. Air and Soil Pollution. 3,195-208.
Rawls, W.J., D.L. Brakensiek, and N. Miller, 1983, Green-Ampt Infiltration Parameters
from Soils Data, J. Hydraulic Div.. ASCE. vol. 109, no. 1, pp.62-70.
Richards, L. A., 1931, Capillary conduction of liquids through porous mediums, Phvsics.
1, 318-333.
Roberts-Schultz, S., 1988, Immiscible-phase Flow in the Subsurface and Implications for
Facilitated Transport, Masters Thesis, University of Texas at Austin.
Rozdestvenskii, B. L. and N. N. Janenko, 1980, Systems of Ouasilinear Eqnations and
Their Dynamics. American Mathematical Society, Providence.
Ryan, J. R., M. L. Hanson, and R. C. Loehr, 1986, Land treatment practices in the
petroleum industry, Land Treatment: A Hazardous Waste Management Alternative. R. C.
Loehr and J. F. Malina eds., The University of Texas at Austin, 1986.
107
-------
Ryan, J., R. Loehr, and R. Sims, 1987, The Land Treatability of Appendix VIII
Constituent Present in Petroleum Refinery Wastes: Laboratory and Modeling Studies,
draft, American Petroleum Institute.
Saraf, D. N. and F. G. McCaffery, 1981, Two- and Three-Phase Relative Permeabilities:
A Review. Petroleum Recovery Institute, Report 81-8.
Schwille, F., 1988, Dense Chlorinated Solvents in Porous and Fractured Media. Lewis,
Chelsea, Michigan, 1988.
Short, TJ2., 1986, Modeling of Processes in the Unsaturated Zone, in Land Treatment A
Hazardous Waste Management Alternative. R.C Loehr and J.F. Malina Jr. eds., Water
Resources Symposium Number Thirteen, Center for Research in Water Resources,
University of Texas, Austin, Texas, pp. 211-240.
Sisson, J. B., A. H. Ferguson, and M. T. van Genuchten, 1980, Simple method for
predicting drainage from field plots, Soil Science Society of Americal Journal. 44,1147-
1152, 1980.
Smith, R. E., 1983, Approximate soil water movement by kinematic characteristics, Soil
Science Society of America Journal. 47,3-8.
Stone, H. L., 1970, Probability model for estimating three-phase relative permeability,
Journal of Petroleum Technology. 20,214-218.
Stone, H. L., 1973, Estimation of three-phase relative permeability and residual oil data,
The Journal of Canadian Petroleum Technology. 12,53-61.
Stumm, W. and J. J. Morgan, 1981, Aquatic Chemistry. Wiley, New York.
Thibodeaux, L. J. and S. T. Hwang, 1982, Landfarming of petroleum wastes - modeling
of the air emissions problem, Environmental Progress. 1,42-46.
Weaver, J.W., 1988, Kinematic Modeling of Multiphase Subsurface Contaminants,
Dissertation, University of Texas, Austin.
Weaver, J. W., 1984, A Critical Review of Immiscible Groundwater Pollutant Transport,
M.S. Thesis, The University of Texas at Austin.
Whitham, G. B., 1974, Linear and Nonlinear Waves. Wiley-Interscience, New York.
Wyllie, M. R. J. and G. H. F. Gardner, 1958, The generalized Kozeny-Carman equation:
Its application to problems of multiphase flow in porous media, World Oil. 146,210-227.
Youngs, E. G. and A. J. Peck, 1964, Moisture profile development and air compression
during water uptake by bounded porous bodies: 1. Theoretical introduction, Soil Science
98, 290-294.
108
-------
APPENDIX A
KOPT USERS GUIDE
109
-------
PROGRAM ICOPT USERS GUIDE
James W. Weaver
Randall J. Charbeneau
Environmental and Water Resources Engineering
Department of Qvil Engineering
The University of Texas at Austin
The program KOPT (Kinematic Oily Pollutant Transport) was developed for
simulation of the fate and multiphase transport of hydrocarbons within the vadose zone.
These hydrocarbons may have been released through leaking tanks and spills or applied
during land treatment of organic wastes. KOPT applies kinematic modeling theory to
simulate the migration of a hydrocarbon phase (referred to as an "oil" ) and a dissolved
constituent as they migrate through the vadose zone under the force of gravity. Water
movement is considered to occur at the average annual infiltration rate. The dissolved
constituent is free to partition within the oil-water-air-soil system according to linear local
equilibrium assumptions. The model includes the processes of volatilization, sorption of
the dissolved constituent, and degradation of both the dissolved constituent and the oil.
The simulation model provides the concentration and hydrocarbon profiles at any time as
well as the hydrocarbon and constituent flux past any horizon as a function of time. KOPT
was developed under funding from both the U.S. Environmental Protection Agency and
the Texas Water Development Board. It is written in Standard Fortran 77 and will run on
PC or mainframe computers.
The Input Data File
The input data for KOPT is contained in a file called "INPUT.KOP" which is free
format and is opened and read by the program. This file must be developed by the user and
stored on the same directory as the executable code. The program creates a file
"OUTPUT.KOP" and saves the program output on it. The output file can be edited or
printed. Each card of the input file is described below.
Card 1-3. Run Title
RUN TITLE (3 LINES OF UP TO 50 CHARACTERS EACH)
Card 4. Matrix Properties
WKS, XLAMB, ETA, RHOS, SWR
WKS saturated hydraulic conductivity (water)
XLAMB pore size distribution index (Brooks & Corey's 1)
ETA porosity
RHOS bulk density of matrix
SWR residual water saturation
110
-------
Card 5. Water Properties
WMU, WRHO, IRT, QW or SWMAX
WMU dynamic viscosity of water
WRHO density of water
IRT rainfall input type: l=flux specified
2=saturation specified
QW or SWMAX constant water flux or saturation
Card 6. Oil Characteristics
PMU, PRHO, SPR, HLO, OLAG, IAT
PMU dynamic viscosity of oil
PRHO oil density
SPR residual (trapped) oil saturation
HLO oil half-life (=0. if non-degrading)
OLAG acclimation rime for biodegradation to begin
IAT oil input type l=flux specified
2=volume/area specified
3=constant head
Card?. (forIAT=l) Oil Flux
QP, TPB, TPE, DPA, IP, HS
QP oil flux (volume/area/time)
TPB oil event beginning time
TPE oil event ending time
DPA depth at which pollutants are applied
IP ponding type 0=excess oil runs off surface
l=variable (see note 4.)
3=constant
HS constant head for IP=3 cases
Card?. (for!AT=2) Oil Volume
PVOL,DPA,DPL
PVOL oil volume/area incorporated into the soil
DPA depth at which pollutants are applied
DPL lower depth of initially polluted zone
Card 7. (forIAT=3) Constant Head
QP, TPB, TPE, DPA, IP, HS
Same as for IAT=1 except QP is meaningless
111
-------
Card 8. Dissolved Constituent
COINI, FCV, HLC, CLAG
COINI initial concentration in oil (see note 5.)
FCV factor for consistent volume units (see note 5.)
HLC constituent half-life (=0. if non-degrading)
CLAG acclimation time for biodegradation to begin
Card 9. Partition Coefficients
XXKO, XXKV, XKOC
XXKO oil-water partition coefficient
XXKV air-water partition coef. (Henry's law)
XKOC soil-water partition coefficient
Card 10. Volatilization Model
DAIR, DWV, EVAP, TEMP, RH, XCF
DAIR diffusion coef. for constituent in air
DWV diffusion coef. for water vapor
EVAP water evaporation rate (vol7area/time)
TEMP temperature in degrees C
RH relative humidity
XCF conversion factor to gm/cc (see note 2.)
Card 11. Simulation Parameters
TM.DM.DTPR
TM simulation ending time
DM maximum solution time step
DTPR minimum time between printed time steps and mass balance checks
Card 12. Profiles and Histories
~~ NTIMES.NZS
NTIMES number of profiles desired (up to 10)
NZS number of histories desired (up to 3)
Card 13. Profiles
PR(NTIMES) profile times
112
-------
Card 14. Histories
HI(NZS) history depths
Notes
1. All variable names are in accordance with fortran naming - names beginning with I
through M are integers, all others are reals.
2. With the following exceptions, all input parameters must be given in a consistent
system of units.
The temperature must be given in degrees C.
The conversion factor from density units to gm/cc (f/XCF = gm/cc;
volatilization model uses an empirical formula for water vapor density).
3. Zeros should be read in fields pertaining to unused values.
4. Variable surface ponding is supported only for cases with no dissolved constituent
present
5. The concentration of the dissolved constituent can be read in any desired units (eg.
mg/L). A conversion to the consistent volume unit must be read in variable FCV in order
for the calculated mass balances to be correct For example, with concentration in mg/L =
gm/m3 and C-G-S units, then the flux must be multiplied by the number of m2/cm2 and
FCV = 0.0001. The output concentration will be in mg/L = grn/m3 while the mass balance
will be in mg/cm2.
6. Detailed information on the program KOPT is contained within the Ph.D.
dissertation of Weaver (1988).
References
Weaver, J.W., 'Kinematic modeling of multiphase subsurface transport', Ph.D.
Dissertation, The University of Texas at Austin, August 1988.
113
-------
APPENDIX B
KOPT PROGRAM LISTING
114
-------
PROGRAM KOPT
C
C
C * *
C * KOPT KINEMATIC OILY POLLUTANT TRANSPORT *
C * *
C * *
C * PURPOSE THIS PROGRAM SIMULATES THE GRAVITY DOMINATED *
C * ONE DIMENSIONAL DOWNWARD TRANSPORT OF AN OILY GROUND- ' *
C * WATER POLLUTANT. THE OIL MAY CONTAIN A DISSOLVED *
C * CONSTITUENT. THE MIGRATION OF THESE OCCURS IN THE *
C * PRESENCE OF A CONSTANT AMOUNT OF PORE WATER. THE CON- *
C * STITUENT IS SUBJECT TO PARITIONING BETWEEN THE OIL, *
C * WATER, SOIL AIR, AND THE SOIL. BOTH THE OIL AND THE *
C * CONSTIUENT MAY BIODEGRADE OCCORDING TO FIRST ORDER *
C * KINETICS.
C * *
C * THE SOLUTION IS OBTAINED VIA A RUNGE-KUTTA-FEHLBERG IN- *
C * TEGRATION OF A SYSTEM OF COUPLED ORDINARY DIFFERENTIAL *
C * EQUATIONS. A GUIDE TO THE REQUIRED INPUT DATA CAN BE *
C * FOUND IN THE SOURCE CODE FOR THE INPUT SUBROUTINE. *
C * OUTPUT IS PRODUCED IN THREE FILES WHICH CONTAIN *
C * THE LOCATION OF THE DEEPEST EDGE OF THE POLLUTANT AT ALL *
C * TIMES, PROFILES AT VARIOUS TIMES AND HISTORIES AT VAR- *
C * IOUS DEPTHS. *
C * *
^
C * WARNING: THE RESULTS FROM THIS CODE HAVE NOT BEEN *
C * VERIFIED BY EITHER LABORATORY OR FIELD STUDIES. THE *
C * SUITIBILITY OF THIS PROGRAM FOR ANY PURPOSE IS NOT IM- *
C * PLIED. THE AUTHOR ACCEPTS NO LIABILITY RESULTING FROM *
C * ITS USAGE. *
\s
C * *
C * AUTHOR JIM WEAVER *
C * THE UNIVERSITY OF TEXAS AT AUSTIN *
C * *
C * *
C * SINGLE PRECISION, ANSI STD. X3.9-1978 FORTRAN 77 *
C * REQUIRED ROUTINES...BLOCK DATA ONE,PFRONT,PCHK,GACHK, *
C * GAHS,RD,CVEL,CONC,OCP,ZMASS, OIL,HIST,SFLUX,INPUT,START,*
C * COUNT, MRKF12,GLQ,INQUAD(INQSET), BINARY, BISEC,WKE (PKE, *
C * SLOPE,KRSET) *
C * OCTOBER 21, 1S87 *
C **************************************************************
£ **************************************************************
DIMENSION NEQ(3),YY(55,3)
- COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /OUTP/ ERCMX,EROMX,NPT,NREG
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /SIMU/ DM,DTPR,TM
EXTERNAL PCHK,PFRONT
DATA ND /55/
NUM(l) = 1
C
CALL INPUT
CALL START (NEQ,YY,ND)
115
-------
c
C ***POLLUTANT AND DISSOLVED CONSTITUENT COMPUTATIONS
CALL MRKF12 (PFPONT,PCHK,TPB,TM,DM,1.E-5,NEQ,NPT,NREG,YY,ND,ET)
C
CALL COUNT
END
BLOCK DATA ONE
COMMON /COUN/ IPR(60) ,NUM(60) ,SEC(60)
COMMON /TRIAL/ DTOLD, DTOLD1, ICC, TEND, ICHAR
COMMON /VAPO/ ABL,AKS,CN,CV(2),CVF(2),DAIR,XJV(2),ZV(2)
COMMON /OUTP/ ERCMX,EROMX,NPT,NREG
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2) ,PZ(2)
COMMON /GAEQ/ CHS(2),DO,EE,FF,GA1,GA2,GA3,GA4,HS,IGA,IP,ZGA
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
DATA NUM /60*0/
DATA DZPMDT, NPT, ISE /O., 2, O/
DATA DPA/DPL,HS,THPZ,TPB,TPE,QP,IP /7*0.,0/
DATA CVF,XJV /4*0./
DATA ERCMX,EROMX /2*0./
DATA DTOLD1,ICC,ICHAR /O.,2*0/
END
SUBROUTINE PFRONT (TT,ZZ,F,ND,NEQ)
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
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
- -*
*
*
*
*
*
*
*
*
*
*
PFRONT POLLUTANT FRONT
PURPOSE THIS SUBROUTINE CONTAINS THE FUNCTIONS F(Y,T)
FOR THE OIL AND CONSTITUENT EQUATIONS THAT ARE USED BY
MRKF12.
THE EQUATIONS ARE DIVIDED INTO 3 GROUPS AS FOLLOWS
GROUP 1. #1 OIL FRONT #2 VOLATILE FRONT
GROUP 2. OIL CHARACTERISTICS
FROM THE SURFACE (#1) TO THE OIL FRONT (UP
TO #51)
GROUP 3. CONSTITUENT CHARACTERISTICS
FROM THE FRONT (#1) TO THE TOP (#11 OR #15)
THIS ARRANGEMENT ALLOWS THE MRKF12 ROUTINE TO VARY THE
NUMBER OF EQUATIONS SOLVED EASILY
INPUT ARGUMENTS
TT = TIME
ZZ(I,J) = DEPTH FOR EQUATION I OF GROUP J
ND = ROW DIMENSION OF ZZ
NEQ(3) = ARRAY OF NUMBER OF EQUATIONS FOR EACH GROUP
OUTPUT ARGUMENTS
F(I,J) = FUNCTION VALUE FOR EQUATION I OF GROUP J
NEQ(3) = NUMBER OF EQUATIONS PER GROUP
-
AUTHOR JIM WEAVER
THE UNIVERSITY OF TEXAS AT AUSTIN
SINGLE PRECISION, ANSI STD. X3. 9-1978 FORTRAN 77
REQUIRED ROUTINES . . . SLOPE, BISEC, GAHS, INQUAD, PKE, OIL, CVEL
RD,CONC
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
it
*
*
*
*
*
*
*
*
*
*
*
*
*
*
116
-------
C * X-21-87
C *
C
DIMENSION F(ND,*),NEQ(*),ZZ(ND,*)
COMMON /CCHR/ CPART,CC(75,2),CO(75),TC(75),ZC(75,2)
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /CONCI/ CLAG/COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2),PZ(2)
COMMON /GAEQ/ CHS(2),DO,EE,FF/GA1,GA2,GA3,GA4,HS,IGA,IP,ZGA
COMMON /MATR/ ETA,OMSOR,OMSWR, SR, SWR
COMMON /OCHR/ DLO,OLAG,SI(75,2),SO(75),TO(75),ZI(75,2)
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /VAPO/ ABL,AKS/CN,CV(2),CVF(2),DAIR,XJV(2),ZV(2)
COMMON /WTER/ IRT,SWMAX,QW,WKS
EXTERNAL GAHS
C
NUM(2) = NUM(2) + 1
PT(2) = TT
ISE = ISE + NEQ(1)+NEQ(2)+NEQ(3)
C
C ***CHARACTERISTICS FOR DEGRADING OILS
DO 10 I=1,NEQ(2)
F(I,2) = 0.
SI (1,2) = SPR
ZI(I,2) = DPA
IF (TT.LT.TO(I)) GO TO 10
XXX = 1.
IF (TT.GT.OLAG) THEN
TAU = AMAX1(OLAG,TO(I))
XXX = EXP(-DLO*(TT-TAU))
END IF
SP = SO(I)*XXX
SI (1,2) = SP
ZI(I,2) = ZZ(I,2)
CALL SLOPE (0,SP,0,SWMAX,SL, O.,0)
F(I,2) = SL/ETA
10 CONTINUE
ZS = DPA + (TT-TPE)*DZPDDT
C
C ***THE OIL FRONT
F(l,l) = 0.
IQF = 1
IF (DLO.EQ.O.) THEN
C ***NON-DEGRADING OILS .
IF (IGA.EQ.l) THEN
DTX = TT - TPB
IF (DTX.GT.DTP) DTX = DTP
IHS = 0
IF (IP.EQ.3) THEN
CHS(2) = HS
IF (TT.GT.TPE) CHS(2) = 0.
ELSE
CHS(2) = QP*DTX - (ZZ(1,1)-DPA)*SPMAX*ETA
IF (CHS(2).LT.O.) CHS(2) = 0.
IF (CHS(l).GT.O..AND.CHS(2).LE.O.) IHS = 1
END IF
END IF
117
-------
IF (ZS.LT.ZZ(1,1),OR.CHS(2).GT.0..OR.IHS.EQ.1) THEN
Z ***CONSTANT MAXIMUM OIL SATURATION
SPF = SPMAX
IF (IGA.EQ.1.AND.CHS(2).GT.O.) THEN
3 ***PRESSURIZED CONDITIONS
IF (IP.EQ.l) THEN
Tl = TT - TPB
IF (Tl.GT.DTP) Tl = DTP
QPF = EE + FF*T1/(ZZ(1,1)-DPA)
ELSE IF (IP.EQ.3) THEN
: ***ANALYTIC SOLUTION FOR FRONT WITH CONSTANT DEPTH PONDING
ISE = ISE - 1
Z2 = ZGA*5.
IF (PZ(l).GE.ZGA) 22 «= PZ(1)*10.
CALL BISEC (GAHS,PZ(1),Z2,0.0001,75,0,TT-TPB,GA1,
1 ZZ(1,1),DT,IE)
IF (ZZ(1,1).GT.DPA) QPF = XKPMX*(l.+CHS<2)/(ZZ (1,1)-DPA))
IQF = 0
END IF
ELSE
: ***GRAVITY DOMINATED FLOW
QPF = XKPMX
END IF
ELSE
: ***GRAVITY DRAINAGE
SL = (ZZ(1,1)-DPA)/(TT-TPE)
CALL INQUAD (WP,50,SPR,1.-SWMAX,SL*ETA,SPF)
CALL PKE (0,SPF,0,SWMAX,QPF,0.,0)
END IF
ELSE IF (DLO.NE.O.) THEN
: ***DEGRADING OILS
CHS(2) = 0.
SL = 0.
SI = SPMAX
LZ = NEQ(2)
DO 20 I=NEQ(2),2,-1
IF (ZI(I,2).GE.ZZ(1,1).AND.ZZ(1,1).GT.ZI(I-1,2)) THEN
LZ = I
GO TO 21
END IF
20 CONTINUE
21 CONTINUE
IF (TT.GT.TO(LZ-l)) SI = SI(LZ-1,2)
IF (ZI(LZ,2).NE.ZI(LZ-1,2))
1 SL =
-------
c
C ***DISSOLVED CONSTITUENTS
C ***CHARACTERISTICS
DO 50 I=1,NEQ(3)
F(I,3) = 0.
ZC(I,2) = DPA
IF (TT.LT.TC(I)) GO TO 50
CALL OIL (ZZ(I,3),TT,SWMAX,NEQ(2),SO)
IF (ABS(TT-TPE).LE.0.0001.AND.I.EQ.NEQ(3)) SO = SPR
CALL CVEL (TT,SO,SWMAX,F(I,3))
XXX = 1.
YYY = 1.
IF (DLC.GT.O..AND.TT.GE.CLAG) XXX = EXP(-DLC*(PT(2)-PT(1)})
IF (DLO.GT.O..AND.SO.GT.O..AND.TT.GT.OLAG) THEN
C ***CORKECTION FOR OIL DECAY
SO2 = SO/EXP( -DLO*(PT(2)-PT(1)) )
YYY = RD(SO2,SWMAX)/RD(SO,SKMAX)
END IF
CC(I,2) = CC(I/1)*XXX*YYY
ZC(I,2) = ZZ(I,3)
50 CONTINUE
F(2,l) = 0.
IF (TT.LT.TPE.OR.XXKV.LE.O..OR.CC(1/2).LT.5.E-5.OR.NEQ(1).LT.2)
1RETURN
C
C ***VOLATILIZATION
CALL OIL (ZZ(2,1),TT,SWMAX,NEQ(2),SO)
IF (ABS(TT-TPE).LT.0.0001.AND.ZZ(2,1).LE.DPA+0.0001) SO = SPR
CALL CONG (TT,ZZ(2,1),SWMAX,SO,NEQ(3),CX,XX,0)
CALL CVEL (TT,SO,SWMAX,F(2,1))
IF (ABS(TT-TPE) .LE.1.E-4.AND.CX.LE.O.) CX = CPART
C ***DIFFUSIVE FLUX
DSOIL = DAIR*( (ETA* (l.-SWMAX-SO) ) **CN) /ETA/ETA
CXV = XXKV*CX
Dl = DSOIL*ABL + DAIR*ZZ(2,1)
IF (Dl.GT.O.) XJV(2) = CXV*DSOIL*DAIR/D1
C ***CUMULATIVE VOLATILE FLUX
CVF(2) = CVF(l) + (TT-PT(1))*0.5*(XJV(2)+XJV(1))*FCV
IF (CX.GT.O.)
1F(2,1) = F(2,l) + XJV(2)/RD(SO,SWMAX)/CX
CV(2) = CX
ZV(2) = ZZ(2,1)
RETURN
END
SUBROUTINE PCHK (TT,ZZ,Z1,ND,NEQ,D1^,DT,DTR,NS,IC,IRR,IR,IRP)
C **************************************************************
C * *
C * PCHK *
C * *
C * *
C * PURPOSE THIS SUBROUTINE CONTAINS A SERIES OF CHECKS *
C * AND CONTROLS FOR THE SOLUTION OF THE KINEMATIC MODEL. *
C * IF A STEP IS NOT REJECTED BY PCHK OR MRKF12, THEN RESULTS*
C * FOR THE CURRENT STEP ARE OUTPUT AND VARIOUS VALUES ARE *
C * ROLLED FOR THE NEXT EVALUATION. *
C * *
C * INPUT ARGUMENTS *
C * TT = SOLUTION TIME JUST OBTAINED * •
119
-------
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
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
**:
ZZ(I,J) - DEPTH FOR EQUATION I GROUP J
Z1(I,J) = PREVIOUS DEPTH FOR EQUATION I GROUP J
ND = ROW DIMENSION OF ZZ()
NEQ(3) = NUMBER OF EQUATION IN EQCH GROUP
DMAX = MAXIMUM TIME STEP
DT = CURRENT TIME STEP
NS = NUMBER OF STEPS
IRP = 1 IF MRKF12 ALREADY HAS REJECTED THE STEP
OUTPUT ARGUMENTS
NEQ(3) = ADJUSTED NUMBER OF EQUATIONS
DMAX = ADJUSTED MAX. TIME STEP
DT = NEW TIME STEP FOR ACCEPTED STEPS
DTR() = SERIES OF PROPOSED REDUCED TIME STEPS
1C = 1 TO END SIMULATION
IRR() = CODES CORESPONDING TO VARIOUS REDUCED STEPS
IR = 1 TO REJECT STEP
AUTHOR JIM WEAVER
THE UNIVERSITY OF TEXAS AT AUSTIN
SINGLE PRECISION, ANSI STD. X3. 9-1978 FORTRAN 77
REQUIRED ROUTINES . . . SFLUX, GACHK, ZMASS, COUNT, HIST
X-21-87
********************************************************
C
C
c
•c
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
DIMENSION DTR(*),IRR(*),NEQ(*),ZZ(ND,*),Z1(ND,*)
COMMON /CCHR/ CPART,CC(75,2) ,CO (75) ,TC(75) ,ZC(75,2)
COMMON /CONCI/ CLAG,COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2),PZ(2)
COMMON /GAEQ/ CHS(2) ,DO,EE,FF,GA1,GA2,GA3,GA4,HS,IGA,IP,ZGA
COMMON /HIPR/ IHIST,NTIMES,PR(10),NZS,HI(5)
COMMON /MATR/ ETA, OMSOR, OMSWR, SR, SWR
COMMON /OCHR/ DLO,OLAG,SI(75,2) ,SO(75),TO(75),ZI (75,2)
COMMON /OUTP/ ERCMX,EROMX/NPT,NREG
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
THPZ/TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /SIMU/ DM,DTPR,TM
COMMON /VAPO/ ABL,AKS,CN,CV(2),CVF(2)/DAIR/XJV(2),ZV(2)
COMMON /WTER/ IRT,SWMAX,QW,WKS
COMMON /TRIAL/ DTOLD, DTOLD1, ICC, TEND, ICHAR
DATA II,12,14,T2 /2,1,1,0./
***SR PCHK PART 1.) VARIOUS RKF CONTROL CHECKS
NUM(3) = NUM(3) + 1
IF (NUM(3).EQ.l) DTOLD = DMAX
***ASSURE AT LEAST 5 STEPS THROUGH OIL APPLICATION
IF (IRR(6).EQ.O.AND.DMAX.GT.0.75*DTP.AND.IAT.NE.2) THEN
DTOLD = DMAX y
DMAX = DTP*.2
DTR(6) = DT
IRR(6) = 1
IR = 1
END IF
120
-------
C ***CHECK FOR POLLUTANT ENDING TIME
IF (IRR(4).EQ.O..AND.IAT.NE.2) THEN
IF (TT.EQ.TPE) THEN
IRR(4) = 1
ELSE IF (TT.GT.TPE) THEN
IR = 1
DTR(4) = TPE - PT(1)
END IF
END IF
C
Tl = TPE
IF (IGA.EQ.1.AND.IP.EQ.1.AND.CHS(2).GT.O.) T1=2.*TT
IF (((COINI.GT.O..AND.I1.LE.NEQ(3)).OR.(DLO.GT.O..AND.II.
1LE.10)).AND.PT(l).LE.T1) THEN
C ***STEP TO THE TIME EACH NEW CHARACTERISTIC STARTS
IF (TC(2).NE.TO(50).AND.IGA.NE.l) STOP 50
IF (IRR(5).EQ.O.) THEN
IF (TC(Il).LT.TT) THEN
C ***n WILL BE INCREMENTED ONLY IF THIS NEW STEP SIZE
DOMINATES
IF (TC(I1).EQ.PT(l)) GO TO 300
IR = 1
DTR(5) = TC(I1)-PT(1)
END IF
300 CONTINUE
END IF
END IF
C
C ***IMPROVE ACCURACY IF CONSTITUENT CHARACTERISTICS LEAVE OIL
305 CONTINUE
IF (COINI.GT.0..AND.IRR(IO).EQ.O.AND.QW.GT.O.) THEN
IF (ZC(I2,2).GE.PZ(2)) THEN
IF (ZC(I2,2).EQ.PZ(2)) THEN
IRR(IO) = 1
GO TO 305
END IF
IR = 1
SL1 = (ZC(I2,2)-ZC(I2,1))/DT
SL2 = (PZ(2)-PZ(1))/DT
TNEW = (PZ(1)-ZC(I2,1)+PT(1)*(SL1-SL2))/(SL1-SL2)
DTR(IO) = TNEW-PT(l)
END IF
END IF
C
C ***MAKE LARGE STEP OUTPUT.LOOK NICE
IF (IRR(ll).EQ.O.AND.PT(2)-PT(1).GE.DTOLD.AND.TT.GT.TPE) THEN
DTR(ll) = FLOAT(IFIX(PT(2)*10000.)/10000) - PT(1)
IF (DTR(ll).GT.O.) IR * 1
END IF
C
IF (DLO.GT.O.) GO TO 310
C
C ***GREEN-AMPT CHECKS IN SR GACHK
IF (IGA.EQ.l) CALL GACHK(TT,ZZ,Z1,ND,NEQ,I1,I2,DMAX,DT,DTR,NS,IC,
1IRR,IR)
C
C ***END OF MAXIMUM OIL SATURATION
IF (IRR(2).EQ.O.AND.((TT.GT.TPE.AND.IGA.EQ.O).OR.
121
-------
1(CHS(1).LE.O.AND.IGA.EQ.l))) THEN
IF (PS(1).LT.SPMAX) IRR(2) = 2
ZS = DPA + (TT-TPE)*DZPDDT
IF (ZZ(1,1).LT.ZS) THEN
SL = DZPMDT
Til = TPB
Zll = DPA
IF (IGA.EQ.l) THEN
SL = (PZ(2)-PZ(1))/(PT(2)-PT(1))
Til = PT(1)
Zll = PZ(1)
ELSE IF (IAT.EQ.2) THEN
Til = TPE
Zll = DPL
END IF
TNEW = (Zll - DPA - T11*SL + TPE*DZPDDT)/(DZPDDT - SL)
DTX = TNEW - PT(1)
IF (DTX.GT.O.) THEN
DTR(2) = DTX
•IR = 1
END IF
END IF
END IF
310 CONTINUE
C
C ***VOLATILIZATION OF ENTIRE MASS OF CONSTITUENT
IF (XXKV.GT.O..AND.IRR(3) .EQ.O) THEN
IF (ZV(2).GT.ZC(1,2)) THEN
IR = 1
DD = PT(2) - PT(1)
SL1 = (ZV(2)-ZV(1))/DD
SL2 = (ZC(1,2)-ZC(1,1))/DD
TNET-7 =
-------
c
C ***RETURN ANY REJECTED STEP
IF (IR.EQ.l.OR.IRP.EQ.l) RETURN
C
C
C
C ***SR PCHK PART IB.) CORRECTIONS FOR CONDITIONS DETECTED ABOVE
C ***ASSURE AT LEAST 5 STEPS THROUGH OIL APPLICATION
IF (IRR(6).EQ.l) THEN
IF (TT.GE.TPE) THEN
DMAX = DTOLD
IRR(6) = 2
END IF
END IF
C
C ***CHECK FOR POLLUTANT ENDING TIME
IF (IRR(4).EQ.l) THEN
IRR(4) = 2
IRR(5) = 1
DT =• 0.0001
DTOLD = DMAX
C IF (DLO.LE.O..AND.IP.NE.1.AND.IGA.NE.1) DMAX = 0.05
IF (XXKV.GT.O..AND.COINI.GT.O.) NEQ(l) = 2
IF (NEQ(l).EQ.O) THEN
NEQ(l) = 1
END IF
END IF
C
C ***STEP TO THE TIME EACH NEW CHARACTERISTIC STARTS
IF (IRR(5).EQ.l) THEN
C ***CORRECT CONCENTRATIONS FOR GREEN AMPT CASES
IF (IGA.EQ.l.AND.QW.GT.O..AND.XXKW.GT.O.) THEN
CALL SFLUX (TC(I1),QO)
CO(II) = QO*COINI/(QO + XXKW*QW)
CC(I1,2) = CO (II)
END IF
II = II + 1
IRR(5) = 0
END IF
C
c ***IMPROVE ACCURACY IF CONSTITUENT CHARACTERISTICS LEAVE OIL
IF (IRR(IO).EQ.l) THEN
IRR(IO) = 0
IF (ABS(ZC(I2,2)-PZ(2)).LT.0.005.OR.DTR(IO).LT.O.) 12=12+1
END IF
C
C ~~***MAKE LARGE STEP OUTPUT LOOK NICE
IF (IRR(ll).EQ.l) THEN
IRR(ll) = 2
ELSE IF (IRR(ll).EQ.2) THEN
IF (PT(2)-PT(1).LT.DTOLD) THEN
IRR(ll) = 0
END IF
END IF
C
C ***GREEN-AMPT CHECKS IN SR GACHK
IF (IGA.EQ.l) CALL GACHK(TT,ZZ,Z1,ND,NEQ,I1,I2,DMAX,DT,DTR,NS,IC,
1IRR,IR)
123
-------
C ***END OF MAXIMUM OIL SATURATION
IF (IRR(2).EQ.l) THEN
IRR(2) = 0
IF (PS(1).LT.SPMAX) IRR(2) = 2
DT = 0.0001
C DMAX = DTOLD
END IF
C
C ***VOLATILIZATION OF ENTIRE MASS OF CONSTITUENT
IF (IRR(3).EQ.l) THEN
IF (ZV(2).GE.ZC(1,2)) THEN
IRR(3) = 2
NEQ(3) = 0
NEQ(l) = 1
ELSE
IRR(3) = 0
END IF
END IF
C
C ***SIMULATION ENDING TIME
IF (IRR(12).EQ.l) 1C = 1
C
C ***PROFILE TIMES
IF (IRR(13).EQ.l) THEN
14 = 14 + 1
IRR(13) = 0
IF (I4.GT.NTIMES) IRR(13) = 2
END IF
C
C ***BIODEGRADATION DURING GREEN AMPT CASES
IF (IRR(8).EQ.l) THEN
IRR(2) = 0
IRR(8) = 3
IF (IP.EQ.l) TPE = TT
DT = 0.0001
IF (ICHAR.EQ.l) THEN
IF (TPE.GT.TC(NEQ(3))) THEN
NEQ(3) = NEQ(3) + 1
ELSE
Nl = NEQ(3)
DO 20 I=N1,1,-1
IF (TPE.GT.TC(I)) THEN
NEQ(3) =1+1
GO TO 21
END IF
20 CONTINUE
STOP 'ICHAR-FAIL'
21 CONTINUE
END IF
TC(NEQ(3)) = TPE
CO(NEQ(3)) = CPART
CC(NEQ(3),2) = CPART
ZC(NEQ(3),2) = DPA
END IF
IF (DO.GT.O.) THEN
C ***TURN ON BIODEGRADATION AFTER PRESSURE FLOW STOPS
DLO = DO
124
-------
NEQ(2) = 51
C ***DEFAULTS FOR CHARACTERISTICS
DZ = (PZ(2)-DPA)*0.1
Z = PZ(2)
DO 50 1=51,41,-!
TO(I) = TPE
ZI(I,1) = Z
ZI(I,2) = Z
ZZ(I,2) = Z
Z = Z - DZ
50 CONTINUE
IF (IP.NE.3) THEN
DO 55 1=1,40
TO(I) = TPE
ZI(I,2) = DPA
ZI(I,1) = DPA
55 CONTINUE
END IF
END IF
END IF
C
C ***SR PCHK PART 2.) DATA OUTPUT FOR ACCEPTED STEPS
VOL = 0.
CVL = 0.
IF (T2+DTPR.LT.PT(2)) THEN
T2 = PT(2)
IZZZ = 0
IF (ABS(TT-2.2461).LE.0.0001) IZZZ = 1
CALL ZMASS (TT,ZZ(1,1),PS(2),NEQ,IZZZ,VOL,CVL)
IF (CVL.LT.1.E-10.AND.VOL.LE.1.E-10) CALL COUNT
ZZ1 = PZ(2)
ZZ2 = ZC(1,2)
ZZ3 = ZC(NEQ(3),2)
CCC = CC(1,2)*XXKW
IF (VOL.LE.l.E-10) ZZ1 = 0.
IF (XXKV.GT.O.AND.ZV(2).GT.ZZ3) ZZ3 = ZV(2)
IF (CVL.LT.l.E-10) THEN
ZZ2 = 0.
ZZ3 = 0.
CCC = 0.
END IF
WRITE (6,9300) NS,TT,ZZ1,CHS(2),PS(2),PF(2),VOL,ZZ2,ZZ3,CCC,CVL,
1 CVF(2)
END IF
C ***HISTORIES
IF (IHIST.NE.O) CALL HIST (ZZ,Z1,ND,NEQ,TT)
C ***PROFILES
DO 25 I=1,NTIMES
Tl = PT(1)
IF (Tl.LT.PR(I).AND.PR(I).LE.TT) THEN
DA = PT(2) - PT(1)
ZP = Z1(1,1) + (PR(D-T1)MZZ(1,1)-Z1(1,1))/DA
SP = PS(1) + (PR(D-T1)*(PS(2)-PS(1))/DA
CALL ZMASS (PR(I),ZP/SP,NEQ,1,VOL,CVL)
END IF
25 CONTINUE
C
C
125
-------
***SR PCHK PART 3.) ROLL VALUES FOR NEXT EVALUATION
PT(1) = PT(2)
PS(1) = PS (2)
PZ(1) = PZ(2)
PF(1) = PF(2)
DO 2 I=1,NEQ(2)
81(1,1) = SI(1,2)
CONTINUE
DO 4 I=1,NEQ(3)
L) = CC(I,2)
L) = ZC(I,2)
4 CONTINUE
CHS(l) = CHS(2)
IF (XXKV.GT.O..AND.TT.GT.TPE) THEN
ZV(1) = ZV(2)
CV(1) = CV(2)
XJV(l) = XJV(2)
CVF(l) = CVF(2)
END IF
IF (DLO.EQ.O.) RETURN
IF (ZZ(NEQ(2)-1,2).GT.ZZ(1,1)) NEQ(2) = NEQ(2) - 1
9300 FORMAT (1X,I5,11F10.4)
RETURN
END
SUBROUTINE GACHK (TT,ZZ,Z1,ND,NEQ, II, I2,DTMAX,DT,DTR,NS, 1C, IRR, IR)
DIMENSION DTR(*),IRR(*),NEQ(*),ZZ(ND,*),Z1(ND,*)
COMMON /CCHR/ CPART,CC(75,2) ,CO(75) ,TC(75) ,ZC(75,2)
COMMON /CONCI/ CLAG,COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2) ,PZ(2)
COMMON /GAEQ/ CHS (2),DO,EE,FF,GA1,GA2,GA3,GA4,HS,IGA,IP,ZGA
COMMON /MATR/ ETA, OMSOR, OMSWR, SR, SWR
COMMON /OCHR/ DLO,OLAG,SI(75,2) ,SO(75) ,TO(75),ZI(75,2)
COMMON /TRIAL/ DTOLD, DTOLD1, ICC, TEND, ICHAR
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT, IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
: ***GREEN-AMPT CHECKS
NUM(4) = NUM(4) + 1
IF (TT.LT.TEND) IRR(8) = 0
IF (IRR(8).EQ.3.0R.PT(2).LT.TPE) RETURN
IF (IP.EQ.3.AND.PT(1).GE.TPE) IRR(8) =1
: ***EXTRA CHARACTERISTICS FOR THE CONSTITUENT IN GREEN-AMPT FLOW
IF (IP.EQ.1.AND.CHS(2).GT.O..AND.COINI.GT.O.) THEN
ICHAR = 1
Nl = NEQ(3)
- DO 10 I=N1,50
IF (TC(NEQ(3)).LT.TT) THEN
NEQ(3) = NEQ(3) + 1
TC(NEQ(3)) = TC(NEQ(3)-1)+DTP*0.1
CO(NEQ(3)) = CPART
CC(NEQ(3),2) = CPART
ZC(NEQ(3),2) = DPA
Z1(NEQ(3),3) = DPA
IR *= 1
DTR(9) = DT
ELSE
GO TO 15
126
-------
END IF
10 CONTINUE
C ***ELIMINATE EVEN NUMBERED CHARACTERISTICS
IE = 2
DO 12 1=3,51,2
TC(IE) = TC(I)
CO(IE) = CO(I)
ZC(IE,2) = ZC(I,2)
ZC(IE,1) = ZC(I,1)
CC(IE,2) = CC(I,2)
CC(IE,1) = CC(I,1)
ZZ(IE,3) = ZZ(I,3)
Z1(IE,3) = Z1(I,3)
IE = IE + 1
12 CONTINUE
NEQ(3) = IE - 1
C ***CORRECT INTEGERS FOR PCHK CHECKS
II = NEQ(3) - 1
IF (MOD(12,2).EQ.O.) 12 = 12 + 1
IS = 1
NN = 1
2 CONTINUE
IF (IS.EQ.I2) GO TO 3
NN = NN + 1
IS = IS + 2
GO TO 2
3 CONTINUE
12 = 12 - NN + 1
15 CONTINUE
END IF
IF (IP.EQ.1.AND.IRR(7).EQ.O.AND.TT.GE.TPE) THEN
C ***APPROXIMATE TIME OF ZERO SURFACE PONDING
IF (CHS(2).LE.O.AND.ICC.EQ.O) THEN
IR = 1
DTOLD1 = DTMAX
DTR(7) = 0.0001
END IF
ELSE IF (IRR(7).EQ.l) THEN
IRR(2) = 5
IRR(7) = 2
ICC = 1
DTMAX = 0.05
END IF
IF (IRR(8) .EQ.O.AND.ICC.EQ.LAND.CHS(2).LE.O.) THEN
C ***TIME OF ZERO SURFACE.PONDING
IF (ZZ(1,1).GE.ZGA) THEN
SL = (PZ(2) - PZ(1))/(PT(2) -
IR = 1
DTR(8) = (ZGA - PZ(1))/SL
TEND = PT(1) + DTR(8)
RETURN
END IF
ELSE IF (IRR(8).EQ.l) THEN
IF (DTOLD1.GT.O.) DTMAX = DTOLD1
END IF
RETURN
END
SUBROUTINE GAHS (IS,ZZ,NN,XY,DT,XX,IE)
127
-------
COMMON /COUN/ I?R(60),NUM(60),SEC (60)
COMMON /GAEQ/ CHS(2),DO,EE,FF,GA1,GA2,GA3,GA4,HS,IGA,IP,ZGA
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
NUM(5) = NUM(5) + 1
C ***ANALYTIC SOLUTION OF THE GREEN AMPT MODEL WITH CONSTANT PONDING
DEPTH
GAS = ALOG(ZZ+GA3)
DT = GA1*(ZZ - GA3*GA5 - GA4 - DPA*(GA5-GA2))
RETURN
END
FUNCTION RD(SO,SW)
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /CONCI/ CLAG,COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /MATR/ ETA, OMSOR, OMSWR, SR, SWR
C ***RETARDATION FACTOR
NUM(6) = NUM(6) + 1
RD = ETA*(SW*XXKW+SO+XSF1)
RETURN
END
SUBROUTINE CVEL (TT,SO,SW,DZCDT)
COMMON /COUN/ IPR(60),NUM(60) ,SEC(60)
COMMON /CONCI/ CLAG,COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2),PZ(2)
COMMON /GAEQ/ CHS (2) ,DO,EE,FF,GA1,GA2,GA3,GA4,HS, IGA, IP,ZGA
COMMON /MATR/ ETA, OMSOR, OMSWR, SR, SWR
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /WTER/ IRT,SWMAX, QW,WKS
C ***CONCENTRATION VELOCITY
NUM(7) = NUM(7) + 1
QWC = QW
IF (IGA.EC.LAND.CHS (2) .GT.O.) THEN
C ***PRESSURE FLOW
QOC = 0.
IF (SO.GT.SPR) THEN
QOC = PF(2)
IF (TT.LT.PT(2)) THEN
IF (TT.LT.PT(l)) STOP 'CVEL-TIME1
SL = (PF(2)-PF(1))/(PT(2)-PT(1))
QOC = PF(1) + (TT-PT(1))*SL
END IF
END IF
ELSE
C ***KENEMATIC FLOW
CALL PKE (0,SO,0,SW,QOC,0.,0)
END IF
DZCDT = (QWC*XXKW + QOC)/RD (SO,SW)
RETURN
END
SUBROUTINE CONC (TT,ZZ,SW,SO,NEQ3,CO,CMSS,IV)
C * *
C * CONC *
C * *
C * CONC DETERMINES THE OIL PHASE CONCENTRATION AT TIME TT *
C * AND DEPTH ZZ. OTHER NEEDED INFORMATION ARE THE WATER *
C * AND OIL SATURATIONS (SW,SO) AND THE NUMBER OF CONC. *
128
-------
C * CHARACTERISTICS (NEQ2) . CMSS IS THE MASS IN ALL PHASES *
C * IF IV=1 VOLATILIZATION IS CONSIDERED. *
C * *
REAL CT(2)
COMMON /CCHR/ CPART,CC(75,2),CO(75),TC(75),ZC(75,2)
COMMON /CONCI/ CLAG,COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /GAEQ/ CHS (2) ,DO,EE,FF,GA1,GA2,GA3,GA4, HS, IGA, IP, ZGA
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2),PZ(2)
COMMON /MATR/ ETA,OMSOR,OMSWR,SR,SWR
COMMON /OCHR/ DLO,OLAG,SI(75,2) ,SO(75) ,TO(75) ,ZI(75, 2)
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZrTPB,TPE,TPR,WP(50)/XKPMX,XK22
COMMON /VAPO/ ABL,AKS,CN,CV(2),CVF(2),DAIR,XJV(2),ZV(2)
COMMON /WTER/ IRT, SWMAX, QW, WKS
NUM(8) = NUM(8) + 1
C
CO = 0.
CMSS = 0.
IF (COINI.LE.O.) RETURN
DT = PT(2) - PT(1)
CT(1) = 0.
CT(2) = 0.
DO 100 IT = 2,1,-1
C ***FIND NUMBER OF CHARACTERISTICS
IF (TT.LT.TPE*1.15) THEN
DO 20 I=2,NEQ3
IF (ZC(I,IT).EQ.DPA) THEN
LZ = I
IF (TT.LT.TPE) LZ = LZ - 1
GO TO 21
END IF
20 CONTINUE
21 CONTINUE
ELSE
LZ = NEQ3
END IF
C ***FIND DEPTH OF DEEPEST & SHALLOWEST CHARACTERISTICS
ZLOW = ZC(LZ,IT)
CLOW = CC(LZ,IT)
ZHIGH = ZC(1,IT)
C ***DETERMINE THE CONCENTRATION
IF (ZZ.GT.ZHIGH) GO TO 100
IF (ZZ.LT.ZLOW) THEN
IF (TT.GT.TPE) THEN
GO TO 100
C "~ ELSE IF (ZVF.GT.DPA.AND.ZVF.LE.ZZ.AND.ZZ.LT.ZLOW) THEN
C ***DURING VOLATILIZATION
C SL = (CLOW-CVF)/(ZLOW-ZVF)
C CO = CVF + (ZZ-ZVF)*SL
ELSE IF (TT.LE.TPE) THEN
c ***DURING LOADING
CC1 = CPART
IF (IGA.EQ.l) THEN
CALL SFLUX (TT,QO)
IF (QO.GT.O..AND.QW.GT.O.) CC1 = QO*COINI/(QO + XXKW*QW)
END IF
129
-------
CT(IT) = CC1 + (ZZ-DPA)*(CLOW - CC1)/(ZLOW-DPA)
ELSE
STOP 'CONC-ZLOW
END IF
IF (IZX.EQ.l) WRITE (7,1510) CT (IT), CLOW, CPART, ZZ,ZLOW
1510 FORMAT (IX,'CONG CT(IT),CLOW,CPART,ZZ,ZLOW,5F10.3)
ELSE
: ***LINEAR INTERPOLATION
CALL BINARY (ZZ,ZC(1,IT),LZ,1,IN2)
Cl = CC(IN2,IT)
C2 = CC(IN2+1,IT)
Zl = ZC(IN2,IT)
Z2 = ZC(IN2+1,IT)
IF (ZC(IN2+1,IT) .LT.PZ(2) .AND.PZ(2) .LT.ZC(IN2, IT) .
1 AND.DLO.GT.O.) THEN
3 ***APPROXIMATE CONCENTRATION NEAR THE OIL FRONT
IF (ZZ.LT.PZ(2)) THEN
Cl = CC(IN2+1,IT)
Zl = PZ(2)
•ELSE
Z2 = PZ(2)
END IF
END IF
SL = (C2-C1)/(Z2-Z1)
CT(IT) = Cl + (ZZ-Z1)*SL
END IF
IF (TT.EQ.PT(2)) THEN
IF (IV.EQ.1.AND.ZZ.LT.ZV(2)) RETURN
CO = CT(2)
GO TO 200
END IF
100 CONTINUE
DD = TT-PT(l)
IF (TT.LT.PT(1)) STOP 'CONC-TIME'
: ***BILINEAR INTERPOLATION ON A QUADRALATERIAL
IF(IZX.EQ. 1) WRITE (7,4330)TT, ZZ, ZC(NEQ3,2),ZC(1,1)
4330 FORMAT(IX,'CONC TT,ZZ,ZC(NEQ3,2) ,ZC(1,1)'.4F10.4)
IF (TT.GT.TPE.AND. (ZZ.LT.ZC(NEQ3,2) .OR.ZZ.LT.ZV(2))) THEN
IF (ZZ.LT.ZV(l).OR.ZZ.LT.ZC(NEQ3,1)) RETURN
Cl = CT(1)
Tl = PT(1)
IF (XXKV.GT.O..AND.IV.EQ.1) THEN
: ***VOLATILE FRONT (TOP)
ZVF = ZV(1) + DD*(ZV(2)-ZV(1))/DT
IF (ZZ.LT.ZVF) RETURN-
C2 = CV(1) + (ZZ-ZV(1))*(CV(2)-CV(1))/(ZV(2)-ZV(1))
T2 = PT(1) + (ZZ-ZV(1))*DT/(ZV(2)-ZV(1))
ELSE
: ***NONVOLATILE FRONT (TOP)
ZF = ZC(NEQ3,1)+DD*(ZC(NEQ3,2)-ZC(NEQ3,1))/DT
IF (ZZ.LT.ZF) RETURN
C2 = CC(NEQ3,1)+DD*(CC(NEQ3,2)-CC(NEQ3,1))/DT
T2 = PT(1) + (ZZ-ZC(NEQ3,1))*DT/(ZC(NEQ3,2)-ZC(NEQ3,1))
END IF
IF (ZZ.GT.ZC(1,1)) THEM
; ***VERY HIGHLY SKEWED QUADS.
ZF = ZC(1,1) + DD*(ZC(1,2)-ZC(1,1))/DT
IF (ZZ.GT.ZF) RETURN
130
-------
Cl = CC(1,1) + DD*(CC(1,2)-CC(1,1))/DT
Tl = PT(1) + (ZZ-ZC(1,1))/(ZC(1,2)-ZC(1,1)
END IF
CO = Cl + (TT-T1)*(C2-C1)/(T2-T1)
ELSE IF (ZZ.GT.ZCd,!)) THEN
C ***LEADING EDGE OF CONSTITUENT
ZF = ZC(1,1) + DD*(ZC(1,2)-ZC(1,1))/DT
IF (ZZ.GT.ZFfO.000001) RETURN
Tl = PT(1) + (ZZ-ZC(1,1))*DT/(ZC(1,2)-ZC(1,1))
Cl = CC(1,1) + (ZZ-ZC(1,1))*(CC(1,2)-CC(1,1))/
CO = Cl
IF (ABS(TT-Tl) .GT.0.00001) CO = Cl + DD* (CT(2)-C1) / (TT-T1)
ELSE
C ***IN CENTER OF DISTRIBUTION
CO = CT{1) + DD*(CT(2)-CT(1))/DT
END IF
200 CONTINUE
CMSS = RD(SO,SW)*CO
RETURN
END
SUBROUTINE OCP (TT,NEQ,A,B, IW, SUM!, SUM2)
DIMENSION NEQ(*),CS(7),CY(7)
COMMON /CONCI/ CLAG,COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT, IAT,PKS,PVOL,QP,SPMAX, SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /WTER/ IRT, SWMAX, QW, WKS
DATA CY /-.949108,-.741531,-.405845,0.,.405845,.741531,.949108/
DATA CS /.129485,.279705, .381830, .417959, .381830,.279705,.129485/
C ***MODIFIED GAUSS-LAGRANGE QUADRATURE FOR OIL AND CONSTITUENT
C ***II-17-87, METHOD FROM ABRAMOWITZ AND STEGUN P. 916
NUM(9) = NUM(9) + 1
XC1 = (B - A)*0.5
XC2 = (B + A)*0.5
C
C ***FUNCTION EVALUATIONS AT GAUSS POINTS
SUM! = 0.
SUM2 = 0.
DO 10 1=1,7
Y = CY(I)*XC1 + XC2
CALL OIL (Y,TT,SWMAX,NEQ(2),SO)
CALL CONC (TT,Y,SWMAX,SO,NEQ(3),CO,CM,1)
IF (SO.LE.O.AND.CO.LE.O.) THEN
C IF (IW.EQ.l) WRITE (9,1340) TT,Y,SO,CO*XXKW
C1340 FORMAT (IX,'OCP TT,Y,SO,CO',4F10.4)
IF (IW.EQ.l) WRITE (9,9100)
9100 FORMAT (/)
RETURN
END IF
IF (IW.EQ.l) WRITE (9,9101) Y,SO,CO*XXKW
SUM1 = SUM1 + SO*CS(I)
SUM2 = SUM2 + CM*CS(I)
10 CONTINUE
SUM! = SUM1*XC1
SUM2 = SUM2*XC1
C
9101 FORMAT (15X,F10.4,5X,3F10.4)
RETURN
131
-------
END
SUBROUTINE ZMASS (TH, ZH, SH,NEQ, IW, VW,CVL)
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
**•*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
**1
**********************************************************
ZMASS
PURPOSE THIS SUBROUTINE CREATES SATURATION AND CON-
CENTRATION PROFILES AT GIVEN TIMES. THE VOLUME/AREA
OR MASS IN THE PROFILE IS OBTAINED BY GAUSS -LAGRANGE
QUADRATURE OF THE PROFILE (ROUTINE GLQ OR OCP) .
INPUT ARGUMENTS
TH = TIME
ZH = MAXIMUM DEPTH OF OIL
SH = OIL SATURATION AT ZH
NEQ = ARRAY CONTAINING NUMBER OF EQUATIONS
IW = WRITE RESULTS INTO FILE PROF IF IW=1
OUTPUT ARGUMENTS
WV = VOLUME/AREA OF OIL
CVL = MASS/AREA OF CONSTITUENT
AUTHOR JIM WEAVER
THE UNIVERSITY OF TEXAS AT AUSTIN
SINGLE PRECISION, ANSI STD. X3. 9-1978 FORTRAN 77
REQUIRED ROUTINES . . . CONC, OIL, OCP, GLQ
X-21-87, VII-18-88
**************************************** *****************^
*•**
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*•**
DIMENSION NEQ(*),ZX(10)
CHARACTER NT*15
COMMON /CCHR/ CPART,CC(75,2) ,CO(75) ,TC{75),ZC(75,2)
COMMON /CONCI/ CLAG,COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2),PZ(2)
COMMON /GAEQ/ CHS(2),DO,EE,FF,GA1,GA2,GA3,GA4,HS,IGA,IP,ZGA
COMMON /MATR/ ETA, OMSOR, OMSWR, SR, SWR
COMMON /NAME/ NT(15)
COMMON /OCHR/ DLO,OLAG,SI(75, 2) ,SO(75),TO(75),ZI(75,2)
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /VAPO/ ABL,AKS/CN,CV(2),CVF(2),DAIR,XJV(2),ZV(2)
COMMON /WTER/ IKT,SWMAX,CW,WKS
NUM(IO) = NUM(IO) + 1
WV = 0.
NN = 1
DO 5 11=1,10
ZX(II) = 0.
CONTINUE
IF (TH.GT.TPE) THEN
SS = SPR
IF (DLO.GT.O..AND.TH.GT.OLAG)
1 SS = SPR*EXP( - DLO*(TH-OLAG))
ELSE
132
-------
SS = SPMAX
END IF
C ***DIVIDE THE PROFILE INTO SUBSECTIONS FOR INTEGRATION
C ***OIL LOCATIONS
NN = NN + 1
ZX(NN) = ZH
IF (SH.GE.SPMAX.AND.DLO.LE.O.) THEN
Tl = TPE
IF (CHS(2).GT.O..AND.IGA.EQ.l) Tl = 2.*TH
IF (TH.GT.T1) THEN
NN = NN + 1 .
ZD = DPA + (TH-TPE) *DZPDDT
ZX(NN) = ZD
NN = NN + 1
ZX(NN) = DPA + 0.3333*(ZD-DPA)
NN = NN + 1
ZX(NN) = DPA + 0.6666*(ZD-DPA)
END IF
ELSE IF (DLO.GT.O.) THEN
IF (TH.GT.TPE.AND.NEQ(2).GT.41.AND.ZI(41,2).GT.DPA
1 .AND.ZH.GT.ZI(41,2))THEN
NN = NN + 1
ZX(NN) = ZI(41,2)
NN = NN + 1
ZX(NN) = DPA + (ZI(41,2)-DPA)*0.5
END IF
END IF
C
IF (COINI.GT.O.) THEN
C ***CONSTITUENT PROFILE LOCATIONS
NN = NN + 1
ZX(NN) = ZC(1,2)
ZCT = ZC(NEQ(3),2)
ZCV = ZV(2)
IF (ZCV.GT.ZCT) ZCT = ZCV
NN = NN -f 1
ZX(NN) = ZCT
END IF
C
C ***ELIMINATE ZEROS
1 = 1
NX = NN
20 CONTINUE
1 = 1 + 1
IF (I.EQ.51) STOP 'ZMASS ELIMINATE ZEROS'
IF (ZX(I).LE.O.) THEN
NX = NX - 1
DO 15 J=I,NX
ZX(J) = ZX(J+1)
15 CONTINUE
1 = 1-1
END IF
IF (I.LT.NX) GO TO 20
C
C ***SORT THE DEPTHS
DO 30 J=l,100
IH = 0
133
-------
DO 25 I=1,NX-1
IF (ZX(I).GT.ZX(I+1)) THEN
IH = 1
ZTEMP = ZX(I)
ZX(I) = ZX(I+1)
ZX(I+1) = ZTEMP
END IF
25 CONTINUE - • .
IF (IH.EQ.O) GO TO 31
30 CONTINUE
STOP 'ZMASS SORTING FAILURE1
31 CONTINUE
NN = NX
C
CVL = 0.
CALL CONC (TH,DPA,SWMAX,SS,NEQ(3),CCC,CXX,1)
IF (IW.EQ.l) THEN
WRITE (9,9100) TH,(NT(II),11=1,15)
WRITE (9,9106)
WRITE (9,9107) DPA,SS,CCC*XXKW
END IF
C
C ***OIL AND CONSTITUENT PROFILE
DO 10 I=1,NN-1
OVOL = 0.
CVOL = 0.
IF (IW.EQ.l) THEN
CALL OIL (ZX(I)*1.000001,TH,SWMAX,NEQ(2),SSS)
CALL CONC (TH,ZX(I)*1.000001,SWMAX,SSS,NEQ(3),CCC,CXX,1)
WRITE (9,9107) ZX(I),SSS,CCC*XXKW
END IF
CALL OCP (TH,NEQ,ZX(I),ZX(I+1),IW,OVOL,CVOL)
IF (IW.EQ.l) THEN
CALL OIL (ZX(I+1)*.999999,TH,SWMAX,NEQ(2),SSS)
CALL CONC (TH,ZX(I+1)*.999999,SWMAX,SSS,NEQ(3),CCC,CXX,1)
WRITE (9,9107) ZX(I+1),SSS,CCC*XXKW
END IF
WV = VW + OVOL*ETA
CVL = CVL + CVOL*FCV
10 CONTINUE
IF (IW.EQ.l) WRITE (9,9115) WV,CVL
C
RETURN
9100 FORMAT (1H1/10X,'SATURATION AND CONCENTRATION PROFILE AT ',F10.4/
2 10X,50(1H*)/
3 3(10X,5A10/)/)
9106 FORMAT (15X, 3X, ' DEPTH', 2X, 5X, 4X, ' SAT.', 2X, ' CONC. (WATER) • /
2 15X,35(1H=))
9107 FORMAT (15X,F10.4,5X,2F10.4, '*•)
9110 FORMAT (15X,F10.4,5X,F10.4//
1 10X,'VOLUME/AREA ',F10.4///)
9115 FORMAT (//10X,'OIL VOLUME/AREA ',F10.4, ' CONSTITUENT MASS/AREA ',
1 F10.4)
END
SUBROUTINE OIL (ZZ,TT,SW,NEQ2,SO)
C * *
C * OIL *
134
-------
c * *
C * FOR TIME TT, (BETWEEN PT(1) AND PT(2)), OIL FINDS THE OIL *
C * SATURATION AT DEPTH ZZ GIVEN THE WATER SATURATION SW *
C * *
o
REAL SS(2)
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2),PZ(2)
COMMON /GAEQ/ CHS(2),DO,EE,FF,GA1,GA2,GA3,GA4,HS,IGA,IP,ZGA
COMMON /MATR/ ETA, OMSOR, OMSWR, SR, SWR
COMMON /OCHR/ DLO,OLAG,SI(75,2) ,SO(75) ,TO(75) ,ZI(75,2)
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /WTER/ IRT,SWMAX,QW,WKS
NUM(ll) = NUM(ll) + 1
SO = 0.
C ***IMMOBILIZED OIL LAYERS
IF (IAT.EQ.2.AND.QP.LE.O.) THEN
IF (DPA.LE.ZZ.AND.ZZ.LE.DPL) THEN
XXX = 1.
IF (TT.GT.OLAG) XXX = EXP(-DLO*(TT-OLAG))
SO = SPMAX*XXX
END IF
RETURN
END IF
IF (ZZ.LT.DPA) RETURN
DT = PT(2) - PT(1)
ZF = PZ(2)
IF (TT.LT.PT(2).AND.PT(2).NE.PT(l)) THEN
IF (TT.LT.PT(l)) STOP 'BAD-OIL1
ZF = PZ(1) + (TT-PT(1))*(PZ(2)-PZ(1))/DT
END IF
IF (ZZ.GT.ZF) RETURN
C
C ***MOBILE OILS
IF (DLO.EQ.O.) THEN
Z22 = DPA + (TT - TPE)*DZPDDT
IHS = 0
IF (IGA.EQ.1.AND.IP.EQ.1.AND.CHS(2).GT.O.) Z22 = DPA
IF (ZZ.GE.Z22.AND.((ZZ.LT.ZF.OR.ABS(ZZ-ZF).LE.l.E-5)
1 .OR. (IP.EQ.LAND.CHS(2) .GT.O. .AND.ZZ.LE.ZF))) THEN
SO = SPMAX
ELSE IF (ZZ.LE.ZF) THEN
SL = (ZZ-DPA)/(TT-TPE)
CALL INQUAD (WP,50,SPR,1.-SWMAX,SL*ETA,SO)
END IF
ELSE
C ***DEGRADING OILS
DO 100 IT=2,1,-1
ZF = PZ(IT)
I-F (PT(l).LE.TPE) THEN
LZ = NEQ2
DO 20 I=NEQ2-1,1,-1
IF (ZI(I,IT).LE.DPA) THEN
LZ = I + 1
GO TO 21
END IF
20 CONTINUE
135
-------
21 CONTINUE
ELSE
LZ = 1
END IF
IF (ZZ.LE.ZI(LZ,IT) .AND.TT.LE.TPE) THEN
SL = 0.
IF (ZI(LZ,IT) .GT.DPA) SL = (SI (LZ,IT) -SPMAX) /
1 (ZI(LZ,IT)-DPA)
SS(IT) = SPMAX + SL*(ZZ-DPA)
ELSE
IX = 0
CALL BINARY (ZZ, ZI (LZ-1,IT) ,NEQ2-LZ+1,1, IX)
IF (LZ.GT.l) IX = LZ + IX - 1
SL = (SI (IX+1, IT) -SI (IX, IT) ) / (ZI (IX+1, IT) -ZI (IX, IT) )
SS(IT) = SI (IX, IT) + (ZZ-ZI (IX, IT) ) *SL
END IF
IF (TT.EQ.PT(IT)) THEN
IF (DPA.LE.ZZ.AND.ZZ.LE.ZF) SO = SS(IT)
RETURN
END IF
100 CONTINUE
SL = (SS(2)-SS(1))/DT
SO = SS(1) + (TT-PT(1))*SL
END IF
RETURN
END
SUBROUTINE HIST (ZZ,Z1,ND,NEQ,TT)
DIMENSION NEQ (*) , SC (15) , TZ (10) , ZZ (ND, * ) , ZI (ND, * )
COMMON /CCHR/ CPART,CC(75,2) ,CO (75) ,TC(75) , ZC(75,2)
COMMON /COUN/ IPR(60) ,NUM(60) , SEC (60)
COMMON /CONCI/ CLAG,COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /FRNT/ ISE,PF (2) ,PS (2) ,PT(2) ,PZ (2)
COMMON /GAEQ/ CHS (2) ,PO,EE,FF,GA1,GA2,GA3,GA4,HS,IGA,IP, ZGA
COMMON /HIPR/ IHIST,NTIMES,PR(10),NZS,HI(5)
COMMON /POLL/ DPA, DPL, DTP, DZPMDT,DZPDDT,IAT,PKS,PVOL,QP, SPMAX, SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /WTER/ IRT,SWMAX,QW,WKS
C
C ***PRODUCES HISTORIES IN FILE HIST (10)
NUM(12) = NUM(12) + 1
II = 1
TZ(1) = TT
II = 1
DO 10 1=1,15
SC(I) = 0.
1C CONTINUE
C
C ***CHECK FOR OIL OR CONSTITUENT FIRST REACHING THE OUTPUT DEPTHS
C ***REASSIGN OUTPUT TIMES IF SO
IF (COINI.GT.O.) THEN
DO 15 1=1, NZS
IF (ZC(1,1).LE.HI(I).AND.HI(I).LT.ZC(1,2)) THEN
II = II + 1
SL = (PT(2)-PT(1))/(ZC(1,2)-ZC(1,1))
II = II
END IF
136
-------
C ***CONSTITU£NT PASSING OUTPUT DEPTHS
IF (ZC(NEQ(3),1).LE.HI(I).AND.HI(I).LT.ZC(NEQ(3),2)) THEN
II = II + 1
SL = (PT(2)-PT(1))/(ZC(NEQ(3),2)-ZC(NEQ(3),1))
TZ(I1) = PT(1) + (HI(I)-ZC(NEQ(3),1))*SL
II = II
END IF
15 CONTINUE
END IF
C WRITE (10,1305) TZ(1),TT
C
C ***OIL FIRST REACHING DEPTHS
DO 16 1=1,NZS
IF (PZ(1).LE.HI(I).AND.HI(I).LT.PZ(2)) THEN
'II =11+1
SL = (PT(2)-PT(1))/(PZ(2)-PZ(1))
TZ(I1) = PT(1) + (HI(I)-PZ(1))*SL
II = II
END IF
16 CONTINUE
C
C ***SORT THE OUTPUT TIMES
C WRITE (10,1305) TZ(1),TT
DO 18 L=1,II
J = 0
DO 17 14=1,11-1
IF (TZ(M+l).LT.TZ(M)) THEN
TEMP = TZ(M)
TZ(M) = TZ(M+1)
TZ(M+1) = TEMP
J = 1
END IF
17 CONTINUE
IF (J.EQ.O) GO TO 19
18 CONTINUE
19 CONTINUE
C
C ***COLLECT INFORMATION FOR EACH OUTPUT TIME
C WRITE (10,1305) TZ(1),TT
DO 30 L=1,II
DO 20 1=1,NZS
J = 3*1-2
CALL OIL (HI(I),TZ(L),SWMAX,NEQ(2),SC(J))
SC(J+1) = 0.
IF (IGA.EQ.1..AND.SC(J).GE.SPMAX) THEN
SC(J+1) = PF(2)
— IF (TZ(L).LT.PT(2)) THEN
IF (TZ(L).LT.PT(l)) STOP 'HIST-TIME'
SL = (PF(2)-PF(1))/(PT(2)-PT(1))
SC(J+1) = PF(1) + (TZ(L)-PT(1))*SL
END IF
ELSE IF (SC(J).GT.O.) THEN
CALL PKE (0,SC(J),0,SWMAX,SC(J-KL),0.,0)
END IF
SC(J+2) = 0.
IF (COINI.GT.O.) CALL CONC (TZ(L),HI(I) ,SWMAX,SC(J),NEQ(3),
1 SC(J+2),X,1)
SC(J+2) = SC(J+2)*XXKW
137
-------
20 CONTINUE
IX = 0
DO 25 IJ = 1,NZS*3
IF (SC(IJ).GT.O.) IX = 1
25 CONTINUE
IF (IX.EQ.l) WRITE (10,9201) TZ(L), (SC(J), J=l,3*NZS)
30 CONTINUE
9201 FORMAT (1X,F9.4,3(3F10.4))
RETURN
END
SUBROUTINE SFLUX (TT,QO)
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2),PZ(2)
COMMON /GAEQ/ CHS(2),DO,EE,FF,GA1,GA2,GA3,GA4,HS,IGA, IP,ZGA
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
C ***RETURNS THE FLUX AT THE SURFACE TO CORRECT SOME CONCENTRATIONS
NUM(13) = NUM(13) + 1
QO = 0.
IF (IGA.EQ.O) THEN
IF (TT.LE.TPE) QO = QP
ELSE
IF (CHS(l).GT.O..OR.(TT.LE.TPE)) THEN
QO = PF(2)
IF (PT(1).LE.TT.AND.TT.LT.PT(2)) THEN
SL = (PF(2)-PF(1))/(PT(2)-PT(1))
QO = PF(1) + (TT-PT(1))*SL
END IF
END IF
END IF
RETURN
END
SUBROUTINE INPUT
C *
C * REQUIRED INPUT DATA
C * (SEE EXPLANATORY NOTES BELOW)
C *
C * CARD 1-3...RUN TITLE...(5A10/5A10/5A10)
C * NT (15) RUN TITLE 3 LINES OF 50 CHARACTERS EACH
C *
C * CARD 4 MATRIX PROPERTIES... (5F10.4)
C * WKS SATURATED HYDRAULIC CONDUCTIVITY (WATER)
C * XLAMB BROOKS AND COREY'S LAMBDA
C * ETA POROSITY
C * RHOS BULK DENSITY OF MATRIX
C * SWR RESIDUAL WATER SATURATION
C *
C * CARDS WATER PROPERTIES... (2F10.4,110, F10.4)
C * WMU DYNAMIC VISCOSITY OF WATER
C * WRHO DENSITY OF WATER
C * IRT RAINFALL INPUT TYPE: 1=FLUX SPECIFIED
C * 2=SATURATION SPECIFIED
C * 3=CONSTANT HEAD
C * QW/SWMAX CONSTANT WATER FLUX OR SATURATION
C *
C * CARDS OIL CHARACTERISTICS...(4F10.4,110)
C * PMU DYNAMIC VISCOSITY OF OIL
138
-------
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
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
PRHO
SPR
HLO
OLAG
IAT
CARD 7. ...
QP
TPB
TPE
DPA
IP
HS
CARD 7. ...
PVOL
DPA
DPL
CARD 7. ...
CARD 8. ...
COINI
FCV
HLC
CLAG
CARD 8
OIL DENSITY
RESIDUAL (TRAPPED) OIL SATURATION
OIL HALF LIFE (=0. IF NON-DEGRADIBLE )
ACCLIMATION TIME FOR BIODEGRADATION TO BEGIN
OIL INPUT TYPE 1=FLUX SPECIFIED
2=VOLUME/AREA SPECIFIED
. (FOR IAT=1) . . .OIL FLUX. . . (4F10. 4, 110, F10. 4) . . .
OIL FLUX (VOLUME/AREA/TIME)
OIL EVENT BEGINNING TIME
OIL EVENT ENDING TIME
DEPTH AT WHICH POLLUTANTS ARE APPLIED
PONDING TYPE 0=EXCESS OIL RUNS OFF SURFACE
1=VARIABLE (SEE NOTE 4.)
3=CONSTANT
CONSTANT HEAD FOR IP=3 CASES
. (FOR IAT = 2) . . .OIL VOLUME (3F10. 4)
OIL VOLUME/AREA INCORPORATED INTO THE SOIL
DEPTH AT WHICH POLLUTANTS ARE APPLIED
LOWER DEPTH OF INITIALLY POLLUTED ZONE
(FOR IAT = 3) CONSTANT HEAD. . (4F10. 4,110, F10 4)
***SAME AS FOR IAT=1, EXCEPT QP IS MEANINGLESS
.DISSOLVED CONSTITUENT. . . (6F10. 4)
INITIAL CONCENTRATION IN OIL (SEE NOTE 5.)
(SEE NOTE 5.)
CONSTITUENT HALF LIFE (=0. IF NONDEGRADING)
ACCLIMATION TIME FOR BIODEGRADATION TO BEGIN
.DISSOLVED CONSTITUENT. . . (6F10. 4)
PARTITIONING COEFFICIENTS:
XXKO
XXKV
XKOC
CARD 10...
DAIR
DWV
EVAP
TEMP
RH
XCF
CARD 11...
TM
DM
DTPR
CARD 12...
NTIMES
NZS
CARD 13...
PR (NTIMES)
LAST CARD.
OIL/WATER (CO = XXKW*CW)
DIMENSIONLESS HENRY'S LAW COEFF.
SOLID/ORGANIC CARBON
.VOLATILIZATION MODEL. . . (6F10. 4)
DIFFUSION COEFFICIENT FOR CONSTITUENT IN AIR
DIFFUSTION COEFFICENT FOR WATER VAPOR
WATER EVAPORATION RATE (VOL . /AREA/TIME )
TEMPERATURE IN DEGREES C
RELATIVE HUMIDITY
CONVERSION FACTOR TO GR/CC (SEE NOTE 2)
.SIMULATION PARAMETERS. . . (3F10.4)
SIMULATION ENDING TIME
MAXIMUM SOLUTION TIME STEP
MINIMUM TIME BETWEEN PRINTED TIME STEPS AND
MASS BALANCE CHECKS
.PROFILES AND HISTORIES. . . (215)
NUMBER OF PROFILES DESIRED (UP TO 10)
NUMBER OF HISTORIES DESIRED (UP TO 3)
. (5F10. 4/5F10. 4)
PROFILE TIMES
. (3F10 . 4 )
139
-------
C * HI(NZS) HISTORY DEPTHS
C *
C *
C * NOTES:
C * 1. ALL VARIABLE NAMES ARE IN ACCORDANCE WITH FORTRAN NAM-
C * .ING CONVENTIONS—NAMES BEGINNING WITH I THROUGH M ARE
C * INTEGERS, ALL OTHERS ARE REALS.
C * 2. WITH THE FOLLOWING EXCEPTIONS, ALL INPUT PARAMETERS MUST
C * BE GIVEN IN A CONSISTENT SYSTEM OF UNITS.
C * -THE TEMPERATURE MUST BE GIVEN IN DEGREES C
C * -THE CONVERSION FACTOR FROM THE CONSISTENT SET OF UNITS TO
C * GR/CC MUST BE GIVEN (PART OF THE VOLATILIZATION MODEL
C * USES AN EMPIRICAL FORMULA FOR VAPOR PRESSURE/DENSITY) .
C * 3. ZEROS SHOULD BE READ IN FIELDS PERTAINING TO UNUSED VALUES
C * 4. VARIABLE SURFACE PONDING IS SUPPORTED ONLY FOR CASES WITH
C * NO DISSOLVED CONSTITUENT PRESENT
C * 5. THE CONCENTRATION OF THE DISSOLVED CONSTITUENT CAN BE READ
C * IN IN ANY DESIRED UNITS (EG. MG/L) . A CONVERSION TO THE
C * CONSISTENT VOLUME UNIT MUST BE READ IN VARIABLE FCV IN
C * ORDER FOR THE CALCULATED MASS BALANCES TO BE CORRECT
C * E.G. CONCENTRATION IN MG/L, C-G-S UNITS, FCV = 0.0001
C * OUTPUT CONCENTRATIONS IN MG/L, MASS BALANCE IN MG/CM2
C *
CHARACTER NT*15
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /CONCI/ CLAG,COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /GAEQ/ CHS(2),DO,EE,FF,GA1,GA2,GA3,GA4,HS,IGA,IP,ZGA
COMMON /HIPR/ IHIST,NTIMES,PR(10),NZS,HI(5)
COMMON /MATR/ ETA, OMSOR, OMSWR, SR, SWR
COMMON /NAME/ NT (15)
COMMON /OCHR/ DLO,OLAG,SI(75,2), SO(75) ,TO(75), ZI(75,2)
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /RLTV/ KRF,XLAMB,A2A, A3A,A4A, A5A,A6A
COMMON /SIMU/ DM,DTPR,TM
COMMON /VAPO/ ABL,AKS,CN,CV(2),CVF(2),DAIR,XJV(2),ZV(2)
COMMON /WTER/ IRT/SWMAX,QW,WKS
EXTERNAL SLOPE
NUM(14) = NUM(14) + 1
C
C ***READ THE INPUT DATA
OPEN (5,FILE=IINPUT.KOPI,STATUS=IOLDI)
OPEN (6,FILE='OUTPUT.KOP',STATUS='UNKNOWN')
OPEN (7,FILE='RES',STATUS='UNKNOWN')
OPEN (8,FILE='PLOTI,STATUS='UNKNOWNI)
.OPEN (9, FILE=' PROF', STATUS=' UNKNOWN')
OPEN (10,FILE='HIST',STATUS='UNKNOWN')
DO 10 1=5,10
REWIND I
10 CONTINUE
C ***RUN TITLE
READ (5,9000) (NT(I),1=1,15)
WRITE (6,9100) (NT(I),1=1,15)
C ***CONSTANTS AND MATRIX PROPERTIES
READ (5,*) WKS,XLAMB,ETA,RHOS,SWR
WRITE (6,9101) WKS,XLAMB,ETA,RHOS,SWR
C ***WATER EVENT CHARACTERISTICS
140
-------
READ (5,*) WMU,WRHO,IRT,QW
IF (IRT.EQ.2) SWMAX = QW
WRITE (6,9102) WMU,WRHO,IRT,QW
C ***POLLUTANT CHARACTERISTICS
READ (5,*) PMU,PRHO,SPR,HLO,OLAG,IAT
WRITE (6,9103) PMU,PRHO,SPR,HLO,OLAG,IAT
IF (IAT.EQ.1.OR.IAT.EQ.3) THEN
C ***POLLUTANT FLUX EVENT CHARACTERISTICS
READ (5,*) QP,TPB,TPE,DPA,IP,HS
WRITE (6,9104) QP,TPB,TPE,DPA,IP,HS
ELSE IF (IAT.EQ.2) THEN
C ***POLLUTANT VOLUME EVENT CHARACTERISTICS
READ (5,*) FVOL,DPA,DPL
WRITE (6,9105) PVOL,DPA,DPL
END IF
OLAG = TPB4OLAG
C ***DISSOLVED CONSTITUENT PARAMETERS
READ (5,*) COTNI,FCV,HLC,CLAG
READ (5,*) XXKO,XXKV,XKOC
WRITE (6,9106) COINI,FCV,HLC,CLAG
WRITE (6,9107) XXKO,XXKV,XKOC
XXKS = XKOC/WRHO
CLAG = TPB+CLAG
C ***VOLATILIZATION PARAMETERS
READ (5,*) DAIR,DWV,EVAP,TEMP,RH,XCF
WRITE (6,9108) DAIR,DWV,EVAP,TEMP,RH,XCF
C ***SIMULATION PARAMETERS
READ (5,*) TM,DM,DTPR
DM = 1.1*DM
WRITE (6,9109) TM,DM,DTPR
C ***OUTPUT PROFILE TIMES AND HISTORY DEPTHS
READ (5,*) NTIMES,NZS
IF (NTIMES.GT.10.OR.NZS.GT.3) STOP 'INPUT PROFILE/HISTORY'
READ (5,*) (PR (I),1=1,NTIMES)
WRITE (6,9110) NTIMES, (PR(I),I=1,NTIMES)
C ***SORT THE OUTPUT TIMES
DO 15 J=l,NTIMES
IH = 0
DO 14 I=1,NTIMES-1
IF (PR(I).GT.PR(I+1)) THEN
IH = 1
PTEM = PR(I)
PR(I) = PR(I+1)
PR(I+1) = PTEM
END IF
14 CONTINUE
IF (IH.EQ.O) GO TO 16
15 ^CONTINUE
16 CONTINUE
READ (5,*) (HI(I),I=1,NZS)
WRITE (6,9111) NZS,(HI(I),I=1,NZS)
WRITE (6,9190)
C ***END OF INPUT DATA
C
SR = SPR + SWR
OMSWR = 1. - SWR
OMSOR = 1. - SPR
PKS = WKS*PRHO*WMU/WRHO/PMU
141
-------
IF (XXKO.EQ.O.) STOP 'INPUT: XXKW CANNOT BE ZERO1
XXKW = l./XXKO
XXKV = XXKV/XXKO
XSF1 = RHOS*XXKS/ETA/XXKO
IHIST = 0
DO 20 1=1,NZS
IF (HI(I).GE.DPA) IHIST = 1
20 CONTINUE
IF (IHIST.EQ.O) WRITE (6,9170)
CALL KRSET
IF (IAT.NE.2) THEN
DPL = DPA
DTP = TPE - TPB
PVOL = DTP*QP
END IF
WRITE (6,9180) PKS
IF (XXKV.GT.O.) THEN
C ***CALCULATION OF ATMOSPHERIC BOUNDARY LAYER
C ***LINSLEY, KOHLER & PAULHUS, 1975, PG. 35
C ***SATURATION VAPOR PRESSURE FORMULA PG. 35
ES = 33.8639*((0.00738*TEMP + 0.8072)**8 - .000019*
1 ABS(1.8*TEMP+ 48.) + 0.001316)
RWV = ES/(TEMP+273)/2876.*.622
ABL = DWV*RWV*XCF*(1.-RH)/2./EVAP/WRHO
WRITE (6,9184) ABL
END IF
IF (HLO.NE.O.) THEN
DLO = - ALOG(0.5)/HLO
WRITE (6,9181) DLO
END IF
IF (HLC.NE.O.) THEN
DLC = - ALOG(0.5)/HLC
WRITE (6,9182) DLC
END IF
RETURN
9000 FORMAT
9001 FORMAT
9002 FORMAT
9003 FORMAT
9004 FORMAT
9005 FORMAT
9007 FORMAT
9100 FORMAT
1
2
3
9101 FORMAT
1
2
3
4
5
9102 FORMAT
1
2
KINEMATIC OILY POLLUTANT TRANSPORT'/
(5A10/5A10/5A10)
(7F10.4)
(4F10.4,I10,F10.4)
(5F10.4,I10)
(2F10.4,I10,F10.4)
(5F10.4)
(1015)
(1H1/10X,50(1H*)/
10X,'KOPT
10X,50(1H*)/
3(10X,5A10/)//
15X,10HINPUT DATA/
15X,10(1H=)/)
(15X,30HCONSTANTS & MATRIX PROPERTIES.,10 (1H.)/
15X,30HSAT. HYDRAULIC CONDUCTIVITY = ,F10.4/
15X,30HPORE SIZE DIST. INDEX
15X,30HPOROSITY
15X,30HBULK DENSITY
15X,30HRESIDUAL WATER SATURATION
(15X,30HWATER EVENT CHARACTERISTICS.
15X,30HDYNAMIC VISCOSITY = ,F10.4/
15X,30HDENSITY = ,F10.4/
= ,F10.4/
= ,F10.4/
= ,F10.4/
,F10.4/)
142
-------
3
4
9103 FORMAT
1
2
3
4
5
6
9104 FORMAT
1
2
3
4
5
9105 FORMAT
1
15X, 30HRAIN TYPE : 1-FLUX 2-SAT. = ,1107
15X,30HWATER FLUX OR SATURATION = ,F10.4/)
(15X,35HPOLLUTANT EVENT CHARACTERISTICS ,5(1H.)/
15X,30HDYNAMIC VISCOSITY = ,F10.4/
15X,30HDENSITY = ,F10.4/
15X,30HRESIDUAL OIL SATURATION = ,F10.4/
15X,30HOIL HALF LIFE = ,F10.4/
15X,30HACCLIMATION TIME = ,F10.4/
15X,30HOIL LOADING TYPE = ,110)
(15X,30HFLUX LOADING RATE = ,F10.4/
15X,30HBEGINNING TIME = ,F10.4/
15X,30HENDING TIME = ,F10.4/
15X,30HDEPTH OF APPLICATION = ,F10.4/
15X,30HTYPE OF SURFACE PONDING = ,1107
15X,30HDEPTH OF PONDING = ,F10.4/)
(15X,30HOIL VOLUME/AREA = ,F10.4/
15X,30HDEPTH: TOP OF OIL LAYER = ,F10.4/
15X,30HDEPTH: BOTTOM OF OIL LAYER = ,F10.4/)
,F10.4/
= ,F10.
= ,F10.
= ,F10.
= ,F10.
= ,F10.
4/
4/
4)
4/
4/
,F10.4/)
2
9106 FORMAT (15X,32HDISSOLVED CONSTITUENT PARAMETERS,8(1H.)/
1
2
3
4
9107 FORMAT
1
2
9108 FORMAT
1
2
3
4
5
6
9109 FORMAT
1
2
3
9110 FORMAT
1
2
3
9111 FORMAT
1
2
9160 FORMAT
9170 FORMAT
9180 FORMAT
1
9181 FORMAT
9182 FORMAT
•9184 FORMAT
9190 FORMAT
END
SUBROUTINE START (NEQ,ZZ,ND)
DIMENSION NEQ(*),ZZ(ND,*)
CHARACTER NT*15
COMMON /CCHR/ CPART,CC(75,2),CO(75) ,TC(75),ZC(75,2)
COMMON /COUN/ IPR(60),NUM(60) ,SEC(60)
15X,30HINITIAL CONC. IN OIL
15X,30HCONVERSION TO CONST. UNITS
15X,30HCONSTITUENT HALF LIFE
15X,30HACCLIMATION TIME
(15X,30HOIL/WATER PARTITION COEF.
15X,30HDIMENSIONLESS HENRY'S LAW
15X,30HSOIL/WATER PARTITION COEF.
(15X,31HVOLATILIZATION MODEL PARAMETERS,9(1H.)/
15X,30HAIR-GAS DIFFUSION COEFF. = ,F10.4/
15X,30HAIR-WATER VAPOR DIFF. COEFF.= ,F10.4/
15X,30HEVAPORATION RATE = ,F10.4/
15X,30HTEMPERATURE DEG. C = ,F10.4/
15X,30HRELATTVE HUMIDITY = ,F10.4/
15X,30HCONVERSION TO GR/CC = ,F10.4/)
(15X, 30HSIMULATION PARAMETERS , 10 (1H.) /
15X,30HSIMULATION ENDING TIME = ,F10.4/
15X,30HMAXIMUM RKF TIME STEP = ,F10.4/
15X,30HMIN. TIME BETWEEN PRINTING = ,F10.4/)
(15X, 30HPROFILES AND HISTORIES ,10 (1H.) /
15X, 30HNUMBER OF PROFILES
15X,30HAT TIMES:
2(15X,5F6.2/))
(15X, 30HNUMBER OF HISTORIES
15X,30HAT DEPTHS:
15Xr5F6.2/)
(15X,30HRELATIVE VAPOR DENSITY
(15X,'***HISTORY DEPTHS ABOVE APPLICATION DEPTH***1)
(15X,21HCALCULATED PARAMETERS,19(1H.)/
15X,30HSAT. OIL CONDUCTIVITY
(15X,30HOIL DECAY RATE
(15X,30HCONSTITUENT DECAY RATE
,1
,F10.4)
= ,F10.
= ,F10.
= ,F10.
(15X,30HATM. BNDRY. LAYER THICKNESS = ,F10
(/15X,'***END OF INPUT DATA***1/)
4)
4)
4)
4)
143
-------
COMMON /CONCI/ CLAG,COINI,DLC,FCV,XSF1,XXKV,XXKW
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2) ,PZ(2)
COMMON /VAPO/ ABL,AKS,CN,CV(2),CVF(2),DAIR,XJV(2),ZV(2)
COMMON /HIPR/ IHIST,NTIMES,PR(10),NZS,HI(5)
COMMON /GAEQ/ CHS(2),DO,EE,FF,GA1,GA2,GA3,GA4,HS,IGA,IP,ZGA
COMMON /MATR/ ETA,OMSOR,OMSWR,SR,SWR
COMMON /NAME/ NT(15)
COMMON /OCHR/ DLO,OLAG, SI(75,2) ,SO(75) ,TO(75), ZI(75,2)
COMMON /OUTP/ ERCMX,EROMX,NPT,NREG
COMMON /POLL/ DPA,DPL,DTP,DZPMDT, DZPDDT,IAT,PKS,PVOL,QP,SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /RLTV/ KRF,XLAMB,EPS,EM1,EM2,CC1,REM1
COMMON /WTER/ IRT,SWMAX,QW,WKS
EXTERNAL SLOPE,PKE,WKE
NUM(15) = NUM(15) + 1
C
DO 1 1=1,51
ZZ(I,2) = DPA
ZZ(I,3) = DPA
1 CONTINUE
ZZ(1,1) = DPL
DO 2 1=1,2
PT(I) = TPB
PZ(I) = DPL
PS(I) = SPR
PF(I) = 0.
2 CONTINUE
C
C ***WATER COMPUTATIONS
C ***CONSTANT WATER SATURATION ABOVE THE WATERTABLE
C ***ASSUMED REPRESENTATIVE OF CLIMATIC CONDITIONS
CALL BISEC (WKE,0.5,1.,1.E-5,104,0,0.5,1.,SLIM,XXX,IE)
SAR = l.-SLIM
IF (IRT.EQ.l) THEN
C ***WATER FLUX SPECIFIED
SWMAX = SWR
IF (QW.GT.O.) SWMAX = SWR + OMSWR*((QW/WKS)**(1./EPS))
ELSE
C ***WATER SATURATION SPECIFIED
IF (SWMAX.GT.l.-SAR) SWMAX = l.-SAR
CALL WKE (0, SWMAX, 0,WKS,QW, 0., 0)
END IF
WRITE (6,9500) SWMAX,OW
IF (QW.GE.0.5*WKS) STOP 'START: WATER FILLS PORE SPACE1
C
C ***PRELIMINARY OIL COMPUTATIONS
NEQ(l) = 1
IGA = 0
IF (IAT.EQ.1.OR.IAT.EQ.3) THEN
SPMAX = l.-SWMAX-SAR
CALL PKE (0,SPMAX,0,SWMAX,XKPMX,0.,0)
CPART = COINI*QP/(QP + XXKW*QW)
ELSE IF (IAT.EQ.2) THEN
THPZ = DPL - DPA
SPMAX = PVOL/ETA/THPZ
IF (SPMAX.LT.SPR) WRITE (6,9505) THPZ,SPMAX
IF (SPMAX.GT.l.-SWMAX-SAR) THEN
SPMAX = l.-SWMAX-SAR
144
-------
PVOL = SPMAX*ETA*THPZ
WRITE (6,9510) SPMAX,PVOL
END IF
CALL PKE (0,SPMAX,0,SWMAX,QP,0.,0)
XKPMX = QP
Z = DPL
DZ = THPZ*0.1
DO 15 1=51,41,-!
ZI(I,2) = Z
ZI(I,1) = Z
ZZ(I,2) = Z
J = 51-1+1
ZC(J,2) = Z
ZC(J,1) = Z
ZZ(J,3) = Z
Z = Z - DZ
15 CONTINUE
CPART = PVOL*COINI/RD(SPMAX,SWMAX)/THPZ
END IF
IF (IE.EQ.O.AND.QP.GT.XKPMX) THEN
QP = XKPMX
PVOL = QP*DTP
CPART = COINI*QP/(QP + XXKW*QW)
END IF
IF (QP.LE.O..AND.IAT.NE.3) THEN
NEQ(l) = 0
NEQ(2) = 0
ELSE IF (QP.LE.XKPMX.AND.IAT.NE.3) THEN
^ ***KINEMATIC OIL FLOW
SL = SPR
IF (SL.GT.SPR) SL = QP/PKS
TOL = l.E-3
45 CONTINUE
CALL BISEC (PKE,SL/1.,TOL,110,1,QP,SWMAX,SPMAX,XX,IE)
CALL PKE (0,SPMAX,0,SWMAX,XKPMX,0.,0)
: ***ASSURE THAT SPMAX GIVES FLUX CLOSE TO LOADING RATE
IF (ABS(XKPMX-QP).GT.5.E-5) THEN
IF (TOL.LT.l.E-10) STOP 'OIL-TOLERANCE1
TOL = TOL*0.1
GO TO 45
END IF
DZPMDT = QP/ETA/SPMAX
ELSE IF (IAT.EQ.3.0R.QP.GT.XKPMX) THEN
: ***PRESSURIZED OIL FLOW
EE = XKPMX*(l.-SPMAX*ETA)
FF = XKPMX*QP
ZZ(1,1) = DPA+0.001
IGA = 1
: co(i) = COINI
CHS(l) = 0.000001
CHS(2) = 0.
IF (IP.EQ.3) THEN
GA1 = SPMAX*ETA/XKPMX
GA2 = ALOG( HS )
GA3 = HS - DPA
GA4= DPA - (HS-DPA)*ALOG(HS)
ZZ(1,1) = DPA
: NEQ(l) = 0
145
-------
END IF
DO = 0.
IF (DLO.GT.O.) THEN
DO = DLO
DLO = 0.
END IF
END IF
IF (IAT.EQ.1.AND.IGA.EQ.O.OR.(IGA.EQ.l.AND.IP.EQ.1)) THEN
WRITE (6,9520) QP,PVOL
IF (COINI.GT.O.) WRITE (6,9525) QP*DTP*COINI*FCV
END IF
CALL SLOPE (0,SPMAX,0,SWMAX,DZPDDT,0.,0)
DZPDDT = DZPDDT/ETA
ZGA = DPA + QP*DTP/SPMAX/ETA
PS(1) = SPMAX
NEQ(2) = 0
IF (COINI.GT.O..AND.DLO.GT.O.) THEN
CODM = XXKW*CPART*RD(SPMAX,SWMAX)/RD(0.,SWMAX)
WRITE (6,9530) CODM
END IF
C
C ***STARTING VALUES FOR CHARACTERISTICS
IL = 40
IF (IAT.NE.2) IL = 51
DO 4 J=l,2
DO 3 1=1,IL
ZI(I,J) = DPA
ZC(I,J) = DPA
SI(I,J) = SPMAX
CC(I,J) = CPART
3 CONTINUE
4 CONTINUE
DT = DTP*0.1
T = TPB - DT
DO 20 1=51,41,-!
T = T + DT
TO(I) = T
SO(I) = SPMAX
J = 51-1+1
TC(J) = T
CO(J) = CPART
NEQ(3) = J
20 CONTINUE
IF (COINI.EQ.O.) NEQ(3) = 0
IF (COINI.GT.O..AND.IGA.EQ.1) THEN
IF (IP.EQ.l) STOP 'FEATURE IP=1 NOT SUPPORTED FOR COINI.GT 0'
C0(l) = COIN!
~ CC(1,2) = COIN!
CC(1,1) = COINI
C ***ADD MORE CONSTITUENT CHARACTERISTICS FOR GREEN-AMPT CASES
DO 40 I=NEQ(3),2,-1
TC(I+4) = TC(I)
CO(1+4) - CO(I)
40 CONTINUE
DT = DT/5.
T = TC(1)
DO 50 1=2,5
T = T + DT
146
-------
TC(I) = T
CO(I) = CPART
50 CONTINUE
NEQ(3) = NEQ(3) + 4
END IF
IF (XXKV.GT.O.) THEN
CN = 10./3.
IF (IAT.EQ.2) NEQ(l) = 2
ZZ(2,1) = DPA
END IF
DS = (SPMAX-SPR)/40.
SO = SPR - DS
DO 65 1=1,40
SO = SO + DS
TO (I) = TPE
SO (I) = SO
SI(1,1) = SO
SI (I, 2) = SO
65 CONTINUE
CALL INQSET (SLOPE,SPR,1.-SWMAX,SWMAX,50,WP)
IF (DLO.GT.O..AND.QP.GT.O.) NEQ(2) = 51
: ***HEADINGS FOR OUTPUT FILES
WRITE (6,9580) (NT(II),11=1,15)
WRITE (6,9585)
IF (IGA.EQ.l) THEN
WRITE (6, 9588) 1,TPB, ZZ (1,1),CHS (2) ,SPMAX,PF (2), 0., ZC (1,2), ZV(2),
1 CC(1,2)*XXKW,0.,0.
9588 FORMAT (1X,I5,12F10.4)
ELSE
WRITE (6,9588)1,TPB,ZZ(1,1),0.,SPMAX,QP,0.,ZC(1,2),ZV(2),
1 CC(1,2)*XXKW,0.,0.
END IF
WRITE (10,9100) (NT(I),1=1,15),(HI(I),I=1,NZS)
IF (NZS.EQ.l) WRITE (10,9102)
IF (NZS.EQ.2) WRITE (10,9103)
IF (NZS.EQ.3) WRITE (10,9104)
9500 FORMAT (15X,30HWATER SATURATION = ,F10.4/
1 15X,30HWATER FLUX = ,F10.4)
9505 FORMAT (15X,30HTHICKNESS OF POLLUTED ZONE = ,F10.4/
2 15X,30HMAXIMUM OIL SATURATION = ,F10.4)
9510 FORMAT (15X,34HOIL SAT. REDUCED TO MAX. ALLOWABLE/
1 15X,30HADJUSTED OIL SATURATION = ,F10.4/
2 15X,30HADJUSTED OIL VOLUME/AREA = ,F10.4)
9520 FORMAT (15X,30HPOLLUANT VOLUME FLUX = ,F10.4/
1 15X,30HTOTAL OIL LOADING, VOL/AREA = ,F10.4)
9525 FORMAT (15X,30HTOTAL CONSTITUENT LOADING = ,F10.4)
9530 FORMAT (15X,30HMAX. CONC. DUE TO OIL DECAY = ,F10.4)
9580 FORMAT (1H1//28X,50(1H*)/
1 33X,'LOCATION OF THE OIL AND CONSTITUENT FRONTS1/
' 2 28X,50(1H*)/
4 3(28X,5A10/)/)
9585 FORMAT (/
1
2
2
3
4
16X, 23X,'OIL',24X,30X,'CONSTITUENT'/
18X,48(1H-),IX,49(1H-)/
IX, ' STEP',4X, ITIMEI,2X,3X, 'DEPTH',2X, 4X, 'HEAD',2X,
1 SATURATION', 4X, ' FLUX', 2X, IX, ' VOL/AREA', IX, 9X, ' DEPTHS'
, 5X, 5X, ' CONC', IX, IX,'MASS/AREA', IX, ' CUMULATIVE'/
147
-------
5 26X, 10X, SOX, 4X, 'LOWER1, IX, 4X, 'UPPER1,20X, 'VOLATILIZATION' /
6 1X,115(1H=)/)
9100 FORMAT (//10X,50(1H*) /
1 10X,' SATURATION AND CONCENTRATION HISTORIES1/
2 10X,50(1H*)/
3 3(10X,5A10/)/
4 1X,10H TIME ,12X,'DEPTH(S)'/
5 10X,3(10X,F10.4,10X))
9102 FORMAT (10X,(4X,'SAT',3X,4X,'FLUX',2X,3X,'CONC',3X)/
1 1X,9(1H=),(1X,29(1H=)))
9103 FORMAT (10X,2 (4X,'SAT', 3X, 4X,'FLUX',2X, 3X,'CONC', 3X) /
1 1X,9(1H=),2(1X,29(1H=)))
9104 FORMAT (10X, 3 (4X, ' SAT' , 3X, 4X, 'FLUX', 2X, 3X, ' CONC', 3X) /
1 IX,9(1H=),3(IX,29(1H=)))
RETURN
END
SUBROUTINE COUNT
CHARACTER NT*15
COMMON /COUN/ IPR(60) ,NUM(60),SEC(60)
COMMON /FRNT/ ISE,PF(2),PS(2),PT(2),PZ(2)
COMMON /NAME/ NT(15)
COMMON /OUTP/ ERCMX,EROMX,NPT,NREG
CHARACTER CN(60)*6
C *-+-+$+—«--*-+-+$ * $ * $ * $ * $ * $—
xx*xxxx
DATA (CN(I),1=1,29) /'MAIN','PFRONT1,'PCHK','GACHK','GAHS',
1 'RD1, 'CVEL','CONC','OCP','ZMASS', 'OIL1,'HIST',
2 'SFLUX', 'INPUT', 'START','COUNT', 'MRKF12', 'GLQ1, 'INQUAD',
3 'BINARY1,'BISEC1,'WKE',' ','PKE1,'SLOPE',' ','KRSET',
4 'INQSET',1 '/
C
NUM(16) = NUM(16) + 1
C ***RUN INFORMATION
WRITE (6,9500) (NT (II) ,11=1,15)
WRITE (6,9501)
Nl = 0
NEND = 28
90 CONTINUE
IF (Nl.GT.NEND) GO TO 100
Nl = Nl + 1
IF (CN(Nl).EQ.1 ') GO TO 90
N2 = Nl + 1
95 CONTINUE
IF (CN(N2) .EQ_.' '.AND.N2.NE.NEND) THEN
N2 = N2 + 1
GO TO 95
END IF
IF (CN(N2).NE.' ') THEN
WRITE (6,9503) N1,CN(N1) ,NUM(N1) ,N2,CN(N2) ,NUM(N2)
ELSE
WRITE (6,9503) N1,CN(N1),NUM(N1)
END IF
Nl = N2
GO TO 90
100 CONTINUE
XI = FLOAT(NUM(2))
X2 = FLOAT(NPT+NREG)
X3 = FLOAT(NPT)
148
-------
WRITE (6,9820)
WRITE (6,9900)
STOP 'END1
NPT+NREG,NPT,NREG,NUM(2),X1/X2,X1/X3,FLOAT(ISE)/Xl
C
C
C
C
C
***OUTPUT FORMATS
FORMAT (1H1//10X,40HUNSTEADY IMMISCIBLE GROUNDWATER POLLUTAN,
14HT INFILTRATION/
10X,20X,15HRUN INFORMATION/
10X,54(1H*)/
3(12X,5A10/)/)
(10X,20HROUTINE CALLS,14X,20HROUTINE CALLS/
10X,20(1H=),14X,20(1H=)/)
(8X,14,3X,A6,3X,16,12X,14,3X,A6,3X,16)
(//15X,1RUNGE-KUTTA-FEHLBERG SOLVER PERFORMANCE'/
10X,54(1H=)//
10X,'NUMBER OF TIME STEPS TAKEN
10X,'ACCEPTED STEPS
10X,'REJECTED STEPS
10X,'FUNCTION EVALUATIONS
10X,'AVERAGE NO. EVALUATIONS PER STEP
10X,'AVERAGE NO. EVALUATIONS PER ACCEPTED STEP
10X,'AVERAGE NO. EQUATIONS SOLVED PER EVAL.
(10X,'HISTORIES ARE IN DEFAULT FILE HIST1/
1OX,'PROFILES ARE IN DEFAULT FILE PROF'/
//10X,50(1H*)/
10X,'SUCESSFUL EXECUTION'/
10X,50(1H*))
9500
1
2
3
4
9501 FORMAT
1
9503 FORMAT
9820 FORMAT
1
2
3
4
5
6
7
8
,F10.4/
,F10.4/
,F10.4/)
9900 FORMAT
1
2
3
4
END
SUBROUTINE WKE (IS,SW,N,XKS,XKR,SV,IE)
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
COMMON /MATR/ ETA,OMSOR,OMSWR,SR, SWR
COMMON /POLL/ DPA,DPL,DTP,DZPMDT,DZPDDT, IAT,PKS,PVOL,QP, SPMAX,SPR,
1 THPZ,TPB,TPE,TPR,WP(50),XKPMX,XK22
COMMON /RLTV/ KRF,XLAMB,EPS,EM1,EM2,CC1,REM1
***RELATIVE PERMEABILITY VIA THE BROOKS-COREY INTEGRATION
***OF THE WYLLIE AND GARDNER BURDINE EQUTIONS
NUM(22) = NUM(22) + 1
XKR = 0.
IF (SW.LE.SWR) RETURN
XKR = XKS*((SW - SWR)/OMSWR)**EPS
RETURN
ENTRY SWKE (IS,SW,N,XKS,SKR,SV,IE)
***DERIVATIVE OF THE RELATIVE PERMEABILITY W.R.T. SW
NUM(23) = NUM(23) + 1
SKR = 0.
.IF (SW.LE.SWR) RETURN
SKR = EPS*XKS*((SW - SWR)**EM1) *CC1
RETURN
ENTRY PKE (IS,SP,N,SW,PK,SVOL,IEND)
***RELATIVE PERMEABILITY FOR OIL
SPF = (SP - SPR)/OMSOR
IF (SPF.LT.O.) SPF = 0.
Tl = (SP + SW - SWR)/OMSWR
IF (Tl.GT.O.) Tl = T1**EM2
T2 = (SW - SWR)/OMSWR
149
-------
IF (T2.GE.O.) T2 = T2**EM2
PK = PKS*SPF*SPF*(Tl - T2)
NUM(24) = NUM(24) + 1
RETURN
C
ENTRY SLOPE (IS,SP,N,SW,DKDSP,SVOL,IEND)
C ***DERIVATIVE OF KR (OIL) WRT SP
T3 = 0.
EM3 = EM2 - 1.
SPF = (SP - SPR)/OMSOR
IF (SPF.LT.O.) SPF = 0.
Tl = (SP + SW - SWR)/OMSWR
IF (Tl.GE.O.) THEN
T3 = T1**EM3
Tl = T1**EM2
END IF
T2 = (SW - SWR)/OMSWR
IF (T2.GE.O.) T2 = T2**EM2
DKDSP = 2.*SPF*(T1-T2)/OMSOR + EM2*SPF*SPF*T3/OMSWR
DKDSP = DKDSP*PKS
NUM(25) = NUM(25) + 1
RETURN
C
ENTRY KRSET
C ***SET UP ROUTINE
NUM(27) = NUM(27) + 1
EPS = (2. + 3.*XLAMB)/XLAMB
EMI = EPS - 1.
EM2 = EPS - 2.
EM3 = EPS - 3.
CC1 = l./(OMSWR**EPS)
REM1 = l./EML
RETURN
END
SUBROUTINE MRKF12 (FNC,CHK,TX,TF,DM,ER,NEQ,NS,NR,YY,ND,ET)
\^
C * *
C * MRKF12 (MULTIPLE) RUNGE KUTTA FELHBERG FIRST-SECOND *
C * ORDER ORDINARY DIFFERENTIAL EQUATION SOLVER *
C * *
C * *
C * PURPOSE THIS SUBROUTINE SOLVES A SYSTEM OF ORDINARY *
C * DIFFERENTIAL EQUATIONS OF THE FORM *
C * *
C * DY/DT = F(Y,T) . *
C * *
C * THE FUNCTIONS F(Z,T) ARE CONTAINED IN THE ROUTINE FNC. *
C .* THE EQUATIONS ARE IN 3 GROUPS, SOME OF WHICH MAY NOT BE *
C * SOLVED AT ALL TIMES DURING THE SOLUTION PROCESS. *
C * ROUTINE CHK CONTAINS CHECKS ON THE SOLUTION AND SWITCHES *
C * TO CONTROL THE EXECUTION. THE FELHBERG FEATURE ALLOWS *
C * THE TIME STEP TO VARY. *
C * *
C * INPUT ARGUMENTS *
C * FNC = SUBROUTINE WHICH CONTAINS FUNCTIONS F(Y,T) *
C * CHK = SUBROUTINE WHICH CONTAINS CHECKS AND CONTROLS *
C * TX STARTING TIME . *
C * TF ENDING TIME *
150
-------
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
*
it
*
*
*
*
*
*
*
*
*
*
*
*
*
*
• *
*
*
*
*
*
*
*
**•
DM = MAXIMUM TIME STEP
ER • = ERROR TOLERANCE PER STEP
NEQ = ARRAY CONTAINING NUMBER OF EQUATIONS PER GROUP
NS = STEP NUMBER (STARTING)
YY = INITIAL CONDITION
ND = ROW DIMENSION OF YY
OUTPUT ARGUMENTS
DM = MAXIMUM TIME STEP
NEQ = NUMBER OF EQUATIONS PER GROUP
NS = STEP NUMBER (ENDING)
NR = NUMBER OF REJECTED STEPS
YY = SOLUTION AT ENDING TIME
ET = ENDING TIME OF SOLUTION
AUTHOR JIM WEAVER
THE UNIVERSITY OF TEXAS At AUSTIN
REFERENCE: FEHLBERG, 1969, NASA TR R-315
SINGLE PRECISION, ANSI STD. X3. 9-1978 FORTRAN 77
REQUIRED ROUTINES... SEE INPUT ARGUMENTS
X-21-87
***********************************************************
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
***
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
C ***CAUTION: THE FIRST DIMENSION OF XKO,XKL,XK2 & YI MUST MATCH ND
DIMENSION DTR(20),IRR(20),NEQ(*) ,XKO(55,3),XK1(55,3),XK2 (55,3),
1YI(55,3),YY(ND,*)
PARAMETER (Al=0.5, B10=A1,B20=l./256.,B21=255./256.,C3=l./512.)
DATA IRR,NDT /20*0,13/
NUM(17) = NUM(17) + 1
IF (ND.NE.55) STOP 'RUNGE-KUTTA'
C
C ***INITIAL VALUES
1C = 0
N = NN
DT = 0.0001
TN = TX
IU =0
C
1 CONTINUE
C ***SAVE INITIAL CONDITION IN CASE PROPOSED TIME STEP IS
C ***TOO LARGE
DO 25 J=l,3
DO 20 1=1,NEQ(J)
YI(I,J) = YY(I,J)
20 CONTINUE
25 CONTINUE
30 CONTINUE
DO 35 1=1, NOT
DTR(I) = 1.E20
35 CONTINUE
C ***INTEGRATE ONE STEP
C ***SET SOME CONSTANTS FOR THIS TIME STEP
AID = DT * Al
B10D = DT * BIO
B20D = DT * B20
151
-------
B21D --= DT * B21
C3D = DT * C3
XT = TN
C
C ***BEGIN SERIES OF FUNCTION EVALUATIONS FOR CERTAIN
C ***VALUES OF TN AND YY
IF (IU.EQ.O) CALL FNC (XT,YI,XKO,ND,NEQ)
C
C ***THE SECOND EVALUATION
TN = XT + AID
DO 45 J=l,3
DO 40 I=1,NEQ(J)
IF (XKO(I,J).EQ.O.) GOTO 40
YY(I,J) = YI(I,J) + B10D*XKO(I,J)
40 CONTINUE
45 CONTINUE
CALL FNC (TN,YY,XK1,ND,NEQ)
C
C ***THE THIRD EVALUATION
TN = XT + DT
DO 55 J=l,3
DO 50 I=a,NEQ(J)
YY(I,J) = YI(I,J) + B20D*XKO(I,J) + B21D*XK1(I,J)
50 CONTINUE
55 CONTINUE
CALL FNC (TN,YY,XK2,ND,NEQ.)
C
C ***DETERMINE NEW TIME STEP FROM THE FUNCTION EVALUATIONS
IE = 0
IE2 = 0
TR = 1000.
DO 75 J=l,3
DO 70 I=1,NEQ(J)
IF (XKO(I, J) .EQ.XK2(I,J)) GO TO 70
IF (XKO(I,J).EQ.O..AND.XK2(I,J).GT.O.) GO TO 70
TE=ABS( C3D*(XKO(I,J) -XK2(I,J)) )
BB = TE / (AMAX1 (1., ABS ( YY(I,J) ) ) )
IF (BB.GT.ER) THEN
C ***STEP REJECTED
IE = 1
DTR(l) = DT*.5
GO TO 76
END IF
TS = DT * ( (ER + ER*ABS(YY(I/J) ) )/TE ) ** 0.25
IF (TS.LT.TR) TR = TS
70 CONTINUE
75 -CONTINUE
76 CONTINUE
IE2 = 0
DTI = DT
CALL CHK (TN,YY,YI,ND,NEQ,DM/DT1,DTR,NS,IC,IRR,IE2,IE)
ET = TN
IF (IC.EQ.l) RETURN
IF (IE2.EQ.1) IU = 0
IF (IE.EQ.1.0R.IE2.EQ.1) THEN
C ***REJECT STEP
NR = NR + 1
TN = XT
152
-------
; ***SELECT SMALLEST STEP SIZE LIMITATION
DT = 1.E19
DO 80 11=1,NOT
IF (DTR(II).LT.DT.AND.DTR(II).GT.O.) THEN
DT = DTK(II)
IN = II
END IF
80 CONTINUE
IRR(IN) = 1
IF (DT.GE.1.E19) STOP 'RKF: DT=1.E19'
GO TO 30
END IF
NS = NS + 1
; ***SET NEW TIME STEP
IU = 1
IF (DT1.LT.DT) TR = DT1*1.25
DT = TR * 0.80
IF (DT.GT.DM) DT = DM
: ***THIRD EVALUATION BECOMES FIRST EVALUATION
DO 95 J=l,3
DO 90 I=1,NEQ(J)
XKO(I,J) = XK2(I,J)
90 CONTINUE
95 CONTINUE
GO TO 1
END
SUBROUTINE GLQ (FNC,A1,NEQ,A,B,IW,SUM)
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
**i
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
-*
*
*
*
it
*
*
*
-------
DATA CS /.129485,.279705, .381830, .417959, .381830,.279705,.1294857
NUM(18) = NUM(18) + 1
Cl = (B - A)*0.5
C2 = (B + A)*0.5
C
C ***FUNCTION EVALUATIONS AT GAUSS POINTS
SUM = 0.
DO 10 1=1,7
Y = CY(I)*C1 + C2
CALL FNC (Y,A1,A2,NEQ,VAL)
IF (IW.EQ.l) WRITE (9,9101) Y,VAL
SUM = SUM + VAL*CS(I)
10 CONTINUE
SUM = SUM*C1
9101 FORMAT (15X,F10.4,5X,F10.4)
RETURN
END
SUBROUTINE INQUAD (X,NX,AMIN,AMAX,VA,RE)
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
**>
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
******* **************************************************
INQUAD
PURPOSE THIS SUBROUTINE PERFORMS INVERSE QUADRATIC
INTERPOLATION ON A SET OF PREDETERMINED NUMBERS TO
INVERT A NONLINEAR FUNCTION. INQSET IS CALLED FIRST
TO SET UP THE ARRAY OF FUNCTION VALUES. SUBSEQUENT
CALLS TO INQUAD GIVE THE INVERSE FOR SPECIFIED VALUES
OF THE DEPENDENT VARIABLE
INPUT ARGUMENTS
X = ARRAY OF FUNCTION VALUES FROM INQSET
NX = NO. OF X'S
AMIN = MINIMUM INDEPENDENT VARIABLE VALUE
AMAX = MAXIMUM INDEPENDENT VARIABLE VALUE
VA = VALUE OF DEPENDENT VARIABLE
OUTPUT ARGUMENTS
RE = VALUE OF INDEPENDENT VARIABLE: VA = F(RE)
AUTHOR JIM WEAVER
THE UNIVERSITY OF TEXAS AT AUSTIN
REFERENCES: ABRAMOWITZ AND STEGUN, PG. 882.
SINGLE PRECISION, ANSI STD. X3. 9-1978 FORTRAN 77
REQUIRED ROUTINES... INQSET
****** ***************************************************
***
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
***
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
REAL X(*)
DATA I /!/
NUM(19) = NUM(19) + 1
DX = (AMAX - AMIN)/FLOAT ( NX - 1 )
CALL BINARY (VA,X,NX,1,1)
IF (I.EQ.l.OR.I.EQ.NX-l) THEN
P= (VA—
154
-------
ELSE IF (I.LT.NX) THEN
A = X(I+1) - 2.*X(I) + X(I-l)
B = X(I+1) - X(I-l)
C = 2.*(X(I) - VA)
Q = B*B - 4.*A*C
IF (Q.LT.O.) STOP 240
Q = SQRT( Q )
P = (Q - B)*0.5/A
ELSE
WRITE (7,1010) I,NX,AMIN,AMAX,VA,X(1),X(NX)
1010 FORMAT (IX,'INQUAD: I,NX,MIN,MAX,VA,X(1),X(NX)•,215,5F10.4)
STOP 241
END IF
RE = AMIN + (FLOAT(1-1) + P) *DX
RETURN
•i
ENTRY INQSET (FNC,AMIN,AMAX,PARA,NX,X)
C
C
C
C
c
C
c
c
c
c
c
c
c
c
c
**1
*
*
*
*
*
*
*
*
it
*
*
*
*
**1
k*******************************************************-J
PURPOSE THIS SUBROUTINE SETS UP ARRAYS FOR INQUAD
INPUT ARGUMENTS
FNC = SUBROUTINE NAME OF A NONLINEAR FUNCTION
AMIN = MINIMUM VALUE OF INDEPENDENT VARIABLE
AMAX = MAXIMUM VALUE OF INDEPENDENT VARIABLE
PARA = REAL PARAMETER PASSED TO FNC
NX = NUMBER OF DEPENDENT VALUES TO PRODUCE
OUTPUT ARGUMENTS
X = ARRAY OF DEPENDENT VARIABLES
****** **************************************************!
t***
*
*
*
*
*
*
*
*
*
*
*
*
*
k***
NUM(28) = NUM(28) + 1
SS = AMIN
DX = (AMAX - AMIN)/FLOAT ( NX - 1 )
DO 1 1=1,NX
CALL FNC (0,SS,0,PARA,X(I),0.,0)
SS = SS + DX
1 CONTINUE
RETURN
END
SUBROUTINE BINARY (U,X,NX,INCX, IN)
C
C * *
C * BINARY • *
C * *
C * *
C * PURPOSE THIS SUBROUTINE PERFORMS A BINARY SEARCH ON *
C * THE ARRAY X TO FIND THE INDEX OF THE X VALUE WHICH LIES *
C * JUST BELOW U *
C * *
C * INPUT ARGUMENTS *
C * U = POINT SEARCHED FOR . *
C * X = ARRAY OF VALUES *
C * NX = NUMBER OF X V7.LUES *
C * INCX = INCREMENT (MUST = 1) *
C * *
C * OUTPUT ARGUMENTS *
155
-------
C * IN INDEX OF X VALUE JUST BELOW U *
C * *
C * *
C * AUTHOR JIM WEAVER *
C * THE UNIVERSITY OF TEXAS AT AUSTIN *
C * *
C * REFERENCE: FORSYTHE, MALCOLM AND MOLER, 1977, PG. 76 *
C * SINGLE PRECISION, ANSI STD. X3.9-1978 FORTRAN 77 *
C * REQUIRED ROUTINES...NONE *
C * *
COMMON /COUN/ IPR(60),NUM(60),SEC(60)
REAL X(*)
DATA I/I/
C
NUM(20) = NUM(20) + 1
IF (INCX.NE.l) STOP 'BINARY-INDEX1
C ***X'S ORDERED BY INCREASING VALUE
IF (X(l).LT.X(NX)) THEN
IF (I.GT.NX) 1=1
IF (U.LT.X(I),OR.X(I+1).LT.U) THEN
IF (U.LT.X(I)) THEN
J = I
1 = 1
ELSE
1 = 1 + 1
J = NX
END IF
2 CONTINUE
K = (I + J)/2
IF (U.LT.X(K)) J = K
IF (U.GE.X(K)) I = K
IF (J.GT.I+1) GO TO 2
END IF
IN = I
IF (IN.EQ.NX) IN = IN - 1
ELSE IF (X(l).GT.X(NX)) THEN
C ***X'S ORDERED BY DECREASING VALUES
IF (I.GT.NX) 1=1
IF (U.LT.X(I+1).OR.X(I).LT.U) THEN
IF (U.LT.X(I+1)) THEN
J = NX
1 = 1 + 1
ELSE
J = I
1 = 1
END IF
5 CONTINUE
K - (I + J)/2
IF (U.LT.X(K)) I = K
IF (U.GE.X(K)) J = K
IF (J.GT.I+1) GO TO 5
END IF
IN = I
END IF
RETURN
END
SUBROUTINE BISEC (FNC,ZB,ZE,CONV,IS, NN,SV,TM,ZM,VM,IE)
156
-------
c * • *
C * BISEC *
C * *
C * *
C * PURPOSE THIS SUBROUTINE USES THE BISECTION METHOD TO *
C * FIND THE ZERO OF A FUNCTION *
C * *
C * INPUT ARGUMENTS *
C * FNC = SUBROUTINE NAME CONTAINING FUNCTION *
C * ZB LOWER LIMIT OF SEARCH *
C * ZE = UPPER LIMIT OF SEARCH *
C * CONV = CONVERGENCE TOLERANCE *
C * IS MAXIMUM NUMBER OF ITERATIONS *
C * NN = INTEGER PASSED TO FNC *
C * SV FUNCTION VALUE TO BE MATCHED BY BISEC ITERATION *
C * TM REAL NUMBER PASSED TO FNC *
C * *
C * OUTPUT ARGUMENTS *
C * ZM BISEC RESULT (ZB.LE.ZM.LE.ZE) *
C * VM FUNCTION VALUE AT ZM (SHOULD MATCH SV) *
C * IE = ERROR CODE IE=-1: BISEC FAILED *
C * *
C * *
C * AUTHOR JIM WEAVER *
C * THE UNIVERSITY OF TEXAS AT AUSTIN *
C * *
C * REFERENCE: FORSYTHE, MALCOLM & MOLER, 1977, PG. 157 *
C * SINGLE PRECISION, ANSI STD. X3.9-1978 FORTRAN 77 *
C * REQUIRED ROUTINES...SEE INPUT ARGUMENTS *
C * THANKS TO SHERYL FRANKLIN *
\^
COMMON /COUN/ IPR(60),NUM(60), SEC(60)
NUM(21) = NUM(21) + 1
Tl = 1.
T2 = 1.
XF = 1.
ZO = ZB .
ZF = ZE
CALL FNC (IS,ZO,NN,TM,V1,SV,IE)
IF (Vl.GT.SV) XF = -1.
Zl = VI - SV
SW = SV "
CALL FNC (IS,ZF/NN,TM,V2,SV,IE)
IF (SIGN(T1,Z1).EQ.SIGN(T2,V2-SV)) THEN
IE = -1
WRITE (7,7020) IS,SIGN(T1,V1-SV),SIGN(T2,V2-SV)
WRITE (7,7025) T1,V1,SW,T2,V2,SV
7025 FORMAT (IX,'BISEC: T1,V1,SW,T2,V2,SV'/ 6G10.4)
7020 FORMAT (IX,'BISEC: IS,SIGNS',I5,2F10.4)
RETURN /
END IF
DO 1 1=1,IS
ZM = ZO + (ZF - Z0)*0.5
CALL FNC (IS,ZM,NN,TM,VM,SV,IE)
XM = XF*(VM - SV)
IF (ABS(XM) .LE.CONV) RETURN
IF (XM.GE.O.) THEN
157
------- |