EPA/600/R-98/046
PARAMETER ESTIMATION OF TWO-FLUID CAPILLARY PRESSURE-
SATURATION AND PERMEABILITY FUNCTIONS
by
Jan. W. Hopmans, Mark E. Grismer, and J. Chen
Department of Land, Air and Water Resources
University of California, Davis, CA 95616
Y.P. Liu
Department of Plant and Soil Science
Alabama A&M University, Normal AL 35762
Cooperative Agreement
CR822204
Project Officer
Bob K. Lien
Subsurface Protection and Remediation Division
National Risk Management Research Laboratory
Ada, Oklahoma 74820
National Risk Management Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Cincinnati, OH 45268
-------
NOTICE
The U.S. Environmental Protection Agency through its Office of Research and Development
partially funded and collaborated in the research described here under Cooperative Agreement No.
CR-822204 to the University of California, Davis. It has been subjected to the Agency's peer and
administrative review and has been approved for publication as an EPA document. Mention of trade
names or commercial products does not constitute endorsement or recommendation for use.
All research projects making conclusions or recommendations based on environmentally related
measurements and funded by the Environmental Protection Agency are required to participate in the
Agency Quality Assurance Program. This project was conducted under an approved Quality
Assurance Project Plan. The procedures specified in this plan were used without exception.
Information on the plan and documentation of the quality assurance activities and results are
available from the Principal Investigator.
-------
FOREWORD
The U.S. Environmental Protection Agency is charged by Congress with protecting the Nation's
land, air and water resources. Under a mandate of national environmental laws, the Agency strives
to formulate and implement actions leading to a compatible balance between human activities and
the ability of natural systems to support and nurture life. To meet these mandates, EPA's research
program is providing data and technical support for solving environmental problems today and
building a science knowledge base necessary to manage our ecological resources wisely, understand
how pollutants affect our health, and prevent or reduce environmental risks in the future.
The National Risk Management Research laboratory is the Agency's center for investigations of
technological and management approaches for reducing risks from threats to human health and the
environment. The focus of the Laboratory's research program is on methods for the prevention and
control of pollution to air, land, water, and subsurface resources; protection of water quality in
public water systems; remediation of contaminated sites and ground water, and prevention and
control of indoor air pollution. The goal of this research effort is to catalyze development and
implementation of innovative, cost-effective environmental technologies; develop scientific and
engineering information needed by EPA to support regulatory and policy decisions; and provide
technical support and information transfer to ensure effective implementation of environmental
regulations and strategies.
Earlier work by the principal investigators has shown that parameter optimization by inverse
modeling can be used to estimate the soil water retention and unsaturated hydraulic conductivity
functions of soils containing air and water only. The research presented in this report focuses on the
application of this method to determine capillary pressure and permeability functions in multi-fluid
soil systems (air-water, air-oil and oil-water) using data from the multi-step outflow method. The
term multi-fluid is used here to indicate what is traditionally defined as multiphase. Whereas soil
water retention and unsaturated hydraulic conductivity are generally used in the soils literature for
air-water systems only, the definition of these relationships in general multi-fluid soil systems
requires use of the capillary pressure and permeability terminology instead. The authors conclude
that the applied inverse model is well posed for the investigated multi-fluid soil systems, and that the
parameter optimization yields accurate capillary pressure-saturation and permeability functions.
Clinton W. Hall, Director
Subsurface Protection and Remediation Division
National Risk Management Research Laboratory
in
-------
ABSTRACT
Capillary pressure and permeability functions are crucial to the quantitative description of
subsurface flow and transport. Earlier work has demonstrated the feasibility of using the inverse
parameter estimation approach in determining these functions if both capillary pressure and
cumulative drainage are measured during a transient flow experiment. However, to date this
method has been applied to air-water systems only, while ignoring the air phase, thereby assuming
that the air phase has a negligible influence on water flow. In this study, we expanded the inverse
parameter estimation method combined with multi-step outflow data, using a modified Tempe cell
for air-water, oil (Soltrol)-water, and air-oil fluid pairs in a Columbia fine sandy loam and a Lincoln
sand. The commonly applied van Genuchten (VG) - Mualem (M) model of capillary pressure and
permeability functions was used in this study. Wetting fluid and oil non-wetting fluid pressures were
measured in the center of a soil sample simultaneously with cumulative outflow of the wetting fluid
as the initially near-saturated soil core is drained by increasing the non-wetting fluid pressure in a
sequence of pressure increments. Results from the multi-step measurements are used to directly
estimate capillary pressure and wetting fluid permeability functions. Uniqueness and stability
analysis indicated that the inverse model is well posed, and that the multi-step transient outflow
experiment provides sufficient information to successfully apply the parameter estimation approach.
This report was submitted in fulfillment of CR-822204 by the University of California, Davis
under the partial sponsorship of the U.S. environmental Protection Agency. This report covers a
period from October 1993 to September 1996, and work was completed as of September 1996.
-------
CONTENTS
NOTICE ii
FORWARD iii
ABSTRACT iv
FIGURES vi
TABLES vii
ACKNOWLEDGEMENTS viii
Chapter 1 Introduction 1
Chapter 2 Materials and Methods 6
Experiments 6
Numerical modeling 10
Governing equations 10
Boundary and initial conditions 12
Constitutive relationships 14
Numerical solution 15
Optimization 17
Scaling method 19
Chapter 3 Results and Discussion 20
Experimental data 20
Air-water 25
Air-oil 24
Oil-water 25
Simplification to single fluid modeling 26
Estimation of capillary pressure functions 28
Estimation of permeability functions 30
Parameter optimization 34
Chapter 4 Conclusions 44
Appendices 46
A. Optimization Algorithm 46
B. Description of Two-fluid Flow Model 49
References 86
-------
FIGURES
1-1. Methodology of inverse parameter estimation approach 2
1-2. Flowchart of parameter optimization approach 4
2-1. Experimental setup for multi-step outflow experiments with water as wetting and
Soltrol as non-wetting fluid 8
2.2. Schematic representation of boundary and initial conditions 13
2.3.Makeup of conductance G and storage D matrices, and H-vector for 4 nodes 16
3-1. Capillary pressure head (hc) and cumulative outflow (Q) as a function of time for air-
water system of (A) Columbia and (B) Lincoln soil 21
3-2. Capillary pressure head (hc) and cumulative outflow (Q) as a function of time for air-
oil system of (A) Columbia and (B) Lincoln soil 22
3-3. Capillary pressure head (hc) and cumulative outflow (Q) as a function of time for oil-
water system of (A) Columbia and (B) Lincoln soil 23
3-4. Changes of oil and water pressure head with time for oil-water system of
Lincoln soil. 27
3-5. Measured capillary pressurehead (hc) and wetting fluid saturation (Se) data for
(A) Columbia and (B) Lincoln soil 29
3-6. Scaled and fitted capillary pressure head versus effective saturation of wetting fluid
for (A) Columbia and (B) Lincoln soil 32
3-7. Measured relative permeability (symbols) with predicted (solid) and fitted
(dashed) curves using parameters in Table 3-2 33
3-8. Comparison of the measured and optimized hc and Q - values of Lincoln soil: (a)
air-water system, (b) oil-water system, (c) air-oil system 35
3-9. Comparison of the measured and optimized hc and Q - values of Columbia soil:
(a) air-water system, (b) oil-water system, (c) air-oil system 36
3-10. Optimized constitutive functions of Lincoln soil corresponding to the parameter
listed in Table 2.2: (a) Individual capillary pressure head functions, (b) Scaled
capillary pressure head functions, (c) Relative permeability functions 40
3-11. Optimized constitutive functions of Columbia soil corresponding to the
parameter listed in Table 2.2: (a) Individual capillary pressure head functions,
(b) Scaled capillary pressure head functions, (c) Relative permeability functions 41
VI
-------
TABLES
2-1. Physical Properties of Experiment Soils 6
2-2. Physical Properties of Fluids at 20°C 7
2-3. Capillary Pressure Head (hc) and Cumulative Outflow (Q) Data at the End of Each
Applied Pressure Step 7
3-1. Experimental Wetting Fluid Content (9), Capillary Pressure Head (hc), and
Cumulative Outflow (Q) Data for Each of the Three Fluid Pairs 20
3-2. Parameters of Individual Capillary Pressure-Saturation Functions (0S, 9r, a, and n)
and of the Combined van Genuchten-Mualem Model, Fitting the Scaled Capillary
Pressure and Permeability Data (0S, 0r, a, n, and k) 28
3-3. Physical Characteristics of Various Porous Materials 31
3-4 Optimized VG-Mualem Parameters of Lincoln Soil for Different Initial Estimates.. 37
3-5 Optimized VG-Mualem Parameters of Columbia Soil for Different Initial
Estimates 38
3-6. Optimized VG-Mualem Parameters of Lincoln and Columbia Soil 39
3.7 Comparison of a Ratio and Interfacial-Tension Ratio 43
Vll
-------
ACKNOWLEDGEMENTS
We are especially grateful to Dr. John Nieber at the University of Minnesota, for supplying us with
his transient one-dimensional two-phase flow model that we used to develop a specific two-fluid
flow simulator for the inverse parameter estimation. We are also indebted to Dr. N.L. Abbott of the
Chemical Engineering Department of the University of California, Davis, for making his laboratory
available for the interfacial tension measurements.
Vlll
-------
Chapter 1
Introduction
Successful environmental protection and remediation strategies associated with hydrocarbon
contamination of soil and ground water requires modeling of multi-fluid flow and transport in
subsurface soil systems. However, the implementation of such models is often hampered by lack
of sufficient information regarding the capillary pressure-saturation and permeability functions for
the different soil materials. Multiphase fluid flow in soils has been studied across disciplines in soil
science, ground water hydrology and petroleum engineering for decades. Petroleum scientists
have focused predominately on brine-oil-natural gas in mostly coarse-textured soils, whereas
hydrologists study mostly water flow in a broad range of unsaturated soils.
Flow and transport in porous media is controlled by interfacial processes between fluid-fluid and
fluid-solid phases, thereby producing an extremely complex flow field that is dominated by
microscopic heterogeneities and discontinuities. Consequently, the macroscopic capillary
pressure - saturation (energy-mass) and permeability functions, which together characterize fluid
storage and flow properties in unsaturated subsurface soils, are highly nonlinear.
Subsurface properties and flows are discontinuously distributed when viewed at the microscopic
scale. Macroscopic continuity is based on the representative element volume (REV) concept and
is derived by volume-averaging, which led to the classical flow concepts and quantitative analyses
used today. In unsaturated flow, Richards (1931) equation was developed by combining the
empirically obtained closed-form Darcy equation with the mass conservation equation. However,
in doing so, all uncertainties of the microscopic processes were embedded in the macroscopic
constitutive functions; e.g., the capillary pressure and permeability relationships. For decades,
scientists who worked in related areas have been working on prediction or estimation of these
constitutive relationships from microscopic processes or macroscopic observations. Despite
considerable progress by Burdine (1953), Brooks and Corey (1964), Mualem (1976), van
Genuchten (1980), and others, the intricate complexity of pore geometry and microscopic
processes augmented by the severe limitation of observation techniques make prediction of these
relationships difficult. At the same time, the increased efforts in environmental investigations and
numerical simulations demand efficient and accurate methods for determination of soil hydraulic
characteristics. For immiscible multi-fluid flow systems, such information is often lacking.
Many laboratory and field methods exist to determine soil capillary pressure and permeability
functions, which can be categorized as measurement and prediction methods. Direct
measurements, including equilibrium and steady-state experimental methods, are often highly
restricted by constrained initial and boundary conditions, and are time-consuming, or otherwise
inconvenient. This is especially so for the permeability measurements. Moreover, both functions
are measured separately, which can cause inconsistent results. With respect to the prediction
methods, they are mostly based on a simplified conceptual soil pore model, such as by Burdine
(1953) and Mualem (1976) for predicting permeability by using pore-size distribution information
-------
obtained from capillary pressure - saturation data. Alternatively, one applies the similarity
assumption to predict unknown capillary pressure functions or permeability by using known
values of easy-to-measure soil or fluid properties (Leverett, 1941; Miller and Miller, 1956). In
contrast with these traditional methods, the inverse modeling approach for estimating the
constitutive properties is based on the concept of system analysis and is receiving increased
attention. Through the measurement of the system's response of a transient experiment and the
simulation of the experimental system, the inverse modeling approach consistently estimates or
calibrates the system's constitutive properties. Transient experimental methods are inherently
faster and more flexible, and as more powerful computers and simulation models become
available, the estimation of the constitutive functions using the inverse method has become more
attractive.
The study of inverse parameter estimation for determination of retention and permeability
functions started in the 1980's, and was further developed in the 1980's and 1990's (Zachmann et
al., 1981, 1982; Hornung, 1983; Kool and Parker, 1988; Russo et al., 1991; Toorman et al.,
1992; Eching et al., 1993, 1994). Inverse parameter estimation is a "gray-box" technique, when
contrasted with a forward problem (white-box) and inverse problem (black-box). As shown in
Figure 1-1, the grayness is used to describe the exposure degree of the constitutive function
knowledge, and is defined as transparent, opaque, and translucent, for the forward, inverse, and
parameter estimation problems, respectively. The fundamental assumption of parameter
estimation is that the constitutive functions can be described by a parametric model for which the
unknown parameters can be estimated by minimization of deviations between observed and
predicted state variables such as flux or
Forward Problem
^. "White-Box" Technique
Inverse Problem
—». "Black-Box" Technique
Parameter Estimation
—». "Gray-Box" Technique
Pc= Pc(Se)
Kr = Kr(S6)
Out = ?
Out
In
Out
6S, 6r, k, a, n, m, I = ?
Figure 1-1. Methodology of inverse parameter estimation approach.
capillary pressure. To determine whether the inverse problem is at all solvable, it must be
"correctly posed" (Carrera and Neuman, 1986b). Ill-posedness of the inverse problem may result
in no solution, nonuniqueness (more than one solution), or instability (solution is sensitive to small
changes in input data). Uniqueness requires an identifiable parameter set with a solution, which is
-------
sensitive to small changes in the parameters. More detailed discussions on parameter estimation
theory can be found in the papers of Carrera and Neuman (1986a), Kool and Parker (1988), and
Russo et al. (1991). In contrast to linear optimization, there is usually a high uncertainty in
nonlinear parameter estimation (Brooke et al., 1992). Therefore, nonlinear parameter estimation
is often recognized as being "ill-posed." In general, there are three concerns related to the
uncertainty: (i) it may be difficult to find a solution; (ii) if a solution is found, it may not be
unique; and (iii) the solution may be instable or excessively sensitive to experimental data, or
insensitive to one or more parameters. Even though theoretically these are unsolved topics, the
particular problem to be solved may become "well-posed" as appropriate and sufficient
experimental information is obtained.
Previous studies have shown that the experimental method plays an important role in determining
whether the parameter estimation problem is well posed, and indicated that a minimum amount of
data must be collected to characterize the simulated flow process. For example, Gardner's (1956)
one-step transient outflow method may be an ill-posed parameter estimation problem, yielding a
non-unique set of parameters for the constitutive relationships. Van Dam et al. (1994) suggested
that cumulative outflow be measured during multiple outflow steps, so as to yield sufficient
information for determination of a unique parameter set using the inverse method. Including
capillary pressure measurements during the outflow experiment further resulted in improved
parameter sensitivity (Toorman et al., 1992). Eching and Hopmans (1993) concluded that the
measurement of capillary pressure in addition to the multiple outflow measurements provided
adequate information for unique solutions of the inverse parameter estimation problem.
Although the limitations with respect to the experiment or modeling are few, the inverse approach
relies on the availability of a universally applicable nonlinear optimization algorithm. Problems
with the parameter optimization technique generally are associated with the difficulty of defining
an objective function, which will yield unique and convergent solutions.
The inverse parameter estimation of soil capillary pressure and permeability functions includes
three functional parts: (1) a controlled transient flow experiment for which the boundary
conditions and additional flow variables such as capillary pressure and cumulative outflow are
accurately measured; (2) a numerical flow model simulating the transient flow regime of the
experiment and which includes the parametric models that describe the constitutive relationships;
(3) and an optimization algorithm, which estimates the unknown parameters through minimization
of the difference between observed and simulated flow variables in the objective function (Figure
1.2). The quality of the final solution of the parameter estimation problem is dependent on the
quality of each of these three components as well as that of their internal relationships. Moreover,
it is tacitly assumed that the formulated constitutive relationships describe the physical behavior of
the soil in question.
-------
Analysis Structure and Flowchart
Outflow
Experiment
"B.C. &I.C.
"Soil Size ^
-Fluid Properties
x
x
Interfacing x **
Rafa Manag^rrtfent
^^
Outflow
Capillary Pressur
i
/
-/Inf
' file
x
X
1 '
)Ut ^
s/ \.
s ^.
Numerical
Simulation
i
Constitutive
parametric
models
Initial
f"\^ oarameters
N !W
paraneters
Capi
Outflov
illarv P
Lw
ressure
Nonlinear
Optimization
Figure 1-2. Flowchart of parameter optimization approach
Although the inverse parameter estimation approach including experimental and analytical
methods has been developed and increasingly applied in recent years, it has been limited to air-
water systems only. Under the traditional Richards' assumption that the air-phase has a negligible
influence on water flow, air-water systems were treated as one-phase (water) systems. In this
study, the inverse parameter estimation method was expanded to two-fluid flow systems including
air-water, air-oil, and oil-water.
It was the objective of this study to expand the multi-step outflow method to two-fluid flow
systems, and to evaluate how the additional complications and constraints would influence the
well-posedness of the inversion problem. We present results of multi-step outflow experiments for
three two-fluid systems: air-water, oil (Soltrol)-water, and air-oil in two soils. These systems are
typical in immiscible organic contaminant research and important to subsurface flow and transport
studies. In this study we also show how capillary pressure- and permeability-saturation
-------
relationships can be estimated directly from the capillary pressure and drainage data obtained from
a multi-step outflow experiment. Moreover, the scaling concept was tested through scaling of the
individual capillary pressure functions from interfacial tension values, and by normalization of the
effective permeability using the fluid-independent intrinsic permeability value of the investigated
soils. The results are expected to be informative for scientists studying inverse parameter
estimation approaches and multi-fluid flow in the subsurface
-------
Chapter 2
Materials and Methods
Experiments
Multi-step outflow experiments were conducted in a constant temperature (20°C) laboratory
using a modified Tempe cell (Figure 2-1). The cell contained a 7.6-cm high brass soil core with
an outside diameter of 6.4 cm and total soil volume of 216 cm3. Colombia fine sandy loam
collected along the Sacramento river near West Sacramento, California, and Lincoln sand
obtained from the EPA R.S. Kerr Environmental Research Laboratory in Ada, Oklahoma, were
used for the experiments (Table 2-1). Soil was air-dried, sieved through a 2-mm screen, and
uniformly packed. Soil texture and bulk densities for both soils are presented in Table 2-1. For
each separate outflow experiment a freshly air-dried soil was packed to the bulk density values
listed in Table 2-1. One-dimensional transient experiments of a draining wetting fluid replaced by
an invading non-wetting fluid were carried out in three two-fluid systems: (1) air-water, (2) air-oil
(Soltrol 130 for Columbia soil and Soltrol 220 for Lincoln soil), and (3) oil-water. Air is
considered the non-wetting fluid in the air-water and air-oil systems, with oil being the non-
wetting fluid in the oil-water system. Soltrol1 is a mixture of isoalkanes (Cio - CB for Soltrol 130
and CB -Cn for Soltrol 220) and has negligible solubility in water. The relevant physical
properties of the three fluids are presented in Table 2-2. Different types of Soltrol were used to
compare results with earlier measurements. We followed the principles of multi-step outflow for
an air-water system described by Gardner (1956) in which an initially water-saturated soil sample
placed on a fully water-saturated ceramic plate was subjected to a series of step increases in air
pressure, resulting in a cumulative volume of water drainage after each incremental increase of air
pressure. In the multi-step experiments, the rate of cumulative drainage and capillary pressure as
a function of time were measured (Table 2-3).
Table 2-1. Physical Properties of Experiment Soils
Soil type
Columbia
Lincoln
Sand
63.2
88.6
Silt
%
27.5
9.4
Clay
9.3
2.0
Bulk density
g/cm3
1.42
1.69
cm/hr
4.2
23.0
1. Saturated hydraulic conductivity, with water as wetting fluid.
1 Phillips Petroleum Company, Bartlesville, OK
-------
Table 2-2. Physical Properties of Fluids at 20°C.
Interfacial Tension (N/m)
'Soltrol 130
2Soltrol 220
Air-Oil1
0.0239
0.0259
Oil1 Wst^r
0.0259
0.0364
Air Wst^r
0.0681
Viscosity (Ns/m2)
'Soltrol 130
2Soltrol 220
Oil
0.00144
0.00392
Air
Water
0.0000181 0.00100
Density (kg/m3)
1.28
1000
'Soltrol 130
2Soltrol 220
762
803
1. Columbia soil
2. Lincoln soil
Table 2-3. Capillary Pressure Head (hj and Cumulative Outflow (Q) Data at the End of Each Applied Pressure Step.
Applied Pressure (cm)
aw ao ow
aw
Q(ml)
ao
Ow
Measured hc (cm)
aw ao ow
Columbia
60
80
120
200
400
700
20
27
40
67
133
233
40
53
80
133
266
466
10.0
16.0
35.5
52.4
60.2
66.1
7.1
14.0
25.3
59.0
68.0
69.1
19.0
29.4
48.4
55.6
63.2
67.5
64.4
83.8
120
198
320
682
24.2
30.2
41.2
62.4
109
163
43.6
54.4
76.6
116
218
345
Lincoln
40
60
80
100
150
200
400
10
15
20
27
33
67
20
32
42
50
73
100
199
13.5
31.0
38.5
45.0
52.5
54.5
57.5
9.0
16.0
26.5
43.5
51.5
57.5
8.8
27.5
40.8
45.2
53.0
55.0
57.0
43.3
60.6
79.8
99.3
148
194
261
13.1
15.8
20.1
25.3
29.0
38.8
23.7
32.3
40.1
47.1
66.5
80.0
90.7
1.Non-wetting fluid entry pressure with air as non-wetting fluid
2.Non-wetting fluid entry pressure with Soltrol as non-wetting fluid
-------
Pnw
I
T2(w)
7.6cm
0.74
T3(nw)
Wetting fluid
Figure 2-1. Experimental setup for multi-step outflow experiments with water as wetting andSoltrol as non-wetting fluid. T denotes a pressure
transducer, and Tl and T2 are the tensiometers for the wetting and non-wetting fluid, respectively.
Using a noninvasive x-ray computed tomography technique, Hopmans et al. (1992) confirmed
previous studies indicating that fluid flow can be properly described only if both wetting and non-
wetting fluids are continuous. To obtain continuous wetting and non-wetting fluids at the onset
of the experiment, the soil core was placed on a high flow rate, 1-bar ceramic plate, and saturated
by the wetting fluid from the bottom upwards. The ceramic plate was 0.74 cm thick with a water-
saturated hydraulic conductivity of 0.048 cm/h. The ceramic plate accepted Soltrol as a wetting
fluid in air-oil systems, just as easy as it absorbs water in air-water systems. Therefore, the
ceramic plate could be used for both water and Soltrol as the wetting fluid, without the need for
additional treatment. The soil sample, saturated with the wetting fluid, was subsequently drained
by applying a suction to the ceramic plate, slightly higher than the non-wetting fluid entry pressure
of the investigated soil for that particular fluid pair. The resulting initial condition was hydraulic
equilibrium for both fluids, after a static capillary pressure profile was achieved and both fluids
were continuous. As positive pressure increments of the non-wetting fluid were applied at the
-------
surface of the soil, cumulative outflow of wetting fluid collected in a burette was monitored as a
function of time by measurement of the fluid pressure, using a 1-psi pressure transducer2
connected to the bottom of the burette (Figure 2-1). When air was the non-wetting fluid, air
pressure was applied directly through a hole at the top of the Tempe cell. However, when oil was
the non-wetting fluid, constant oil pressure for a given pressure increment was maintained using a
Marriotte siphon controlled bottle filled with oil connected to pressurized air at one side with the
other tube connected to the top of the flow cell and completely filled with oil. By stepwise
adjustment of the air pressure, a stepwise constant oil pressure at the upper boundary of the flow
cell was attained.
Ceramic tensiometers, connected to transducers, were installed vertically with the tip of the
tensiometer placed at the 3.8-cm depth below the soil surface to measure pressure changes of the
wetting and non-wetting fluids during drainage of the wetting fluid. Only one tensiometer was
needed for the measurement of the wetting fluid pressure in the air-water and air-oil system, since
the air pressure is equal to the applied air pressure across the sample at all times. If oil was the
non-wetting fluid, both a hydrophilic and a hydrophobic tensiometer were used to monitor
pressure changes of the water and oil fluid, respectively. Hydrophobic tensiometers were
obtained using the treatment methods described by Lenhard and Parker (1987) and Busby et al.
(1995). Since Soltrol is corrosive to the transducer membrane, oil pressures were measured by
filling the transducers with water, after which they were connected to the Soltrol-filled Teflon
tubing. All transducers were multiplexed and connected to a datalogger for automatic data
acquisition of pressures and cumulative outflow during transient drainage at a measurement
frequency of 12 readings per minute. Tensiometers were rigid and they were completely filled
with either the wetting (water) or non-wetting (Soltrol) fluid, thereby allowing their use as a fast-
response tensiometer. The immediate response of the transducer to changing fluid pressures was
independently determined, and a response time of less than 10 seconds was adequate for the
transient type of experiments described here.
Applied air pressures in the air-water system were 60, 80, 120, 200, 400, and 700 cm above
atmospheric pressure for the Columbia soil. These pressure steps were chosen based on previous
work by Eching et al. (1994). Selected air pressure steps for the Lincoln soil were 40, 60, 80,
100, 150, 200 and 400 cm above atmospheric pressure. These steps were smaller than that for
Columbia soil because of the Lincoln soil's coarser texture. The pressure steps for the other fluid
pairs (Table 2-3) were determined using these pressures and the scaled relationships of the
interfacial tension between the tested fluid pair and air-water (Table 2-2) in an effort to apply
pressure steps that yielded approximately equal drainage volumes for each wetting fluid within
each pressure step for an investigated soil. A pressure increment lasted from 5 to 36 hours, and
varied according to the pressure of the wetting fluid. After the wetting fluid pressure head
(defined by hw = Pw/pn2og, where Pw is the wetting fluid pressure, and pH2o and g are the water
density and gravitational acceleration constant, respectively) remained approximately constant
(i.e., variations of less than 1 cm), the non-wetting fluid pressure was incrementally increased for
the next drainage step. Final wetting fluid saturation was determined by oven-drying the soil
: MICRO SWITCH, Freeport, IL 61032.
-------
following the last pressure step for the air-water and air-oil systems. Since the oven-drying
method can not be used for determination of the final wetting fluid saturation for an oil-water
system, the final wetting fluid saturation value for oil-water systems was determined from the
average value of the air- water and air-oil systems. The saturated wetting fluid content and initial
saturation were calculated based on the final degree of saturation and cumulative outflow.
Detailed information on procedures for the multi-step experiment as applied to air-water systems
can be found in Eching and Hopmans (1993).
The original outflow experiments for an air-water system assumed a constant air pressure in the
draining soil, since air is continuous and its viscosity and density values are relatively low as
compared to water. As such, a positive air pressure in the soil sample is assumed to be equivalent
to a negative water pressure applied at the lower soil boundary, with the air pressure in the soil at
atmospheric pressure (Kool et al. 1985a). This allows use of a single-fluid Richards' equation
model for simulation of the draining soil core. If the non-wetting fluid is oil, however, the
assumption of a constant non-wetting fluid pressure is not necessarily valid, since the viscosity
and density of the oil fluid are similar to that of water (Table 2-2). Thus, a hydrophobic
tensiometer was used to monitor the oil pressure during water drainage of the oil-water system.
Values of the interfacial tension between air and water are documented in handbooks and
textbooks. But reported interfacial tension values of air-Soltrol and Soltrol-water vary widely.
Moreover, since different types of Soltrol (Soltrol 130, Soltrol 170, and Soltrol 220) have been
presented for flow and transport investigations in porous media, it was difficult to find exact
values from the literature. Therefore, we independently measured the interfacial tensions of air-
water and air-oil (Soltrol 130 and Soltrol 220) using the plate method (Adamson, 1990), while the
interfacial tension between oil (Soltrol 130) and water was measured by the ring method
(Adamson, 1990). The drained pore water at the completion of the air- water experiments was
used for measurement of air-water surface tension, thereby including possible interactions of the
soil material with the draining fluid. The interfacial tension between Soltrol 220 and water was
taken from the work of Schroth et al. (1995), who concluded that the interfacial tension for an
oil-water interface is time-dependent due to oil contamination of water at oil-water interfaces.
Numerical modeling of two-fluid phase flow
Governing Equation
Under the assumption that both fluids and the porous medium are incompressible, the most
general form of the two-fluid flow equations without source-sink terms is described by the two-
fluid volume-averaged momentum and continuity equations (Whitaker, 1986):
qw = -PW + Pwg + w,nw -qnw
10
-------
+Pnwg + nw>w -q
M'nw
In Eqs. (2-1), the subscripts w and nw denote the wetting and non-wetting fluids, respectively; P;
(i = w, nw) denotes pressure (N/m2); S; (i = w, nw) is the degree of fluid saturation relative to the
porosity ((>; q; is flux density vector (m/s); (J,; (i = w, nw) denotes fluid dynamic viscosity (Ns/m2);
k; (i = w, nw) is the effective permeability tensor (m2) = kr; k, where k is the intrinsic permeability
(m2) and kri = kri(S;) is the relative permeability. In Eq. (2-1), the cross term k;j denotes the
viscous drag tensor (Whitaker, 1994) representing the influence of the viscous drag that exists
between the flowing wetting and non-wetting fluids. The cross term describes the coupling
between the nw- and w- fluids, which appears to be highly dependent on the viscosity ratio of
both fluids (Whitaker, 1986 and 1994), and the magnitude of the interfacial areas between the two
fluids. Since these cross terms are either undetectable or unimportant or both for flow in complex
soil systems, the current practice is that they are of secondary importance. Moreover, as stated by
Whitaker (1986), in flow systems with air as the non-wetting fluid, its viscosity is small relative to
that of the wetting fluid, thereby justifying the insignificance of the coupling parameters. Also,
Bentsen (1994) presented experimental evidence that removing these terms introduced little error.
Eqs. [2-1] also hold if air is the non-wetting fluid as in multi-step outflow experiments, if it is
assumed that incremental changes in applied air pressure occur instantaneously across the soil
sample, and that variations in soil air pressure between pressure increments have a negligible
influence on the air density. Otherwise, the continuity equation of the air phase should include a
pressure-dependent air density (Celia and Binning, 1992). Consequently, for one-dimensional
vertical flow systems with z denoting vertical position, Eq. (2-1) is simplified to:
(2-2a)
dt dz
k,,, SP
(2-2b)
k 8P
g) (2-2c)
di dz
For an incompressible porous medium, Eq. (2-2) is supplemented by:
11
-------
Sw + Snw = 1 (2-3)
Substituting (2-3) into (2-2a), and using the definitions of capillary pressure (Pc = Pnw - Pw) and
fluid capacity (C = - (|)dSw/dPc), one obtains the following governing equations,
P ) d k dP
)_=o^L_^ +
az |j,w dz
/0 ,x
(2-5>
which can be solved simultaneously for the unknown pressures Pw and Pnw, and thus for Pc.
In addition, since only non-wetting fluid enters the soil core at the top, and only wetting fluid
leaves the core through the bottom, volume balance considerations across the cell at any time (t)
requires that:
(2-6)
Boundary and Initial Conditions
The boundary conditions (EC's) and initial conditions (IC's) are determined by the experimental
conditions of the transient multi-step outflow experiment during which the wetting fluid is drained
from an initially slightly-unsaturated soil core. Figure 2.2 shows an overview of all IC's and EC's.
At the upper boundary of soil core, the wetting fluid density flux is zero, and the non-wetting fluid
pressure is prescribed by the imposed series of multi-step pressures, described by the stepwise
function Pj(Tj):
qw(ztop,t) = 0 (2-7a)
POT(ztop,t) = Pro(Tj), j = l,2,...,M (2-7b)
In Eq. (2-7b), M is the total number of pressure steps used in the multi-step experiment, and
Pnw(Tj) is the pressure value applied to the non-wetting fluid at ztop during the time period Tj.
Using this notation, TI is the time period when pressure step one is applied, T2 is the time period
when the pressure step two is applied, and so on. At the lower boundary of the flow system, the
flux of the non-wetting fluid is zero, since the entry value of the ceramic plate for the non-wetting
fluid is higher than any of the imposed pressures. The wetting fluid pressure at Zbottom is the
12
-------
prescribed pressure condition determined by the height of the wetting fluid in the burette, which is
a function of time, as the wetting fluid drains into the burette:
Top -
B.C.
= Pnw(Tj), j = 1, .... M
Pnw - Pnw(to, z)
Pw = Pw(to, z)
qnw = 0, Pw = PW (Zbottom, tj), i = 1, .... N
Bottom
- B.C.
Figure 2.2 Schematic representation of boundary and initial conditions.
13
-------
qnw(zbottom,t) = 0 (2-8a)
pw (zbottom, t) = Pwghw|outflow (t) (2-8b)
In Eq. (2-8b), pw is the density of the wetting fluid and hw|outflow (t) is the height of the wetting fluid
in the receiving burette above the bottom of the ceramic plate (Figure 2.1). Because limited
observations in time were used in the numerical modeling, the lower boundary-pressure
measurements are a discrete function of time t;, i = 1, 2, ..., N.
At time zero, both fluids are in static equilibrium. Because the soil is slightly unsaturated with
respect to the wetting fluid, the initial conditions are described by the hydrostatic pressure of each
fluid:
Pw (z, t0) = Pw (zbottom, t0) - pwgz (2-9a)
Pnw (z, t0) = Pnw (ztop, t0) + pnwg(ztop - z) (2-9b)
where z is assumed positive upwards and z = 0 at the bottom of the ceramic plate. Pw(zbottom, t0)
and Pnw(ztop, t0) denote the boundary pressure values at the start of the transient outflow
experiment.
Constitutive Relationships
The governing Eqs. (2-4) and (2-5) require apriori knowledge of the capillary pressure
function hc(Sw) and permeability function k;(Sw) with i = w, nw, which are defined by functional
parametric models:
Sw=Sw(hc,b) (2-10a)
kw=kw(Sw,b) (2-10b)
knw=knw(Sw,b) (2-10c)
where b denotes the vector containing the parameters of the assumed functions. In this study, we
used the van Genuchten model (1980) to characterize the capillary pressure function, which is
used with Mualem's model (Mualem, 1976, Parker et al., 1987; Luckner et al., 1989) to describe
the permeability functions:
14
-------
n\m (2-11 a)
_ kw _ i i „ 2
7 nw /i O \ ^ n O i2w? /o 11 \
"Ynw = = (1 ~ ^ew) U ~ ^ew"1 J (2-1 1C)
where, kr is the relative permeability and k denotes the intrinsic permeability (L2); Sew is the
effective saturation of the wetting fluid with Sew = (Sw-SW;r)/(l-SW;r), where Sw = 9W/(|) is the
saturation of wetting fluid, and Sw>r is the residual saturation of the wetting fluid; a and n are
unknown parameters which are inversely proportional to the non-wetting fluid entry value and the
width of pore-size distribution, respectively, and m was assumed to equal to m = 1-1/n (van
Genuchten, 1980); the parameter 1 is related to the tortuosity of the soil and is here assumed
equal to 0.5 for both the wetting and non-wetting fluid.
Numerical Solution
The governing Eqs. (2-4) and (2-5), the boundary and initial conditions (2-7, 2-8, and 2-9), and
the constitutive relationships in Eq.(2-ll) combined make up the mathematical model of the
experimental system. The mathematical model has no analytical solution available because of the
nonlinearity of the constitutive functions. Therefore, a numerical model was adapted to simulate
the two-fluid flow regime.
The adapted two-phase numerical model in this study was developed from a two-phase model of
Dr. John Nieber at the University of Minnesota (personal communication). We used the same
numerical scheme, which includes a modified Picard linearization algorithm of the mixed-form
governing equation and a lumped finite element approximation (Celia et al., 1990 and 1992).
Celia et al. (1990) showed that the numerical solution based on the mixed form equation was
inherently mass conservative and that the lumped finite element approximation eliminated
oscillations.
Briefly, by using total head H instead of pressure (Hw and Hnw for wetting and non-wetting
PI Pi
phases: H= = — + —~z), the numerical approximation of governing Eqs. (2-4) and (2-
PH20§ PH20
5) becomes:
,,-H.) A(^^}
Tt—=—^— <2'12)
15
-------
and
-C
(Hnw-Hw)
> x
nw __ nw
Ziz
(2-13)
Applying the modified Picard and lumped finite element approximation, the finite element matrix
equation can be written as
(At
D)Httl=Ziq-DHt
(2-14)
where, G is conductance matrix derived from a combination of individual permeability terms, D is
the storage matrix with elements accounting for the capacity term C, and A6 is a vector
describing fluid-content changes for a time increase At. Figure 2.3 demonstrates the makeup of G
and D, from elemental matrices (sub-matrices El, E2, and E3) and nodes (Nl, N2, N3, and N4).
Figure 2.3 also demonstrates that each node has
N1
N2
N3
N4
* - *
* - *
* -
*
*
*
*
*
*
- *
*
. *
*
- *
* . *
* *
HI w
H-l n
H2iw
H2,n
Hs.w
H 3 n
H 4 w
H4,n
H
Figure 2.3 Makeup of conductance G and storage D matrices, andH- vector for 4 nodes.
N = node, E = element. Dashes indicate zero-value entries.
both wetting and non-wetting fluid contributions in all arrays and vectors of Eqs. (2-12) and (2-
13). The first subscript of total head H vector represents the index of the finite element nodes.
The entries out of the diagonal line in Matrix D represent the coupling between the two-fluid
pressures of each node. From the makeup of matrices in Figure 2.3, it is clear that the left-hand
side of Eq.(2-14) consists of a matrix with a total bandwidth of five. The quintic-diagonal
algorithm (Vemuri and Korplus, 1981) was used to solve Eq.(2-14). The air-water flow option of
this two-fluid flow model was tested by comparing simulations with the model used by Eching
and Hopmans (1993).
16
-------
Optimization
The inverse parameter estimation is cast as a nonlinear optimization problem; i.e., a vector b
containing the unknown parameters of the constitutive relationships in Eq.(2-ll) is estimated by
minimizing an objective function O(b), containing deviations between observed and predicted
system response variables. In general, an optimization procedure includes a formulation of the
objective function, a solution algorithm, and an analysis of convergence and uncertainty.
Theoretically, the objective function can be formulated from a maximum likelihood (ML)
consideration (Carrera and Neuman, 1986a; Kool and Parker, 1988). For a given predictive
model, which is assumed to be true with no error, ML yields parameter estimates with non zero
values of the objective function being attributed to measurement errors.
Using the notation v which is the vector containing the elements of measured flow variables
(capillary pressure, outflow) with the subscript m and s denoting measured and simulated values,
respectively, the observation error is equal to e = vm - vs, with the assumption that model errors
are zero. When observation errors (n is the total number of observations) are assumed to describe
multivariable normal distribution with zero-mean and covariance matrix
V = E(vm - vs)(vm - vs)T, the maximum likelihood estimation (Bard, 1974) formulates the
objective function O(b) as:
O(b) = -ln27r + -ln(detV) + -eTV-1e (2-15)
and leads to the least squares (LS) problem:
O(b) = eTV 'e (2-16)
Thus, for a general covariance matrix V, Eq. (2-16) is a general least-squares problem (GLS),
where the inverse of the error covariance matrix denotes the weighting matrix, and includes
measurement accuracies and their correlations. If these errors are independent, but the variance
varies among observation type, Visa diagonal matrix leading to a weighted least squares problem
(WLS). For the case that measurement errors are assumed independent with constant variances,
V is an identity matrix of size n x n, and Eq. (2-16) reduces to an ordinary least squares problem
(OLS). In this study, the assumption is made that the errors e are independent and that their
variances are proportional to the magnitudes of the mean values of particular measurements (Kool
and Parker, 1988). Weighting factors are chosen to be inversely proportional to their mean
measured values of cumulative drainage (Q), capillary pressure head (hc) and initial water content
9(hc, to), yielding a WLS problem of the form:
17
-------
(2-17)
+W02^{cdk[9m(hck,tk)-i
k=l
or using W = V"1,
0(b) = [vm - vs(b)]TW[vm - vs(b)] = eTWe
(2-18)
In Eq. (2-17), N, M, and L denote the number of observations of Q, Pc, and 9, respectively and
O(b) is normalized by the weighting factors:
WQ=1
(2-19a)
(2-19b)
Wa =
(2-19c)
where L = 1 with 0m(hc,0, t0) corresponding to the initial volumetric water content value after
exceeding the non-wetting fluid entry value. In addition, co;,
-------
We applied the Levenberg-Marquardt optimization algorithm (More, 1977) to minimize (2-20).
The details of Levenberg-Marquardt method and the convergence and uncertainty analysis are
included in the Appendix. In this study, an existing optimization program (Eching and Hopmans,
1993) originally developed by Kool et al. (1985a) was modified to interface with the two-fluid
flow model.
Scaling Method
The traditional method of estimating the capillary pressure function for one fluid pair from
another is based on the scaling method as first introduced by Leverett (1941). The theoretical
basis stems from the similarity theory, when identical soils with different fluid-pairs are considered
to be similar systems. The basic assumptions include that the solid matrix is rigid with negligible
solid-fluid interactions, fluids are held in the porous matrix by capillary forces only, and air-water
and oil-water interfaces act independently. Capillary pressure as a function of fluid saturation is
determined by interfacial properties such as interfacial tension, contact angle and interfacial
curvature. Whereas the interfacial tension and contact angle are dependent on the particular solid
and fluid materials, the interfacial curvature is dependent on the pore geometry of the soil matrix.
Consequently, Leverett (1941) introduced the dimensionless Leverett's function J(Sew) to
describe the similarity relationship between soil systems 1 and 2 by,
pci k, i PC 2 k, i
where a and k denote interfacial tension and intrinsic permeability, respectively. In Eq. (2-21),
(kA|))"1/2 (length unit, L) denotes the microscopic length for which the size depends on the pore
geometry of soil medium only. Thus, for the same soil matrix but different fluid pairs 1 and 2, Eq.
(2-21) reduces to:
) = L = 1 (2-22)
or, he2(Sm) = hel(Sm) (2-23)
<*\
Eq. (2-23) states that from a known soil's capillary pressure head-saturation relationship for a
particular fluid-pair (e.g., air- water), the soil's capillary pressure function for another fluid-pair
can be predicted according to the corresponding interfacial tension values. Additional assumptions
for proper application of Eq. (2-23) are: (i) the soil is completely water-wet; (ii) the oil fluid is
present between water and air, and (iii) no fluid trapping occurs. Also, it is assumed that the
contact angle is independent of fluid type, which may not be correct (Demond and Roberts,
1991). The optimized capillary pressure PC(SW) function from the proposed inverse parameter
estimation of the air-water, air-oil, and oil-water systems are scaled to compare the scaling
relationships, and to provide a means to test the optimized solutions.
19
-------
Chapter 3
Results and Discussion
Experimental data
Measured values of capillary pressure and cumulative outflow at the end of each pressure step,
together with the applied non-wetting fluid pressure are listed in Table 2-3. Figures 3-1,3-2, and
3-3 present all the collected cumulative outflow and capillary pressure data as a function of time
for air-water, air-oil, and oil-water, respectively, for both soils. Table 3-1 summarizes
experimental observations of saturated wetting fluid content, initial and final capillary pressures,
and cumulative outflow for the two soils.
Table 3-1. Experimental Wetting Fluid Content (9), Capillary Pressure Head (h^ and Cumulative Outflow (Q) Data for
Each Fluid Pair of Investigated Soils.
Columbia Soil Lincoln Soil
Air-Water Air-Oil Oil-Water Air-Water Air-Oil Oil-Water
0,1 (m3/m3)
Initial hc2 (cm)
9f3 (m3/m3)
Final hc (cm)
Total O fmn
0.45
23.0
0.14
682
66.1
0.43
10.0
0.11
163
69.4
0.444
18.0
0.13
345
67.5
0.32
16.0
0.066
261
56.7
0.33
7.2
0.059
39
57.5
0.334
10.0
0.059
91
57.0
1 Saturated wetting fluid content
2 Under suction
3 Final wetting fluid content
4 Estimated value
As shown by Hopmans et al. (1992) for air-water systems, the initial capillary pressure head (hc)
must exceed the non-wetting fluid entry pressure, so as to achieve non-wetting fluid continuity at
the onset of the transient drainage experiment. The corresponding wetting fluid saturation is
estimated from the saturated (9S) and cumulative wetting fluid drainage at static equilibrium when
the initial capillary pressure head value has been attained, thereby yielding the first point of the
capillary pressure curve. The final wetting fluid content (9f) was determined from oven-drying
upon completion of each experiment. The saturated wetting fluid content (0S) was calculated
from total cumulative drainage (total Q) and the final wetting fluid content.
20
-------
dd
1000
2000
3000 4000 5000
Time (minutes)
6000
200 T
T60
O
0
1000 2000 3000 4000 5000 6000 7000 8000 9000
Time (minutes)
Fig. 3-1. Capillary pressure head (h,,) and cumulative outflow (Q) as a function with time for air-water
system of (A) Columbia and (B) Lincoln soil.
21
-------
dd
120 T
1000
2000
3000 4000
Time (minutes)
5000
6000
O
50 T
T60
2000
4000
Time (min)
6000
O
Fig. 3-2. Capillary pressure head (hc) and cumulative outflow (Q) as a function with time for air-oil system of (A)
Columbia and (B) Lincoln soil.
22
-------
dd
120 T
1000
2000
3000 4000
Time (minutes)
5000
6000
50 T
T60
2000
4000
Time (min)
6000
8000
Fig. 3-3. Capillary pressure head (hc) and cumulative outflow (Q) as a function witii time for oil-water system of
(A) Columbia and (B) Lincoln soil.
23
-------
However, for the oil-water system, the 9s-value was assumed to be the average of measured 0S-
values for the air-water and air-oil systems for each soil.
Since the Lincoln soil is coarser than the Columbia soil (Table 2-1), the Lincoln soil has smaller
initial capillary pressure head values than the Columbia soil. Due to the scaling of the non-wetting
fluid pressure at the top boundary of the soil sample, total outflow values are approximately equal
for the three fluid pairs for each soil.
Air-water system
The measured cumulative outflow and capillary pressure head of the wetting fluid as a function of
time for the Columbia and Lincoln soils are presented in Figures 3-la and 3-lb, respectively. The
capillary pressure is computed from the difference between the non-wetting and wetting fluid
pressures, and is expressed in capillary pressure head, hc (cm of water). For the air-water
systems, the soil air pressure was assumed to be equal to the applied air pressure everywhere
within the sample, whereas the water pressure was measured by the hydrophilic tensiometer at the
center of the soil core. Outflow rates are relatively high at the beginning of each pressure
increment and decrease toward zero, as the cumulative outflow (Q) approaches a constant value.
After the drainage rate is reduced to near zero, the non-wetting fluid pressure was incrementally
increased for the next drainage step. The initial high flow rates are the result of the larger capillary
pressure gradients when air pressure is increased, and the relatively high effective permeability.
As time passes, the capillary pressure head gradient is reduced and the outflow rate decreases
toward zero during each pressure step. The capillary pressure head changes rapidly along with
the cumulative outflow data for the first few pressure increments at relatively high wetting fluid
saturations. However, as the capillary pressure is further increased, drainage rates decrease
because of decreasing wetting fluid permeability as the soil desaturates, thereby increasing the
time period required to achieve hydraulic equilibrium. As is usually done for estimation of
capillary pressure-saturation functions under equilibrium flow conditions, we used simultaneously
measured capillary pressure and drainage volume data towards the end of each pressure increment
to estimate this function. The data in Table 2-3 also clearly demonstrates that as the applied
pressure increases the difference between the applied pressure head and measured capillary
pressure head at the end of each pressure increment becomes larger. Note that at true hydraulic
equilibrium the measured capillary pressure in the center of the soil core should be equal to the
applied air pressure. This tendency of not reaching equilibrium within a reasonable measurement
time period had been determined earlier by Eching and Hopmans (1993) for the Oso Flaco sand,
and is caused by the low wetting fluid permeability at low saturations, especially for coarse-
textured soils. We should also point out that the estimation of the soil's capillary pressure and
permeability data does not require equilibrium soil-water conditions, but assumes a constant
capillary pressure gradient in the soil sample. This assumption is satisfied if the soil samples
approach hydraulic equilibrium.
24
-------
Air-oil system
As in the air-water system, air pressure across the vertical profile of the sample is equal to the
applied air pressure for each pressure step in the air-oil system. This system differs from the air-
water system only by the physical properties of the wetting fluid. Based on the scaling theory of
Leverett (1941), the capillary pressure-saturation curve for an air-oil system can be determined
from that for an air-water system using the scaling relationship (see also Eq. 2-23):
hc(ao) = (— )hc(aw (3-1)
c(ao)
aw
where hc denotes the capillary pressure head (cm of water), a is the interfacial tension (N/m), and
the subscripts ao and aw indicate air-oil and air-water systems, respectively. Using this relation,
the range of measured capillary pressure heads in the air-oil system is approximately one third of
the air- water systems (Table 2-2) for the same range in degree of wetting fluid saturation.
Though pressure increments for the air-oil system experiments were adjusted with respect to the
ratio of interfacial tensions, the drainage curve differs from that of the air-water system, because
times at which the applied air pressure was changed were different for the two systems.
However, the measured capillary pressure and cumulative oil outflow versus time data in the air-
oil systems for the Columbia and Lincoln soils as shown in Figures 3-2a and 3-2b, respectively,
are similar to those for the air-water systems in Figures 3-la and 3-lb, with the exception of less
distinct plateau values for a few pressure steps in the Lincoln soil.
Oil-water system
Shown in Figures 3-3a and 3-3b are the outflow and capillary pressure head changes with time for
the oil-water system in the Columbia and Lincoln soils, respectively. Again, the shapes of the
curves in Figures 4a and 4b are similar to those of Figures 3-la and 3-lb, while the range of
capillary pressure heads for the oil-water system is between those of air-water and air-oil systems
as determined by the adjusted applied air pressures.
Although we assumed a constant non-wetting fluid pressure profile for the air-water and air-oil
systems, such a profile was not expected to be the case for oil-water systems, since the viscosity
and density values of the non-wetting fluid (Soltrol) are similar to that of water (Table 2-2). As an
example, Figure 3-4 shows the pressure head changes of the non-wetting (h0;i = Poii/pmog) and
wetting fluids (hwater = Pwater/Pmog) with time for the Lincoln soil. As expected, the soil-water
pressure first increases in response to the increased oil pressure and subsequently decreases with
time as the soil water drains, approaching near zero (wetting fluid pressure below ceramic plate)
as hydraulic equilibrium is established. However, as the soil desaturates with respect to the
wetting fluid (water) at the higher applied oil pressures, the wetting fluid permeability becomes
limiting and equilibrium is not achieved within the measured time period. With an increase of the
25
-------
oil pressure at the upper boundary of the soil core, the oil pressure in the center of the sample
increases immediately to a value equal to the incremented pressure and thereafter remains
constant with time. We postulate that the constant oil pressure in the soil sample is a consequence
of the imposed boundary conditions. In our experiments, the oil pressure is equal to the applied oil
pressure corrected for static oil pressure at any time for both soils, thereby simplifying the
physical description of wetting fluid drainage in the oil-water system.
Simplification to single fluid flow modeling
The one-dimensional Darcy flow equation combined with the continuity equation is used to
describe single-fluid transient flow processes in soils. For a homogeneous soil, the flow and
continuity equations for the wetting fluid (as indicated by the subscript w) are (see also Eqs. 2-2a
and 2-2b):
(3-2)
(3-3)
where qw is Darcy's flux (L/T), kw is the effective permeability (L2) of the wetting fluid, Pw is the
wetting fluid pressure (M/L T2) , pw is the density of the wetting fluid (M/L3), |j,w is viscosity of
the wetting fluid ( M/T L) , g is the gravitational acceleration constant (L/T2), (() is porosity, Sw =
9/9s (-) is degree of saturation of the wetting fluid, t is time (T), and z is the vertical coordinate
(L, positive upwards). The effective permeability is defined as:
(3-4)
where k is the intrinsic permeability (L2) and kr is the relative permeability of the wetting fluid,
which is a function of degree of saturation and is related to the unsaturated hydraulic conductivity
K(SW) by
K(S ) = g2o (3-5)
26
-------
200 T
-50
2000 4000 6000
Time (minutes)
8000
Figure 3-4. Change of oil and water pressure head as a function of time for oil-water system of Lincoln soil
where |J,W is the viscosity of the wetting fluid. Combining Eqs. (3-2) and (3-3) yields:
di
(3-6)
where Pc=Pnw-Pw is the capillary pressure. While Eq. (3-6) has two unknowns, Pnw and Pw, we
already noted that changes of the non-wetting fluid pressure with time were negligible in our
experiments, hence Eq. (3-6) can be rewritten as:
f~\ w ___[
dt dz
d
(3-7)
where C = (|)(dSw/dPc) is the slope of the capillary pressure-saturation curve. Hence, it seems to
be adequate to apply a single-liquid flow model to simulate pressure changes of the wetting fluid
in an oil-water system, much the same as in predicting water pressure changes in an air-water
system (Kool et al., 1985b). Thus, the assumption of the time-independent non-wetting fluid
pressure, equal to the applied non-wetting fluid pressure and augmented with the hydrostatic oil
pressure in the soil core, simplifies the wetting-fluid phase flow to the original Richards equation.
27
-------
Estimation of capillary pressure function
The outflow and capillary pressure data in Figures 3-1, 3-2, and 3-3 show that both measured
capillary pressure and cumulative outflow change rapidly at the beginning of each pressure step.
In most experiments, drainage flow rate approaches zero at the end of each pressure step when
the capillary pressure head in the center of the soil core reaches a constant value. At the end of
each pressure increment then, the capillary pressure profile can be assumed to vary linearly with
vertical position in the 7.6-cm tall soil core (Eching et al. 1994). Therefore, the measured outflow
and capillary pressure data at the end of each pressure step can be used to estimate the capillary
pressure-saturation function. The volumetric wetting fluid content (9) is an average value of the
soil core, just prior to increasing the non-wetting fluid pressure, calculated from cumulative
outflow and initial wetting fluid content. The capillary pressure measured simultaneously in the
center of the soil sample is thus assumed to correspond to the sample-average capillary pressure.
The capillary pressure-saturation function in Eq. (2-11 a) of van Genuchten equation (1980) was
fit to the experimental data shown in Fig. 3-5. Through curve fitting, using nonlinear least
squares optimization, provided in the spreadsheet program EXCEL® (Wraith and Or, 1997), the
parameters of the van Genuchten capillary pressure functions (Table 3-2) were obtained for each
fluid pair. Figure 3-5 shows the measured (data points) and fitted capillary pressure-saturation
functions (lines) for the three fluid pairs and both soils.
Table 3-2. Parameters of Individual Capillary Pressure-Saturation Functions (9S, 9r, a, and n) and of the Combined Van
Genuchten-Mualem Model, Fitting the Scaled Capillary Pressure and Permeability Data (9S, Q,, a, n, and k).
Air-
Water
Columbia
Air-
Oil
Oil-
Water
Van Genuchten-
Mualem model
Air-
Water
Lincoln
Air-
Oil
Oil-
Water
Van Genuchten-
Mualem model
a (cm"1)
k (cm )
0.44
0.120
0.009
3.197
0.44
0.112
0.024
4.423
0.43
0.132
0.02
3.352
0.44
0.112
0.009
2.873
2.0 x 10-9
0.32
0.06
0.019
3.686
0.31
0.001
0.051
4.180
0.32
0.056
0.033
4.633
0.32
0.062
0.019
5.712
7.0 x 10'!
28
-------
•s
^» s
•- u
es ^^
800
600
400
200
0
A. Columbia
• air-water
• air-oil
A oil-water
0.000
0.200
0.400
0.600
0.800
1.001
Se of wetting fluid
300
0
0.000
0.200 0.400 0.600
Se of wetting fluid
0.800
1.000
Figure 3-5. Measured capillary pressure head (h^ and wetting fluid saturation (Se) data for (A) Columbia and (B) Lincoln
soil
29
-------
The 9r-values reported in Table 3-2 are obtained from fitting capillary pressure data to the van
Genuchten Eq. [3-8], and were not independently measured as done by Demond and Roberts
(1991).
Using the interfacial tension ratios, the estimated capillary pressure-saturation data of air-oil and
oil-water systems for each soil were scaled relative to the air-water capillary pressure data, and
these are plotted versus effective saturation in Fig. 3-6. Demond and Roberts (1991)
demonstrated that the coalescing of capillary pressure data using the interfacial tension-scaling
concept was more successful if the capillary pressure is plotted against effective saturation of the
wetting fluid, rather than simple saturation. Although not presented, we found little difference in
scaling results when hc was plotted against wetting fluid saturation. Nevertheless, the presented
scaling results in Fig. 3-6 demonstrate excellent agreement between scaled and measured
relationships for both soils, with spreading between scaled and measured capillary pressure data
increasing at lower values of effective saturation. Similar findings were reported by Demond and
Roberts (1991), who attributed discrepancies to the uncertainty of residual saturation values and
the possible need to include the contact angle in the original Leverett scaling relationship. The
fitted curve in Figure 3-6 was obtained by the combined fitting of all scaled capillary pressure and
permeability data to the van Genuchten-Mualem relationship, which will be discussed later. Also
the 9r-values used in the presentation of the scaled capillary pressure data in Fig. 3-6 were
estimated from the fitting of the combined van Genuchten-Mualem model.
Estimation of permeability functions
Measured outflow and capillary pressure data were used to estimate relative permeability
functions of the wetting fluid for the different soil and fluid pair systems. Since changes in
outflow are relatively high at the beginning of each pressure step, only data from the time periods
immediately following an increase in non-wetting fluid pressure were used to estimate
permeability functions. The principle of permeability estimation is the same as the direct K(9)-
method as discussed by Eching et al. (1994) from calculation of Pw,toP at the soil-plate interface.
However, since the reported experiments include wetting fluids other than water, the Darcy flux
equation, Eq. [3-2], was rearranged to yield
P = P — d
w,top w,bottom
[3-8]
Thus, the wetting fluid pressure at the soil-plate interface was estimated from the effective
permeability (kw) and thickness (d) of the saturated porous plate in combination with known
measured values of the wetting fluid pressure at the bottom of the plate (PW!bottom) and the
measured drainage rate (qw), after a pressure increment was applied. Although qw is decreasing
with time within one pressure step, we assumed a constant average flux for the small time interval
30
-------
used at the beginning of each pressure step. The effective permeability of the soil was
subsequently estimated using Eq. [3-2] by solving for kw(Se), after substituting the average
drainage rate and the assumed Pw-gradient in the soil core using the measured wetting fluid
pressures in the center of the core and Pw,toP. The estimation of the average wetting fluid pressure
gradient in the soil sample is based on a linear distribution of wetting fluid pressure in the bottom
half of the soil core. While this assumption appears to hold for the finer-textured Columbia soil
(Eching et al., 1994), it may not be satisfactory for the larger applied non-wetting fluid pressures
in the coarser-textured Lincoln soil. Possible errors caused by the constant gradient assumption
can be largely reduced by placing the wetting-fluid tensiometer nearer to the outflow end of the
soil sample, thereby providing a more accurate wetting-fluid pressure gradient representative for
the measured drainage rate. Moreover, the need for a separate calculation of the wetting-fluid
pressure at the soil-plate interface is eliminated if the thick porous ceramic plate is replaced by a
thin nylon porous membrane with a non-wetting fluid entry pressure large enough for the intended
capillary pressure range. The low resistance nylon membrane3 (Table 3-2) causes only minor
differences in wetting fluid pressure across the membrane, and is now routinely used in outflow
experiments.
Finally, the intrinsic permeability was estimated using the pore-size distribution model of Mualem
(1976) to predict the permeability-saturation function (Eq. 2-1 Ib). The scaled capillary pressure
(Fig. 3-6) and estimated permeability data were fitted to the combined van Genuchten-Mualem
model using the same Excel spreadsheet program used for fitting capillary pressure-saturation
data. As in the estimation of the capillary pressure-saturation curve, average saturation was
obtained from cumulative outflow and initial saturation data. The fitting parameter values of 9r, a,
n, and k (Table 3-2) together describe the capillary pressure and permeability functions for all
three fluid pairs, and characterize the pore structure and pore-size distribution of both soils. The
fitted curves in Figures 3-6 (scaled capillary pressure-saturation) and Figure 3-7 (permeability
functions) were obtained using an 1-value of 0.5 and m-values of m=l-l/n. The plotted relative
permeability data for the three fluid pairs were computed using the fitted intrinsic permeability
values. As is shown in Figure 3-7, the kr(Se) estimates for the three fluid-pair combinations for
both soils are well described by a single relationship, thereby indicating that the relative
permeability is independent of the fluids present being a function of the porous medium properties
only.
Table 3-3. Thickness (d), Saturated Hydraulic Conductivity (Ks), Plate Resistance (Rp), and Air Entry Pressure for
Various Porous Materials.
Ceramic
Plastic
Stainless Steel
Nylon
t(cm)
0.74
0.05
0.1
0.01
KS;Water (cm/hr)
0.0498
0.0166
0.0265
0.025
Rp (hr) Air-entry pressure
14.86
3.01
3.77
0.4
1000
400
250
1700
3 MSI Inc, P.O. Box 1046, Westborough, MA 01581-6046
31
-------
1000 n A. Columbia
Curve
• air-water
• air-oil
A oil-water
0.000 0.200 0.400 0.600 0.800
Se of wetting fluid
1.000
300-1
250
B. Lincoln
0.200
0.400 0.600
Se of wetting fluid
0.800
1.000
Figure 3-6. Scaled and fitted capillary pressure head versus effective saturation of wetting fluid for (A) Columbia and (B)
Lincoln soil
32
-------
l.E+00
^- l.E-01
>,
« l.E-02
-------
Permeability functions were also estimated using a parameter optimization procedure (Eching et
al., 1994) for the air-water systems of soil cores having identical bulk densities. In this approach,
the capillary pressure and permeability functions were estimated indirectly by numerical solution
of the water flow equation (Eq. [3-7] with water as the wetting fluid) for the imposed
experimental boundary and initial conditions. Parameters for the capillary pressure and
permeability functions (van Genuchten-Mualem model for air-water system) were optimized by
minimization of an objective function containing the sum of squared deviations between the
measured and simulated flow variables (capillary pressure and cumulative outflow). The resulting
optimized relative permeability functions are included in Fig. 3-7. For both soils, the estimated
permeability data for the air-water system matched the independently-estimated permeability
functions using the parameter optimization approach quite well.
Parameter optimization
Using the described methodology, the inverse parameter estimation procedure was carried out for
all 3 two-fluid (air-water, air-oil, and oil-water) flow systems of the Lincoln and Columbia soil.
Figures 3-8 and 3-9 compare the measured (circles) and optimized (line) cumulative outflow
Q(cm3) and capillary pressure head (hc = Pc/pn2og cm equivalent head of water), for the air-water
(aw), air-oil (as), and oil-water (ow) systems of the Lincoln (Figure 3-8) and Columbia (Figure 3-
9) soils. Overall, the optimized simulations show an excellent match with the corresponding
measurements, indicating that the optimized parameters captured the main features of the
measured flow procedure.
The RMSSR values (Eq. 2-20), the optimization iteration numbers, with respect to the different
initial parameter (IP) values, are listed in Tables 3-4 and 3-5 for each two-fluid flow system for
the Lincoln and Columbia soil, respectively. The initial parameters were chosen to cover a
relatively broad range of soils, especially for the most sensitive parameters a and n. The choice of
the final selected parameters was based on two considerations: (1) selection of the parameter set
with minimum RMSSR value, or (2) selection of the mean parameter set if there are several
parameter sets with identical RMSSR values. The obtained capillary pressure - saturation and
permeability - saturation curves, based on the optimized parameters are summarized in Table 3-6,
and are shown in Figure 3-10 and Figure 3-11 for the Lincoln and Columbia soil, respectively. As
described in Eq. 2-1 la, the capillary pressure function is dependent on and n only, where is
inversely proportional to the non-wetting-fluid entry value, and thus varies among fluids and as
their corresponding interfacial tension values. The n-parameter is inversely proportional to the
soil's pore-size distribution and determines the slope of the capillary pressure curve. Relative
permeability functions (Eqs. 2-1 Ib and 2-1 Ic ) are only dependent on n, if 1 is fixed to a value of
0.5.
Several important concerns are involved in the above results: the choice of the adjustable or
optimizing parameters, the well-posedness of the relevant inverse problem, and the scaling
relationship of the resulting inverse parameter estimation. In total, the parametric models of
34
-------
180
60
90
120
150
150
° Observed
— Optimized
180
60
150
90
120
150
Time (hour)
Figure 3-8. Comparison of the measured and optimized hc and Q - values of Lincoln soil: (a) air-water system, (b) oil-
water system, (c) air-oil system.
-------
i
E
o
100
100
60
90
120
150
Soooooo oooo
° Observed
— Optimized
60
90
120
150
Time (hour)
Figure 3-9. Comparison of the measured and optimized hc and Q - values of Columbia soil: (a) air-water system, (b) oil-
water system, (c) air-oil system.
36
-------
Table 3-4. Optimized VG-Mualern Parameters of Lincoln Soil for Different Initial Estimates.
System
Air-
Water
Oil-
Water
Air-
Oil
Set Parameter
1 6r(cm3/cm3)
k (10"9 cm2)
a ( cm"1 )
n (-)
2 er
k
a
n
3 er
k
a
n
i er
k
a
n
2 er
k
a
n
3 er
k
a
n
i er
k
a
n
2 er
k
a
n
3 er
k
a
n
IP1 Itera RMSSR2
tions
0.11 8 1.656
5.08
0.04
2.00
0.05 7 1.656
2.00
0.02
3 .00
0.05 12 1.656
8.00
0.05
4.00
0.11 10 2.663
5.08
0.04
2.00
0.05 9 2.640
2.00
0.03
1.80
0.10 7 2.640
4.00
0.05
4.00
0.11 10 3.760
6.11
0.04
2.00
0.05 8 3.830
2.00
0.02
3 .00
0.01 12 3.760
10.00
0.06
4.00
FP3
0.0212
24.8
0.0189
2.8124
0.0211
24.7
0.0189
2.8120
0.0210
24.7
0.0189
2.8098
l.lle-5
9.4
0.0366
2.9472
1.5E-7
9.3
0.0365
2.9660
0.0143
9.0
0.0350
3.3070
2.02e-5
10.9
0.0519
3.0590
1.85e-5
8.4
0.0525
3.0360
3.2e-6
11.3
0.0538
2.9450
Si4
0.003
0.708
0.000
0.040
0.003
0.706
0.000
0.040
0.003
0.704
0.000
0.040
0.008
0.445
0.001
0.079
0.008
0.361
0.001
0.097
0.008
0.488
0.001
0.082
0.011
0.170
0.001
0.080
0.011
0.128
0.001
0.008
0.010
0.211
0.001
0.071
CVi5(%)
14.151
7.9266
0.0000
1.4223
14.218
7.9236
0.0000
1.4225
14.286
7.9048
0.0000
1.4236
.*
13.404
2.7322
2.6805
-
9.9730
2.7397
3.2704
-
16.267
2.8571
2.4796
-
17.375
1.9268
2.6152
-
16.929
1.9048
0.2635
-
17.653
1.8587
2.4109
Fixed parameters: 6S = 0.33 (cm3/cm3), m = 1-1/n, 1 = 0.5
1 IP : Initial Parameter
2 RMSSR : Root of Mean Sum Square Residual
3 FP : Final Parameter obtained from the inverse parameter estimation
4 Si: Standard deviation of parameter bi
5 CVi: Coefficient Variance of parameter bias defined by CVi = (Si/ FPi)-100
37
-------
Table 3-5 Optimized VG-Mualern Parameters of Columbia Soil for Different Initial Estimates.
System Set
Air-
Water
Oil-
Water
Air-
Oil
Parameter IP1
1 6r(cm3/cm3)
k (10~9 cm2)
a ( cm"1 )
n (-)
2 er
k
a
n
3 er
k
a
n
i er
k
a
n
2 er
k
a
n
3 er
k
a
n
i er
k
a
n
2 er
k
a
n
3 er
k
a
n
Iterations RMSSR2 FP3
0.11 12 2.919
5.08
0.04
2.00
0.05 8 2.919
2.00
0.02
3 .00
0.05 10 2.919
4.00
0.05
4.00
0.11 6 3.155
5.08
0.04
2.00
0.10 5 3.155
4.00
0.02
1.20
0.01 8 3.155
8.08
0.05
4.00
0.11 7 3.854
5.08
0.04
2.00
0.05 7 3.854
2.00
0.02
3 .00
0.05 5 3.949
7.00
0.05
4.00
s*
0.0917
5.3
0.0099
2.1474
0.0914
5.3
0.0099
2.1447
0.0912
5.2
0.0100
2.1433
0.0724
8.0
0.0239
2.0372
0.0723
8.0
0.0239
2.0367
0.0723
8.0
0.0239
2.0368
0.0612
6.2
0.0253
2.6939
0.0613
6.2
0.0253
2.6945
0.0485
5.6
0.0248
2.6175
CVf
0.008
0.186
0.000
0.049
0.008
0.184
0.000
0.049
0.008
0.184
0.000
0.049
0.006
0.306
0.000
0.031
0.006
0.307
0.000
0.031
0.006
0.306
0.000
0.031
0.006
0.172
0.000
0.048
0.006
0.172
0.000
0.048
0.004
0.037
0.000
0.038
(%)
8.7241
9.7546
0.0000
2.2818
8.7527
9.6689
0.0000
2.2847
8.7719
9.6832
0.0000
2.2862
8.2873
10.5765
0.0000
1.5217
8.2988
10.6243
0.0000
1.5221
8.2988
10.5996
0.0000
1.5220
9.8039
10.4299
0.0000
1.7818
9.7879
10.4198
0.0000
1.7814
8.2474
2.9152
0.0000
1.4518
Fixed parameters: 6S = 0.45 (cm /cm ), m = 1-1/n, 1 = 0.5
1 IP : Initial Parameter
2 RMSSR : Root of Mean Sum Square Residual
3 FP : Final Parameter obtained from the inverse parameter estimation
4 Si: Standard deviation of parameter bi
5 CVi: Coefficient Variance of parameter bias defined by CVi = (Si/FPi)-100
38
-------
Table 3-6. Optimized VG-Mualem Parameters of Lincoln and Columbia Soil
^\
Lincoln
Columbia
Parameter
9r (cm3/cm3)
k(10'9cm2)
a (cm"1)
n(-)
9r
k
a
n
Air-Water
0.0210
24.8
0.0189
2.8111
0.0913
5.3
0.0100
2.1451
Oil-Water
0.0072
9.4
0.0358
3.1365
0.0723
8.0
0.0239
2.0369
Air-Oil
0.00001
10.9
0.0529
3.002
0.0613
6.2
0.0253
2.6942
Eq. (2-11) requires seven parameters (0S, 9r, k, a, n, m, 1). The choice of how many and which
parameters to optimize was based on a number of considerations. First, it must be recognized
that as the number of optimized parameters increases, the parameter estimation procedure will
generally lead to better matching of optimized with measured flow variables. However, this
occurs at the expense of the uniqueness of the optimized solution and a corresponding increase in
the uncertainty of the parameter estimates. Thus, it is preferred to minimize the total number of
free parameters. If any of the parameters can be measured independently with relatively high
accuracy, it should be fixed rather than optimized. This is certainly the case for the saturated
wetting-fluid content (9S), which was obtained experimentally for each soil core from cumulative
outflow, initial fluid saturation and oven drying of the soil core. Also, the intrinsic permeability
(k) can be used as a fixed parameter in the parameter optimization procedure. However, fixing its
value appears to overly constrain the permeability relationship, allowing it to be a function of n
(or) m only. Moreover, differences in pore geometry between the larger pores (defining soil
structure) and soil matrix pores (defined by soil texture) preclude the use of relationships (2-lib
and c) across the whole saturation range, using a physically-based intrinsic permeability value as
measured from a saturated soil (Demond and Roberts, 1993). Thus, the largest soil pores
(macropores) are not viewed as being predictive of the hydraulic properties of the bulk soil
matrix. Following Mualem (1976), we also reduced the number of optimized parameters by
setting 1 = 0.5 in the permeability function (2-1 Ib and c) for both the wetting and non-wetting
fluids. Finally, we used the relationship between n and m to reduce the number of optimized
parameters further. As determined by Toorman et al. (1992), this final set of four free parameters
is sensitive to the experimental conditions of the multi-step outflow experiment, when cumulative
outflow measurements are supplemented with capillary pressure measurements. Also
39
-------
Air-Water
Soltrol-Water
Air-Soltrol
o>
E
o
3
(A
(A
2!
Q.
o.
re
O
0.4 --
0.2 -•
--0.8
-- 0.6
-- 0.4
--0.2
0.2 0.4 0.6 0.8
Sew (Effective wetting fluid saturation)
Figure 3-10. Optimized constitutive functions of Lincoln soil corresponding to the parameter listed in Table 2.2: (a)
Individual capillary pressure functions, (b) Scaled capillary pressure functions, (c) Relative permeability functions.
40
-------
Air-Water
Soltrol-Water
Air-Soltrol
0 0.2 0.4 0.6 0.8 1
Sew (Effective wetting fluid saturation)
Figure 3-11. Optimized constitutive functions of Columbia soil corresponding to the parameter listed in Table 2.2: (a)
Individual capillary pressure functions, (b) Scaled capillary pressure functions, (c) Relative permeability functions.
41
-------
Finsterle and Pruess (1995) showed that k, a, and n are the most sensitive parameters in two-fluid
flow modeling. In Tables 3-4 and 3-5, the terms S; and CV; are the standard deviation and
coefficient variance of the parameter b;, which were used to measure the estimation accuracy as
well as the sensitivity of parameter b; in the inverse parameter estimation. Of the four adjustable
VG-M parameters (9r, k, a and n), a and n are the most sensitive and accurate parameters.
The optimization results were obtained by imposing constraints on the allowable ranges of the
adjustable parameters. These ranges were set to allow the maximum possible flexibility of the
parameter values, yet limit them so that their values remain to have physical meaning. The
limitation on determining a global feasible region stems from the intrinsic limitation of global
nonlinear optimization theory. The determination of the feasible range for each optimized
parameter is done a priori by a trial-error procedure.
The well-posedness analysis of the inverse parameter estimations was performed by carrying out
the optimizations for different initial parameter estimations for each of the two-fluid flow systems.
According to Russo et al. (1991), the evaluation of the well-posedness of an inverse problem can
be obtained by solving the problem several times with different initial parameter estimates. The
restriction on a priori evaluation of the posedness of an inversion problem stems from the fact that
it is generally impossible to examine whether the converged solution corresponds to global
minimum of the objective function. The final parameters in Tables 3-5 and 3-6, obtained from the
optimizations started at different initial parameter estimates, to test the uniqueness of the solution.
The choice of the initial parameter estimations was based on the consideration to cover the largest
feasible range of the adjustable parameters, especially for the most sensitive parameters (a and n)
of the constitutive function models. The tests show that the feasible initial parameter estimates
were more limited for the single-fluid case, indicating a higher constraint on the two-fluid flow
constitutive relationships than for single fluid flow. The low estimation errors (RMSSR values)
also indicated an acceptable fit of the data, considering that the minimum value of the RMSSR is
1.0, if the objective function is a WLS problem. Consequently, we conclude that the presented
two-fluid flow inverse parameter estimation of a transient multi-step outflow experiment is well-
posed.
For the Leverett assumption to apply to each of the two-fluid flow systems, the n-value for the
three fluid pairs should be identical for each soil, since it is determined by soil pore geometry only.
The permeability functions in Figures 3-10c and 3-1 Ic, which are dependent only on the n-values,
show indeed that they are similar for the three two-fluid flow systems. Furthermore, scaling
requires that the a-values be inversely proportional to the interfacial tension values. Indeed, we
found that the ratios of interfacial tension values (Table 2-2) were quite similar to the a-ratios
from scaling (Table 3-7).
42
-------
Table 3.7 Comparison of a. Ratio and Inter facial-Tension Ratio.
Lincoln Soil Systems
Columbia Soil Systems
a-Ratio a-Ratio
/Gaw = 0.380 0,aw/0,ao = 0.373
/Gaw = 0.534 aaw/a0w = 0.514
a-Ratio a-Ratio
/Gaw = 0.351 aaw /Otao = 0.400
/Gaw = 0.380 aaw/a0w = 0.417
To further examine the scaling relationships between the two-fluid flow systems of the same soil,
we included the scaled capillary pressure functions in Figure 3-10b (Lincoln) and 3-1 Ib
(Columbia), using the capillary pressure curve of the air-water system as a reference. The scaling
factor was simply determined from the corresponded interfacial tensions by using the
dimensionless Leverett's function (2-23):
c,aw(Sew) = hc,aw(Sew)
(3-9a)
(3-9b)
hc,awWew)l3 ~~
>-hc,ow(Sewi
(3-9c)
The subscripts i = 1, 2, and 3 denote the scaled air-water curves hcaw(Sew)|;, obtained from the
optimized air-water, air-oil and oil-water curves, respectively. The scaled curves coalesce well in
the high saturation range. The deviations in the low saturation ranges coincide with the findings
of Demond and Roberts (1991), who indicated that deviations might be caused by limitations of
the traditional Leverett's (1941) scaling function. Demond and Roberts (1991) included other
measurements such as intrinsic contact angle and roughness to correct the scaling factor, thereby
improving the match at the low saturation range. Nevertheless, the close agreement between our
scaling ratios with the interfacial tension ratios attests to the accuracy of the parameter
optimization method.
43
-------
Chapter 4
Conclusions
Multi-step outflow experiments have been successfully used for indirect estimation of capillary
pressure and permeability functions for air-water systems. Here, we extend the application of the
multi-step method to two-fluid phase systems in general, and present results for air-water, air-
Soltrol and Soltrol-water systems in a Columbia fine sandy loam and Lincoln sandy loam for
direct estimation. Subsequently, the experimental data were used in a parameter optimization
algorithm using an inverse technique, thereby providing an indirect method for estimating capillary
pressure and permeability functions. Because of the transient nature of the multi-step outflow
experiments, capillary pressure and permeability functions can be obtained much faster than by
conventional equilibrium methods. This aspect is especially useful when several capillary
pressure and permeability functions are needed to characterize heterogeneous contaminated sites
with large soil spatial variability.
When the experimental capillary pressure-drainage data for the three two-fluid systems are
compared, the interfacial tension of each fluid pair is important because the capillary pressure
value at a given degree of saturation decreases with decreasing interfacial tension. Therefore, non-
wetting fluid pressure increments for each fluid pair were based on the ratio of the interfacial
tension relative to that of air-water, so as to obtain approximately equal amounts of wetting fluid
drainage for each fluid pair combination. The oil-water experiments showed that the oil pressure
remained constant within a pressure increment and was equal to the applied oil pressure, adjusted
for its hydrostatic pressure in the soil core. We therefore conclude from the presented data that a
single-fluid flow model is sufficient for the description of water flow in multi-step outflow
experiments with oil present as a non-wetting fluid.
Using the multi-step outflow method, capillary pressure-saturation and permeability functions are
directly estimated from the experimental data; i.e, from capillary pressure measured in the draining
soil core and from cumulative drainage of the wetting fluid. The accuracy of the directly
measured capillary pressure functions was determined from the success by which interfacial
tension ratio values coalesced individual capillary pressure-saturation curves to a single curve.
Assuming that the scaling factor is only dependent on the interfacial tension of each fluid pair,
scaled capillary pressure-saturation curves were in good agreement with the modified Leverett's
scaling relationship, especially at high saturation values, though some discrepancies occurred at
low wetting-fluid saturations. The combined relative permeability data of both soils coalesced to a
single van Genuchten-Mualem permeability function, thus confirming their independence of the
fluids present in the porous medium. It appears that the employed assumptions with regard to the
wetting-fluid pressure gradient in the draining soil core could be largely removed by inserting the
wetting-fluid tensiometer nearer to the outflow end of the soil core and by replacing the ceramic
porous plate with a thin, nylon porous membrane.
44
-------
This study also demonstrated the feasibility of using the inverse parameter estimation for two-fluid
flow systems. Using the proposed indirect approach the governing flow equations for both the
wetting and non-wetting fluid are solved, so that no assumptions are needed with regard to the
influence of the non-wetting fluid on outflow or wetting fluid pressure gradients. Even though it is
impossible to remove the inherent uncertainty that stems from the high nonlinearity and
complexity of soil systems, the posteriori analysis of the presented inversion problem indicates
that the inverse parameter estimation of the constitutive functions from transient multi-step
outflow experiments of a two-fluid soil system is a well-posed problem. Of the four adjustable
VG-M parameters 9r, k, a and n, the parameters a and n are the most sensitive for the inverse
parameter estimation approach. The selection of proper initial parameter estimate has been shown
to be important for successful optimization. The comparison of the results with those from using
the Leverett scaling method provided a means to test the presented solutions. The advantage of
the proposed method is (1) that the measurements are simple and accurate, (2) simultaneous
estimation of capillary pressure and permeability functions are obtained from a single sample, and
(3) parameter estimation using a flow simulator is consistent with computer modeling flow
simulations.
45
-------
Appendix A.
Optimization Algorithm
Levenberg-Marquardt Method
Generally, a solution algorithm searches the solution for the objective function:
O(b) = eTV 'e = eTWe (A-l)
where, O(b) denotes the objective function, e = vm - vs is the observation error vector with vm
and vs denoting the measured and simulated variables, and W = V"1 is the weighting matrix, with
V the covariance matrix of the error e defined by V = E(vm - vs)(vm - vs)T . The Levenberg-
Marquardt (LM) method has become a standard solution algorithm for nonlinear least-squares
problems. Basically, the LM method is a Newton-type minimization method in which the objective
function is locally approximated to a quadratic form (Press et al., 1992):
0(b) = 0(bi) + (b-bi)-VO(bi) + |(b-bi)-H-(b-bi) (A-2)
or VO(b) = VO(bi) + H-(b-bi) (A-3)
where, H is the Hessian matrix, the second derivative of objective function O(b), with
^2n
H = - . In Newton's method, we set the left-hand side of (A-3) equal to zero in order to
11 «.• «.• " \ / 1
dbtdb.
determine the next iteration point, or:
(A-4)
with b[+l=b[+Ab (A-5)
Substitution of (A-l) into (A-4) and defining Ve = J (Jacobian or sensitivity matrix), with Jy =
de;/dbj, one obtains:
H Ab = - JTWe (A-6)
Since the Hessian matrix is difficult to calculate, one uses a Hessian approximation
H = JrWJ (Kool and Parker, 1988) in (A-6) to yield the Gauss-Newton algorithm:
JTW J • Ab = - JTWe (A-7)
46
-------
By inspection of (A4), for the optimization to proceed in a descending direction, H or its
approximation JTWJ must be positive-definite. This is ensured in the Levenberg-Marquardt
algorithm by adding a positive quantity to JTWJ:
(JTWJ + ?J)TD)Zib = - JTWe (A-8)
where A, is a positive scalar or Levenberg parameter, and D is a diagonal scaling matrix with
elements equal to the norms of the corresponding columns of J (More, 1977). Marquardt (1963)
developed an effective strategy to update A,. When far from the minimum, A, is given a large
value, yielding a step in the steepest descent direction, whereas, when approaching the optimum,
A, is given a small value, so that (A-8) converges to a Gauss-Newton step.
The parameter optimization problem formulated by the objective function (A-l) and solved by the
iteration solution (A-8) is referred to as unconstrained optimization. Usually, there are physical
and mathematics constraints for the estimated parameters. With the added parameter bounds, the
corresponding optimization is therefore referred to as bounded optimization. When the
parameters are outside the constrains, they will be forced back to the bounded region. Kool and
Parker (1988) reported on a bounded iteration strategy.
Convergence and Reliability
The iteration solution of optimization (A-8) is terminated by convergence criteria. The commonly
used stopping criteria include two types of tests. The first one is based on the magnitude of the
RMSSR as defined by:
4RMSSR
or
RMSSR + sa '
whereas the second criterion was the relative change in parameter values:
"; ' or — ^2 (A-ll)
where sa and s\, are small values ensuring that the denominators are not equal to zero, TI, TI' and
i2 are convergence accuracy tolerances. Usually, the second criterion (A-ll) is tested after the
first criterion (A-10) is satisfied. Only using the second criterion (A-ll) or meeting a small step
Ab does not guarantee that the solution is at a minimum, since a large value for A will also
produce a very small step Ab. The accuracy tolerance T/ is equal to 0 in order for the RMSSR
47
-------
to change toward to the decrease direction. The accuracy tolerance T2 is often problem dependent
and is a compromise between estimation accuracy and computational expense. For example,
when the level of uncertainty in input data is increased, objective function has flat minimum, the
parameters will tend to wander around near the minimum. In that case, the smaller convergence
criterion has little effect on estimation accuracy, but leads to additional iterations, thereby
significantly increasing computational expense. An accuracy tolerance T2 = 0.01 is chosen in our
study.
The uncertainty measurement of the estimated parameters is expressed by their confidence region,
which is derived from linear regression analysis under the assumption of normality and linearity.
The normality assumption is that the distribution of a sum of random variables always tends
towards normal if the sample size is sufficiently large. The implication is that the measurement
errors are dependent on a linear combination of a large number of small random factors. The
linearity assumption is that nonlinear functions of parameters b can be approximated by a
linearization within the confidence region. For a maximum likelihood estimator, the parameter
covariance matrix is asymptotically given by :
C = s02H(b)-1 (A-12)
T
or under Hessian approximation of H = J WJ , get:
C = s02(JTWJ)"1 (A-13)
, eTWe
with s02 = ^-^ (A-14)
n- m
where, the circumflex indicates a posterior value, n is the total number of observations, m is the
number of parameters to be estimated, and s02 is the estimated residual variance at the optimum.
Consequently, the standard deviation s; of the parameter b; and the confidence region can be
described by (Kool and Parker, 1988):
where, a is the probability that hypothesis is rejected even though it is true, Q; is the parameter
variance, t = tv,i-o5a is the value of Student's t-distribution for confident level 1-a and v-degrees of
freedom. For example, with a 95% confidence level, the boundaries of the confidence region
are:
bi,min=bi,opt-tv,o.975Si (A- 17)
bi,max=bi,opt+tv,o.975Si (A- 18)
where b;,opt is the optimized value of the parameter b;.
48
-------
Appendix B
Description of Two-fluid Flow Model (TF-OPT)
Introduction
What is TF-OPT and when is it useful?
TF-OPT is a software specifically developed for estimation of the constitutive functions (capillary
pressure and permeability) of two-fluid flow in a soil core from multi-step outflow experiments.
TF-OPT consists of two parts: a one-dimensional two-fluid flow model with a finite element
scheme and an optimization algorithm using Levenberg-Marquardt (LM) method.
Which steps are needed for using TF-OPT?
• Install TF-OPT: TF-OPT is written in FORTRAN and has been run on PC-compatibles and
UNIX environments. The software can be downloaded via anonymous FTP from
theis.ucdavis.edu (128.120.37.3) by using the user's E-mail address as password. It is
located in the directory /pub/out/tf-opt. The files in tf-opt/ are:
README Explanation file
example.in Example input file
example.out Example output file
TF-OPT.for TF-OPT FORTRAN code file
tfcombk.dat Include file of TF-OPT
• Prepare input file:
The efforts have been made to make the input file user-friendly and self-explainable. However, the
setup of the system conditions (1C and BC and choice of finite element nodes), appropriate
selection of the parametric model of capillary pressure and permeability functions, choice of initial
parameter estimates and parameter limits are determined by user.
• Search for successful optimization solution:
"As well been known, it is generally more difficult to find solutions to nonlinear
optimization problems than to linear ones. For difficult nonlinear problems, users usually
have to pay much more attention to apparently inconsequential details than they would like. "
— GAMS Release 2.25
It has been demonstrated that the inverse parameter estimation of the multi-step outflow
experiment of soil column is a "well-posed1 problem. After preparing the input file, changes might
be needed. For example, you may need to change the parameter-limits, or initial parameter
49
-------
estimates, or adjust the data-type weighting factors. A proper choice of the parametric model for
the capillary pressure and permeability function is also needed to obtain a successful solution.
Example
A-2-1 Input file
INPUT FILE
PART I: MODELING DATA
Total-Nodes Plate-nodes Obs-node Soil-L Plate-L Core-Dia. (L in cm)
(Plate-node and Plate-L is set as 0 if a thickless membrane is used)
60 6 33 7.64 0.74 6.
Ceramic-Plate Hydraulic-Property (any value for a thickless membrane)
Theta-s Ks
0 . 506 0.0477
Fluid Property:
RhoW RhoNW MuW MuNW (Rho:density in g/cm3; Mu:viscosity in dyn.s/cm2)
1. 0. l.E-2 1.81E-4
Tmax(hr) DTmax IDT(=1,CONST.DT) INIDT Err
170. 5. 1 0.1 0.1
IEQ (1-VGM, 2-VGB, 3-BCM, 4-BCB, 5-BRB, 6-GDM, 7-LNM)
1
MODE (l:TF_Simu., 2:TF_opti.; -1:ESF_Simu., -2:ESF_Opti.):
(TF-two_fluid flow; ESF-Equivalent single_fluid flow)
2 ~ ~
hnwO (Initial NW total head at top), HBO (Initial outflow height)
20.6 0.
Step number of Multi-step pressure: ITIM
8
B.C. : time(hr),applied pressure (pb(i),i=l,itim,in cm water)
0.0 25.5
5.67
21.25
45.3
69. 05
92.97
122.6
150.9
PART II: OPTIMIZATION DATA
MAXTRY MIT
15 15
NTOB(Time Points), IPAR (Parameter Number),ITYP (Data-Type Number)
229 8 3
Data-type weighting factor
IWHC IWQ IWTHETA
1 1 30
PARAMETER NAMES
THETAS(cm3/cm3)
THETAR(cm3/cm3)
INTRIK(cm2)
ALPHA(-)
N(-)
M(-)
L(-)
NONE
PVALI(I) ,PMINI(I) ,PMAXI(I) ,IOPT(I) ,1 = 1,I PAR
(Initial estimate, Limits(MIN, MAX), Fixed:O/Free:1)
0.32 0.25 0.362 0
0.11 0.000 0.3 1
0.000000014399 0.00000000141723 0.000000141723 1
0. 04 0.001 0.5 1
2.0 1.1 10.1
1. 1. 1.0
50
-------
0.5 1.
Observation data points:
Time (hr ) ,
1 .
4
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
6
6
6
6
6
6
6
6
6
6
T
T
T
•-j
a
8
9
9
10
11
12
12
14
16
19
20
21
21
21
21
21
21
21
21
21
21
21
0167
0334
0500
0667
0834
1167
2667
5334
1667
4334
6834
7000
7167
7334
7500
7667
7834
8000
8167
8500
8834
9167
9500
9834
0167
0667
1500
^334
3000
3834
4834
6334
7500
9000
1167
3667
6500
9500
3334
7500
2667
6167
2500
3167
2667
7334
7167
^334
4167
6167
2667
2834
3000
3167
3334
3500
3667
3834
4000
4167
4334
he
2 4
25
2 6
2 6
27
27
28
28
29
29
33
34
35
35
36
36
37
37
38
38
39
39
39
40
40
40
41
41
41
4 2
4 2
4 2
4 2
4 2
4 2
43
43
43
43
43
43
43
43
43
43
43
43
43
43
43
43
46
47
48
49
49
50
51
51
52
52
cm] , Q ( cm3
2000
4400
1500
6000
0000
4900
4200
9600
6200
6200
4500
2500
1800
7600
3800
8700
3600
7600
1600
8300
3200
7600
7600
1600
4800
8300
3700
8100
9000
2100
4300
4800
7400
8800
9700
0600
2300
3200
3200
3200
3200
3200
3200
3200
3200
3200
3200
3200
3200
3200
3200
2000
4000
4200
2200
9800
1500
0900
6200
0700
5100
1
1
1
1
1
1
2
2
2
2
2
2
^
^
^
^
^
^
4
4
4
4
5
5
5
5
6
6
6
6
6
~->
~->
~->
~->
8
8
8
8
8
9
9
9
9
10
10
10
10
10
11
11
12
12
12
13
13
13
14
14
14
1. 0
, HB(cm)
8300
1500
3400
5200
6600
8000
9800
0700
2100
3000
5300
6800
9200
1500
3000
4500
6100
8400
9900
2200
5300
6800
9100
1500
3700
6000
8300
0600
3000
5200
7500
9900
2200
4500
6000
9100
1400
3700
6000
8300
9900
2900
5200
7500
9000
2100
3600
5200
9000
9800
5000
8300
3300
6600
9900
3100
6400
9700
2200
4600
7100
.1500 !' °
.2100
.2400
.2700
.3000
.3200
.3600
.3700
.4000
.4100
.4600
.4800
.5200
.5700
.5900
. 6200
.6500
.6900
.7200
.7600
.8100
.8400
.8800
. 9300
. 9700
1.0100
1.0500
1 . 0900
1.1300
1.1700
1.2200
1 . 2600
1.3000
1.3400
1.3700
1.4200
1.4600
1.5100
1.5500
1.5900
1.6200
1.6700
1.7100
1.7500
1.7800
1.8400
1.8700
1 . 8900
1. 9600
1. 9800
2.0700
2.1300
2.2200
2.2800
2.3400
2. 4000
2.4500
2.5100
2.5600
2. 6000
2.6500
51
-------
94
94
95
97
101
102
102
104
117
122
122
122
123
123
124
124
125
126
139
146
150
151
151
152
153
155
157
162
167
168
16.
.3000
.6000
.0000
7334
.2000
.5167
.3667
.7500
.7334
.1167
.5167
. 6667
.8000
.0834
.4834
. 2667
.7000
.4500
.4167
5500
.2334
.0834
. 9000
.1000
. 2500
.6000
5500
5500
. 7834
.8834
. 5500
. 9167
0000
120.
122 .
125.
ill:
134 .
135.
143.
144.
144.
145.
148.
148.
148.
149.
151.
154.
159.
161.
164.
168.
181 .
188.
192.
193.
198.
199.
205.
210.
219.
244 .
257 .
261.
^
3400
9600
9900
§^88
2600
7700
3800
1300
4900
7400
0000
4000
6100
0600
2800
3000
1500
4200
8900
5800
6100
4200
2000
8900
9500
3800
6400
0800
0700
7900
7900
6900
2000
080
47
47
48
U
48
48
49
49
49
50
50
50
50
50
50
50
51
51
51
51
51
51
52
52
53
53
53
53
53
54
7400
9900
2400
1188
8200
9900
1600
4900
6600
0000
0000
0000
2000
3900
5900
7800
0400
1600
2900
5500
6800
9300
2000
2000
4300
5^00
8300
0800
3^00
5500
8000
Q500
2000
8
8
8
8
8
8
8
8
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
5900
6400
6800
7900
8200
8500
9100
9400
0000
0000
0000
0400
0700
1100
1400
1900
2100
2300
2800
3000
3500
4000
4000
4400
4500
5100
5500
6000
6400
6800
7100
7600
52
-------
Line Variable
Description
5 NNP Total finite element nodes
IFNODE Interfacial node between the soil core and ceramic plate
( = 0 if a nylon membrane is used instead of a ceramic plate)
IOBSNODE Observation node
SLENTH Soil core length
PLTLENTH Ceramic plate thickness (= 0 for nylon membrane case)
DIAM Soil core diameter
8 PLTPROP PLTPROP( 1) = K, and PLTPROP(2) = 6S of ceramic plate
11 RHOW Wetting fluid density (g/cm3)
RHONW Non-wetting fluid density (g/cm3)
CMUW Wetting fluid viscosity (dyne sec/cm2)
CMUNW Non-wetting fluid viscosity (dyne sec/cm )
13 TMAX Maximum simulation time (user defined unit)
DTMAX Maximum time step for used in the numerical scheme
IDT Index of time step, equal to 1 for the constant time step case
INIDT Initial time step
ERR Numerical solution accuracy
15 IEQ Index of the constitutive functions (hc-Se and kr-Se) model:
1: van Genuchten - Mualem model (VGM)
2: van Genuchten - Burtine model (VGB)
3: Brooks-Corey -Mualem model (BCM)
4: Brooks-Corey-Burdine model (BCB)
5: Brutseart-Burdine model (BRB)
6: Gardner-Mualem model (GDM)
7: Lognormal-Mualem model (LNM)
18
20
MODE Index of calculation type:
1: Two-fluid flow forward simulation
2: Two-fluid flow inverse parameter estimation
hnwO Initial non-wetting fluid total head at top
HBO Initial outflow height
ITIM Total number of multi-step pressure step
22
24 - 31 TIM(I), PB(I), 1=1,..., ITIM, time moment and value of each pressure step
34 MAXTRY
MIT
36 NTOB
IPAR
ITYP
39 IWHC
IWQ
IWTHETA Weighting factor for the 6w(hc) data
Maximum running number in each optimization iteration
Maximum iteration number
Total observation time points
Total parameter number (maximum is 8)
Total number of data type
Weighting factor for capillary pressure hc data
Weighting factor for cumulative outflow Qw data
41-43 THETAS, THETAR, INTRIK: Fixed names for parameter 1, 2 and 3
44-48 ALPHA, N, M, L, NONE: User defined name for parameter 4, 5, 6, 7, and 8
51-58 PVALI(I), PMIN(I), PMAX(I), IOPT(I), I = 1,..., IPAR, initial estimates,
minimum, maximum, and adjustable index (1: free, 0: fixed) of parameters
61-End Time(I), hc(I), Q(I), HB(I), I = 1,..., NTOB, Observation data
End-line: Initial hc and 6W data when soil core is in equilibrium state
53
-------
Output file
INITIAL OBS-CAL FITTING:
TIME
0.
0.
0.
0.
0.
0.
0.
0.
4.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
7.
7.
7.
7.
8.
8.
9.
9.
10.
11.
12.
12.
14.
16.
19.
20.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
OBS-HC
,017
,033
,050
,067
,083
,117
,267
,533
,167
,433
,683
,700
,717
,733
,750
,767
,783
,800
,817
,850
,883
,917
,950
,983
,017
,067
,150
,233
,300
,383
,483
,633
,750
,900
,117
,367
,650
,950
,333
,750
,267
,617
,250
,317
,267
,733
,717
,233
,417
,617
,267
,283
,300
,317
,333
,350
,367
,383
,400
,417
,433
,450
24.
25.
26.
26.
27.
27.
28.
28.
29.
29.
33.
34.
35.
35.
36.
36.
37.
37.
38.
38.
39.
39.
39.
40.
40.
40.
41.
41.
41.
42.
42.
42.
42.
42.
42.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
46.
47.
48.
49.
49.
50.
51.
51.
52.
52.
52.
CAL-HC DFHC OBS-Q
,200
,440
,150
,600
,000
,490
,420
,960
,620
,620
,450
,250
,180
,760
,380
,870
,360
,760
,160
,830
,320
,760
,760
,160
,480
,830
,370
,810
,900
,210
,430
,480
,740
,880
,970
,060
,230
,320
,320
,320
,320
,320
,320
,320
,320
,320
,320
,320
,320
,320
,320
,200
,400
,420
,220
,980
,150
,090
,620
,070
,510
,870
25.
25.
25.
25.
25.
25.
26.
27.
29.
29.
29.
30.
30.
30.
31.
31.
31.
32.
32.
33.
33.
34.
34.
35.
35.
36.
37.
38.
38.
39.
39.
40.
41.
41.
42.
42.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
44.
44.
44.
45.
45.
45.
46.
46.
46.
,039
,192
,338
,479
,615
,870
,793
,903
,672
,657
,747
,039
,410
,830
,236
,627
,999
,351
,689
,318
,897
,433
,932
,401
,839
,441
,298
,051
,599
,200
,813
,588
,081
,607
,190
,675
,050
,308
,506
,620
,683
,694
,680
,643
,609
,582
,511
,484
,436
,400
,416
,566
,795
,098
,431
,770
,107
,435
,749
,053
,344
,622
0.
0.
0.
1.
1.
1.
1.
1.
0.
0.
3.
4.
4.
4.
5.
5.
5.
5.
5.
5.
5.
5.
4.
4.
4.
4.
4.
3.
3.
3.
2.
1.
1.
1.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
2.
3.
4.
4.
5.
5.
5.
5.
6.
6.
6.
CAL-Q DFQ
,839
,248
,812
,121
,385
,620
,627
,057
,052
,037
,703
,211
,770
,930
,144
,243
,361
,409
,471
,512
,423
,327
,828
,759
,641
,389
,072
,759
,301
,010
,617
,892
,659
,273
,780
,385
,180
,012
,186
,300
,363
,374
,360
,323
,289
,262
,191
,164
,116
,080
,096
,634
,605
,322
,789
,210
,043
,655
,871
,017
,166
,248
0.
1.
1.
1.
1.
1.
1.
2.
2.
2.
2.
2.
2.
3.
3.
3.
3.
3.
3.
4.
4.
4.
4.
5.
5.
5.
5.
6.
6.
6.
6.
6.
7.
7.
7.
7.
8.
8.
8.
8.
8.
9.
9.
9.
9.
10.
10.
10.
10.
10.
11.
11.
12.
12.
12.
13.
13.
13.
14.
14.
14.
14.
,830
,150
,340
,520
,660
,800
,980
,070
,210
,300
,530
,680
,920
,150
,300
,450
,610
,840
,990
,220
,530
,680
,910
,150
,370
,600
,830
,060
,300
,520
,750
,990
,220
,450
,600
,910
,140
,370
,600
,830
,990
,290
,520
,750
,900
,210
,360
,520
,900
,980
,500
,830
,330
,660
,990
,310
,640
,970
,220
,460
,710
,960
0.
0.
0.
0.
0.
0.
1.
1.
2.
2.
3.
3.
3.
3.
4.
4.
4.
4.
4.
5.
5.
5.
5.
6.
6.
6.
7.
7.
7.
7.
8.
8.
8.
8.
8.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
10.
10.
10.
10.
10.
10.
10.
11.
11.
11.
,075
,168
,257
,343
,427
,583
,140
,791
,789
,780
,066
,372
,645
,892
,118
,328
,524
,705
,877
,186
,467
,723
,957
,173
,372
,640
,016
,335
,561
,806
,052
,355
,544
,743
,960
,138
,273
,365
,435
,474
,496
,498
,492
,477
,464
,453
,426
,416
,397
,382
,709
,948
,139
,305
,451
,581
,700
,809
,909
,004
,093
,177
0,
0,
1,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
3,
3,
3,
3,
3,
.755
.982
.083
.177
.233
.217
.840
.279
.579
.480
.536
.692
.725
.742
.818
.878
.914
.865
.887
.966
.937
.043
.047
.023
.002
.040
.186
.275
.261
.286
.302
.365
.324
.293
.360
.228
.133
.995
.835
.644
.506
.208
.028
.273
.436
.757
.934
.104
.503
.598
.791
.882
.191
.355
.539
.729
.940
.161
.311
.456
.617
.783
54
-------
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
23.
23.
23.
23.
24.
24.
25.
26.
27.
27.
27.
27.
28.
29.
30.
32.
33.
34.
36.
37.
41.
44.
45.
45.
45.
45.
45.
45.
45.
45.
45.
,483
,500
,533
,550
,567
, 600
, 633
, 650
, 683
,700
,750
,783
,817
,850
,883
,933
,983
,033
,067
,083
,100
,117
,133
,150
,167
,183
,217
,283
,350
,417
,533
, 650
,767
,900
,050
,233
,467
,750
,117
, 650
, 600
,167
,150
,433
, 600
,900
,483
,567
,717
,017
,033
,117
,067
, 650
,150
, 633
,317
,333
,367
,383
,400
,417
,450
,483
,517
53.
54.
54.
54.
55.
55.
55.
56.
56.
56.
57.
57.
57.
57.
57.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
59.
59.
59.
59.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
64.
65.
66.
67.
67.
68.
69.
70.
70.
, 670
,030
,510
,780
,000
,490
,940
,120
,430
, 610
,050
,230
,490
, 670
,940
,160
, 610
,740
,740
,740
,740
,740
,740
,740
,740
,740
,740
,740
,740
,830
,360
,540
,720
,900
,160
,250
,340
,520
,520
,520
,520
,520
,520
,520
,520
,520
,520
, 610
, 610
, 610
, 610
, 610
, 610
, 610
, 610
, 610
,370
,220
, 640
,310
,840
,370
,260
,110
,780
47.
47.
47.
48.
48.
48.
49.
49.
49.
49.
50.
50.
50.
51.
51.
51.
52.
52.
52.
53.
53.
53.
53.
53.
53.
53.
53.
54.
54.
55.
55.
56.
56.
57.
57.
57.
58.
59.
59.
60.
60.
60.
61.
61.
61.
61.
61.
61.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
61.
61.
61.
61.
62.
,148
,395
,866
,089
,307
,723
,118
,308
, 675
,853
,359
, 683
,993
,293
,582
,994
,385
,758
,998
,116
,232
,346
,459
,568
, 676
,783
,987
,369
,728
,067
, 610
,102
,550
,012
,475
,967
,498
,023
,552
,100
, 667
,853
,027
,050
,056
,057
,047
,017
,979
,934
,895
,851
,783
,736
, 685
, 648
, 641
, 667
,830
,944
,083
,243
, 610
,992
,370
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
5.
5.
5.
5.
5.
5.
5.
5.
4.
4.
4.
4.
3.
3.
3.
3.
2.
2.
2.
1.
1.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
3.
4.
5.
6.
6.
7.
7.
8.
8.
,522
, 635
, 644
, 691
, 693
,767
,822
,812
,755
,757
, 691
,547
,497
,377
,358
,166
,225
,982
,742
, 624
,508
,394
,281
,172
,064
,957
,753
,371
,012
,763
,750
,438
,170
,888
, 685
,283
,842
,497
,968
,420
,147
,333
,507
,530
,536
,537
,527
,407
,369
,324
,285
,241
,173
,126
,075
,038
,729
,553
,810
,366
,757
,127
, 650
,118
,410
15.
15.
15.
16.
16.
16.
16.
17.
17.
17.
17.
17.
18.
18.
18.
18.
19.
19.
19.
19.
20.
20.
20.
20.
21.
21.
21.
21.
21.
22.
22.
22.
22.
23.
23.
23.
23.
24.
24.
24.
24.
25.
25.
25.
25.
26.
26.
26.
26.
27.
27.
27.
27.
28.
28.
28.
28.
28.
29.
29.
29.
29.
29.
30.
30.
,370
,530
,860
,110
,270
, 600
,850
,000
,330
,420
,830
,990
,320
,480
,730
,970
,300
,470
, 630
,880
,130
,370
,710
,950
,110
,190
,450
, 690
,930
,100
,510
, 670
,920
,160
,410
, 660
,910
,150
,400
, 640
,970
,130
,380
,700
,870
,110
,360
, 610
,850
,100
,350
,590
,930
,090
,340
,500
, 680
,940
,180
,370
,490
, 670
,800
,040
,280
11,
11,
11,
11,
11,
11,
11,
11,
12,
12,
12,
12,
12,
12,
12,
12,
12,
12,
12,
12,
13,
13,
13,
13,
13,
13,
13,
13,
13,
13,
13,
13,
13,
13,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
14,
15,
15,
15,
15,
15,
15,
15,
15,
.333
.406
.544
. 609
. 673
.792
.905
.960
.063
.113
.254
.344
.429
.510
.589
. 699
.802
.900
.963
.993
.023
.051
.079
.106
.133
.159
.210
.305
.394
.478
. 611
.730
.837
.946
.054
.168
.289
.407
.525
. 646
.770
.809
.846
.849
.849
.849
.846
.839
.829
.818
.809
.798
.782
.771
.758
.750
.944
.060
.217
.286
.346
.402
.498
.584
. 661
4.
4.
4.
4.
4.
4.
4.
5.
5.
5.
5.
5.
5.
5.
6.
6.
6.
6.
6.
6.
7.
7.
7.
7.
7.
8.
8.
8.
8.
8.
8.
8.
9.
9.
9.
9.
9.
9.
9.
9.
10.
10.
10.
10.
11.
11.
11.
11.
12.
12.
12.
12.
13.
13.
13.
13.
13.
13.
13.
14.
14.
14.
14.
14.
14.
,037
,124
,316
,501
,597
,808
,945
,040
,267
,307
,576
, 646
,891
,970
,141
,271
,498
,570
, 667
,887
,107
,319
, 631
,844
,977
,031
,240
,385
,536
, 622
,899
,940
,083
,214
,356
,492
, 621
,743
,875
,994
,200
,321
,534
,851
,021
,261
,514
,771
,021
,282
,541
,792
,148
,319
,582
,750
,736
,880
,963
,084
,144
,268
,302
,456
, 619
55
-------
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
46.
46.
46.
46.
46.
47.
47.
50.
53.
53.
55.
57.
61.
64.
67.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
70.
70.
71.
74.
74.
79.
80.
87.
92.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
,533
,567
, 617
, 667
,717
,767
,817
,867
,900
,950
,967
,000
,183
,417
, 667
,900
,233
,983
,867
,150
,367
,117
,533
,233
,100
,883
,033
,067
,083
,117
,150
,183
,217
,233
,267
,317
,400
,483
,533
, 617
,700
,833
,983
,217
,933
,983
,050
,850
,367
,500
, 633
,967
,000
,017
,033
,067
,100
,133
,167
,233
,317
,383
,467
,550
, 633
71.
71.
72.
73.
73.
74.
74.
75.
75.
75.
75.
76.
77.
78.
78.
78.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
80.
80.
81.
82.
83.
84.
84.
85.
86.
88.
89.
90.
91.
91.
92.
93.
95.
97.
98.
98.
99.
99.
99.
99.
99.
99.
100.
100.
101.
102.
103.
104.
105.
107.
108.
110.
111.
112.
,130
,760
,560
,400
,980
,470
,960
,400
, 630
,940
,940
,340
,230
,120
,700
,870
,230
,410
, 630
, 630
, 670
, 670
,720
,720
,720
,720
,810
,060
,590
, 610
, 680
, 610
,460
,900
, 620
, 640
,060
,310
,020
,000
,800
,870
,850
,180
,140
,290
,870
,050
,270
,270
,270
,270
,430
,100
,500
,390
,330
,210
,020
, 620
,440
,730
,200
,530
,730
62.
62.
63.
63.
64.
64.
65.
65.
65.
66.
66.
66.
67.
68.
69.
70.
71.
74.
77.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
80.
80.
80.
80.
81.
81.
82.
82.
82.
83.
84.
85.
87.
90.
93.
94.
97.
97.
98.
98.
98.
98.
98.
98.
98.
98.
99.
99.
99.
100.
100.
101.
101.
,555
,916
,431
,910
,358
,777
,173
,547
,787
,131
,244
,462
,543
,742
,869
,800
,965
,032
,988
,096
,157
,483
, 680
,753
,750
,723
,708
,714
,716
,746
,812
,922
,068
,149
,335
, 648
,197
,732
,041
,531
,991
, 668
,356
,304
, 639
,201
,575
,484
,366
,722
, 679
,837
,834
,836
,839
,858
,897
,962
,053
,337
,797
,204
,747
,301
,853
8.
8.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
8.
8.
7.
5.
1.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
1.
2.
3.
4.
4.
5.
5.
6.
7.
7.
8.
8.
9.
9.
9.
9.
8.
5.
4.
1.
1.
0.
0.
0.
1.
1.
2.
3.
4.
4.
6.
7.
8.
9.
10.
10.
,575
,844
,129
,490
, 622
, 693
,787
,853
,843
,809
, 696
,878
, 687
,378
,831
,070
,265
,378
, 642
,534
,513
,187
,040
,033
,030
,003
,102
,346
,874
,864
,868
, 688
,392
,751
,285
,992
,863
,578
,979
,469
,809
,202
,494
,876
,501
,089
,295
,566
,904
,548
,591
,433
,596
,264
, 661
,532
,433
,248
,967
,283
, 643
,526
,453
,229
,877
30.
30.
30.
30.
31.
31.
31.
31.
31.
31.
32.
33.
33.
33.
33.
34.
34.
34.
34.
34.
35.
35.
35.
35.
35.
35.
36.
36.
36.
37.
37.
37.
38.
38.
38.
38.
39.
39.
39.
39.
40.
40.
40.
40.
41.
41.
41.
41.
41.
42.
42.
42.
42.
43.
43.
44.
44.
44.
44.
45.
45.
45.
45.
46.
46.
,410
,530
,780
,960
,140
,270
,510
, 640
,700
,820
,990
,300
,480
, 660
,910
,090
,210
,390
,580
,700
,010
,130
,320
,560
, 690
,940
,000
,270
, 620
,070
,430
,780
,050
,230
,500
,760
,120
,390
, 660
,920
,090
,370
, 640
,910
,250
,430
,790
,790
,960
,220
,490
,500
,990
,750
,990
,330
,570
,830
,990
,320
,490
,820
,990
,240
,490
15.
15.
15.
15.
16.
16.
16.
16.
16.
16.
16.
16.
16.
16.
16.
17.
17.
17.
18.
18.
18.
18.
18.
18.
18.
18.
18.
18.
18.
18.
18.
18.
18.
18.
18.
18.
19.
19.
19.
19.
19.
19.
19.
19.
19.
20.
20.
20.
20.
20.
20.
20.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
, 698
,765
,858
,944
,023
,097
,166
,232
,274
,334
,354
,390
,576
,784
,975
,129
,319
, 644
,231
,387
,395
,440
,467
,476
,476
,471
,469
,577
, 634
,709
,769
,819
,864
,884
,922
,972
,045
,110
,146
,201
,252
,327
,402
,508
,771
,053
,409
,501
,786
,820
,910
,925
,078
,122
,158
,215
,264
,307
,345
,412
,484
,536
,595
, 648
, 698
14.
14.
14.
15.
15.
15.
15.
15.
15.
15.
16.
16.
16.
16.
16.
16.
16.
16.
16.
16.
16.
16.
16.
17.
17.
17.
17.
17.
17.
18.
18.
18.
19.
19.
19.
19.
20.
20.
20.
20.
20.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
22.
22.
23.
23.
23.
23.
23.
24.
24.
24.
24.
24.
,712
,765
,922
,016
,117
,173
,344
,408
,426
,486
, 636
,910
,904
,876
,935
,961
,891
,746
,349
,313
, 615
, 690
,853
,084
,214
,469
,531
, 693
,986
,361
, 661
,961
,186
,346
,578
,788
,075
,280
,514
,719
,838
,043
,238
,402
,479
,377
,381
,289
,174
,400
,580
,575
,912
, 628
,832
,115
,306
,523
, 645
,908
,006
,284
,395
,592
,792
56
-------
93.
93.
93.
94.
94.
94.
95.
95.
96.
97.
101.
102.
102.
104.
117.
122.
122.
122.
123.
123.
124.
124.
125.
126.
132.
139.
146.
150.
151.
151.
152.
153.
155.
157.
162.
167.
168.
16.
733
850
967
100
300
600
000
617
733
200
517
367
750
733
117
517
667
800
083
483
267
700
450
417
550
233
083
900
100
250
600
550
550
783
883
550
917
000
ITERATION
0
1
2
3
4
5
6
7
0,
0,
0,
0,
0,
0,
0,
0,
114.070
115. 620
117.000
118.340
120.340
122.960
125.990
129.590
134.260
135.770
143.380
144.130
144.490
145.740
148.000
148.400
148. 610
149.060
151.280
154.300
159.150
161.420
164.890
168.580
181. 610
188.420
192.200
193.890
198.950
199.380
205. 640
210.080
219.070
227.790
244.790
257. 690
261.200
0.308
102
103
103
104
105
107
109
111
115
116
127
128
129
132
142
144
144
144
145
146
148
149
151
154
165
173
179
183
183
183
185
187
193
198
209
218
220
0
505
239
942
708
790
280
072
533
383
825
143
691
348
377
770
812
855
908
174
018
348
659
769
178
220
547
793
191
140
242
435
972
455
991
850
553
952
287
SSQ THETAR(cm3/cm3)
.1159E+02
. 6138E+01
.1968E+01
.1675E+01
.1659E+01
.1657E+01
.1656E+01
.1656E+01
MEET MAXIMUM NET, NO
0
0
0
0
0
0
0
0
.1100 0
.0068 0
.0146 0
.0115 0
.0172 0
.0197 0
.0209 0
.0212 0
11.565 46
12.381 47
13.058 47
13.632 47
14.550 47
15.680 47
16.918 48
18.057 48
18.877 48
18.945 48
16.237 49
15.439 49
15.142 49
13.363 50
5.230 50
3.588 50
3.755 50
4.152 50
6.106 50
8.282 50
10.802 51
11.761 51
13.121 51
14.402 51
16.390 51
14.873 51
12.407 52
10.699 52
15.810 52
16.138 52
20.205 52
22.108 53
25.615 53
28.799 53
34.940 53
39.137 53
40.248 54
0.021
.740
.080
.240
.490
.740
.990
.240
.490
.820
.990
.160
.490
. 660
.000
.000
.000
.200
.390
.590
.780
.040
.160
.290
.550
. 680
.930
.200
.200
.430
.520
.830
.080
.320
.550
.800
.950
.200
21.
21.
21.
21.
22.
22.
22.
22.
22.
22.
23.
23.
23.
23.
24.
24.
24.
24.
24.
24.
24.
24.
24.
25.
25.
25.
25.
25.
26.
26.
26.
26.
26.
26.
26.
27.
27.
INTRIK(cm2) ALPHA (-)
.0000000144
.0000000149
.0000000230
.0000000211
.0000000233
.0000000245
.0000000252
.0000000253
0.0400
0.0155
0.0203
0.0192
0.0191
0.0190
0.0190
0.0189
2
2
2
2
2
2
2
2
753
813
869
928
Oil
123
257
437
708
806
450
538
574
739
260
355
431
494
578
661
780
833
914
004
405
682
875
974
090
134
334
427
583
724
980
169
218
N(
.0000
.2677
.5337
. 6716
.7513
.7906
.8079
.8124
24.
25.
25.
25.
25.
25.
25.
26.
26.
26.
25.
25.
26.
26.
25.
25.
25.
25.
26.
26.
26.
26.
26.
26.
26.
26.
26.
26.
26.
26.
26.
26.
26.
26.
26.
26.
26.
-)
987
267
371
562
729
867
983
053
112
184
710
952
086
261
740
645
769
896
012
119
260
327
376
546
275
248
325
226
340
386
496
653
737
826
820
781
982
FURTHER REDUCTION IN SSQ.
CONVERGENCE .
RSQ=
ssq=
0.998517455981560
1. 65646876697055
Correlation
1
2
3
4
1
Matrix:
1.000
0.777
-0.166
0.892
NON-LINEAR
2
1.
0.
0.
3
000
018
699
LEAST-SQUARES ANALYSIS:
4
1.000
-0.483
1.
000
FINAL RESULTS
VARIABLE
VALUE
S.E.COEFF.
95% CONFIDENCE LIMITS
LOWER UPPER
57
-------
THETAR(cm3/
INTRIK(cm2)
ALPHA(-)
N(-)
THETAS(cm3/cm3) =
THETAR(cm3/cm3)=
INTRIK(cm2)
ALPHA(-)
N(-)
M(-)
L(-)
NONE
Ksw(cm/hr)=
Ksnw(cm/hr)=
0.021 0.003
0.000 0.000
0.019 0.000
2.812 0.040
0.320000000000000
2.115814578887867E-002
= 2.531770278228916E-008
= 1.893833661150859E-002
2.81236060923382
= 0.644426821824804
= 0.500000000000000
1.00000000000000
.93208554159162
493.485389038211
0.015
0.000
0.019
2.734
0.028
0.000
0.019
2.891
FINAL
TIME
0.
0.
0.
0.
0.
0.
0.
0.
4.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
7.
7.
7.
7.
8.
8.
9.
9.
10.
11.
12.
12.
14.
16.
OBS-CAL FITTING:
OBS-HC
,017
,033
,050
,067
,083
,117
,267
,533
,167
,433
, 683
,700
,717
,733
,750
,767
,783
,800
,817
,850
,883
,917
,950
,983
,017
,067
,150
,233
,300
,383
,483
, 633
,750
,900
,117
,367
, 650
,950
,333
,750
,267
, 617
,250
,317
,267
,733
,717
,233
24.
25.
26.
26.
27.
27.
28.
28.
29.
29.
33.
34.
35.
35.
36.
36.
37.
37.
38.
38.
39.
39.
39.
40.
40.
40.
41.
41.
41.
42.
42.
42.
42.
42.
42.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
CAL-HC DFHC OBS-Q
,200
,440
,150
, 600
,000
,490
,420
,960
, 620
, 620
,450
,250
,180
,760
,380
,870
,360
,760
,160
,830
,320
,760
,760
,160
,480
,830
,370
,810
,900
,210
,430
,480
,740
,880
,970
,060
,230
,320
,320
,320
,320
,320
,320
,320
,320
,320
,320
,320
24.
25.
25.
25.
25.
26.
27.
28.
29.
29.
30.
30.
31.
32.
32.
33.
33.
34.
34.
35.
36.
36.
37.
37.
38.
39.
39.
40.
41.
41.
42.
42.
42.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
43.
,819
,115
,386
, 638
,871
,276
,514
, 663
, 675
, 656
,185
,871
,514
,113
, 671
,198
, 694
,158
, 600
,398
,120
,770
,361
,899
,388
,030
,906
, 626
,114
, 618
,097
, 646
,957
,250
,524
,703
,799
,834
,832
,808
,773
,743
, 690
, 637
, 601
,569
,506
,481
0,
0,
0,
0,
1,
1,
0,
0,
0,
0,
3,
3,
3,
3,
3,
3,
3,
3,
3,
3,
3,
2,
2,
2,
2,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
CAL-Q DFQ
. 619
.325
.764
.962
.129
.214
.906
.297
.055
.036
.265
.379
. 666
. 647
.709
. 672
. 666
. 602
.560
.432
.200
.990
.399
.261
.092
.800
.464
.184
.786
.592
.333
.166
.217
.370
.554
. 643
.569
.514
.512
.488
.453
.423
.370
.317
.281
.249
.186
.161
0.
1.
1.
1.
1.
1.
1.
2.
2.
2.
2.
2.
2.
3.
3.
3.
3.
3.
3.
4.
4.
4.
4.
5.
5.
5.
5.
6.
6.
6.
6.
6.
7.
7.
7.
7.
8.
8.
8.
8.
8.
9.
9.
9.
9.
10.
10.
10.
,830
,150
,340
,520
, 660
,800
,980
,070
,210
,300
,530
, 680
,920
,150
,300
,450
, 610
,840
,990
,220
,530
, 680
,910
,150
,370
, 600
,830
,060
,300
,520
,750
,990
,220
,450
, 600
,910
,140
,370
, 600
,830
,990
,290
,520
,750
,900
,210
,360
,520
-0.
0.
0.
0.
0.
0.
1.
1.
2.
2.
2.
3.
3.
3.
4.
4.
4.
5.
5.
6.
6.
6.
7.
7.
8.
8.
9.
9.
10.
10.
10.
11.
11.
11.
11.
11.
11.
12.
12.
11.
11.
11.
11.
11.
11.
11.
11.
11.
,116
,023
,151
,273
,386
,587
,224
,846
,417
,406
,746
,154
,543
,913
,263
,598
,917
,220
,510
,042
,529
,974
,382
,757
,099
,552
,174
, 689
,040
,403
,749
,146
,372
,584
,783
,913
,982
,007
,006
,988
,962
,940
,901
,862
,836
,813
,767
,748
0,
1,
1,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
1,
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
3,
3,
3,
3,
3,
4,
4,
4,
4,
4,
3,
3,
3,
3,
2,
2,
2,
2,
1,
1,
1,
1,
.946
.127
.189
.247
.274
.213
.756
.224
.207
.106
.216
.474
. 623
.763
.963
.148
.307
.380
.520
.822
.999
.294
.472
. 607
.729
.952
.344
. 629
.740
.883
.999
.156
.152
.134
.183
.003
.842
. 637
.406
.158
.972
. 650
.381
.112
.936
. 603
.407
.228
58
-------
19.
20.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
21.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
22.
23.
23.
23.
23.
24.
24.
25.
26.
27.
27.
27.
27.
28.
29.
30.
32.
33.
,417
, 617
,267
,283
,300
,317
,333
,350
,367
,383
,400
,417
,433
,450
,483
,500
,533
,550
,567
, 600
, 633
, 650
, 683
,700
,750
,783
,817
,850
,883
,933
,983
,033
,067
,083
,100
,117
,133
,150
,167
,183
,217
,283
,350
,417
,533
, 650
,767
,900
,050
,233
,467
,750
,117
, 650
, 600
,167
,150
,433
, 600
,900
,483
,567
,717
,017
,033
43.
43.
43.
46.
47.
48.
49.
49.
50.
51.
51.
52.
52.
52.
53.
54.
54.
54.
55.
55.
55.
56.
56.
56.
57.
57.
57.
57.
57.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
58.
59.
59.
59.
59.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
60.
,320
,320
,320
,200
,400
,420
,220
,980
,150
,090
, 620
,070
,510
,870
, 670
,030
,510
,780
,000
,490
,940
,120
,430
, 610
,050
,230
,490
, 670
,940
,160
, 610
,740
,740
,740
,740
,740
,740
,740
,740
,740
,740
,740
,740
,830
,360
,540
,720
,900
,160
,250
,340
,520
,520
,520
,520
,520
,520
,520
,520
,520
,520
, 610
, 610
, 610
, 610
43.
43.
43.
44.
45.
45.
46.
47.
47.
48.
48.
49.
49.
50.
51.
51.
52.
52.
53.
53.
54.
54.
55.
55.
56.
56.
57.
57.
57.
58.
58.
58.
59.
59.
59.
59.
59.
59.
59.
59.
59.
60.
60.
60.
60.
61.
61.
61.
61.
61.
61.
61.
61.
61.
61.
61.
61.
61.
61.
61.
61.
61.
60.
60.
60.
,435
,394
,961
, 666
,345
,997
, 622
,217
,791
,342
,867
,374
,862
,327
,198
, 610
,386
,753
,110
,772
,386
, 678
,226
,486
,192
, 627
,028
,400
,744
,205
, 615
,980
,205
,312
,414
,511
, 603
, 690
,772
,850
,994
,238
,444
, 619
,858
,035
,167
,273
,352
,406
,435
,435
,412
,368
,306
,271
,223
,196
,173
,131
,068
,003
,955
,909
,866
0,
0,
0,
1,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
1,
1,
1,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
.115
.074
. 641
.534
.055
.423
.598
.763
.359
.748
.753
. 696
. 648
.543
.472
.420
.124
.027
.890
.718
.554
.442
.204
.124
.858
. 603
.462
.270
.196
.045
.005
.240
.465
.572
. 674
.771
.863
.950
.032
.110
.254
.498
.704
.789
.498
.495
.447
.373
.192
.156
.095
.915
.892
.848
.786
.751
.703
. 676
. 653
. 611
.548
.393
.345
.299
.256
10.
10.
11.
11.
12.
12.
12.
13.
13.
13.
14.
14.
14.
14.
15.
15.
15.
16.
16.
16.
16.
17.
17.
17.
17.
17.
18.
18.
18.
18.
19.
19.
19.
19.
20.
20.
20.
20.
21.
21.
21.
21.
21.
22.
22.
22.
22.
23.
23.
23.
23.
24.
24.
24.
24.
25.
25.
25.
25.
26.
26.
26.
26.
27.
27.
,900
,980
,500
,830
,330
, 660
,990
,310
, 640
,970
,220
,460
,710
,960
,370
,530
,860
,110
,270
, 600
,850
,000
,330
,420
,830
,990
,320
,480
,730
,970
,300
,470
, 630
,880
,130
,370
,710
,950
,110
,190
,450
, 690
,930
,100
,510
, 670
,920
,160
,410
, 660
,910
,150
,400
, 640
,970
,130
,380
,700
,870
,110
,360
, 610
,850
,100
,350
11.
11.
12.
12.
13.
13.
14.
14.
15.
15.
15.
16.
16.
16.
17.
17.
18.
18.
18.
19.
19.
19.
20.
20.
20.
21.
21.
21.
21.
22.
22.
22.
22.
23.
23.
23.
23.
23.
23.
23.
23.
23.
23.
23.
23.
24.
24.
24.
24.
24.
24.
24.
24.
24.
24.
24.
24.
24.
24.
24.
24.
24.
24.
23.
23.
,715
, 685
,217
,747
,247
,727
,186
, 622
,041
,443
,824
,192
,544
,879
,500
,795
,341
, 601
,851
,311
,735
,936
,310
,487
,962
,253
,520
,765
,991
,292
,557
,792
,936
,005
,069
,131
,188
,242
,294
,343
,433
,586
,715
,824
,972
,082
,163
,228
,276
,308
,325
,325
,309
,282
,243
,221
,190
,172
,158
,131
,091
,050
,020
,990
,963
0,
0,
0,
0,
0,
1,
1,
1,
1,
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
2,
2,
2,
2,
3,
3,
3,
3,
3,
3,
3,
3,
3,
3,
3,
2,
2,
2,
2,
2,
2,
1,
1,
1,
1,
1,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
1,
1,
1,
1,
2,
2,
2,
3,
3,
.815
.705
.717
.917
.917
.067
.196
.312
.401
.473
. 604
.732
.834
.919
.130
.265
.481
.491
.581
.711
.885
.936
.980
.067
.132
.263
.200
.285
.261
.322
.257
.322
.306
.125
.939
.761
.478
.292
.184
.153
.983
.896
.785
.724
.462
.412
.243
.068
.866
. 648
.415
.175
.091
.358
.727
.909
.190
.528
.712
.979
.269
.560
.830
.110
.387
59
-------
34.
36.
37.
41.
44.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
45.
46.
46.
46.
46.
46.
47.
47.
50.
53.
53.
55.
57.
61.
64.
67.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
69.
70.
70.
71.
74.
74.
79.
80.
87.
,117
,067
, 650
,150
, 633
,317
,333
,367
,383
,400
,417
,450
,483
,517
,533
,567
, 617
, 667
,717
,767
,817
,867
,900
,950
,967
,000
,183
,417
, 667
,900
,233
,983
,867
,150
,367
,117
,533
,233
,100
,883
,033
,067
,083
,117
,150
,183
,217
,233
,267
,317
,400
,483
,533
, 617
,700
,833
,983
,217
,933
,983
,050
,850
,367
,500
, 633
60.
60.
60.
60.
60.
64.
65.
66.
67.
67.
68.
69.
70.
70.
71.
71.
72.
73.
73.
74.
74.
75.
75.
75.
75.
76.
77.
78.
78.
78.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
80.
80.
81.
82.
83.
84.
84.
85.
86.
88.
89.
90.
91.
91.
92.
93.
95.
97.
98.
98.
99.
99.
99.
99.
, 610
, 610
, 610
, 610
, 610
,370
,220
, 640
,310
,840
,370
,260
,110
,780
,130
,760
,560
,400
,980
,470
,960
,400
, 630
,940
,940
,340
,230
,120
,700
,870
,230
,410
, 630
, 630
, 670
, 670
,720
,720
,720
,720
,810
,060
,590
, 610
, 680
, 610
,460
,900
, 620
, 640
,060
,310
,020
,000
,800
,870
,850
,180
,140
,290
,870
,050
,270
,270
,270
60.
60.
60.
60.
60.
61.
61.
62.
63.
64.
64.
65.
66.
67.
68.
68.
70.
71.
71.
72.
73.
74.
74.
75.
75.
75.
76.
78.
78.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
79.
80.
81.
82.
83.
84.
84.
85.
86.
87.
89.
89.
90.
91.
93.
94.
95.
97.
98.
99.
99.
99.
98.
98.
,821
,761
,717
, 680
, 645
,016
, 651
,950
,557
,144
,715
,782
,769
, 681
,117
,936
,044
,038
,931
,735
,460
,112
,515
,064
,238
,556
,880
,012
,773
,213
,576
,879
,955
,925
,913
,853
,820
,780
,750
,715
, 689
,874
,264
,274
,360
,335
,212
, 623
,392
,430
,914
,196
,898
,925
,826
,053
,173
,493
, 681
,725
,043
,040
,025
,994
,935
0,
0,
0,
0,
0,
3,
3,
3,
3,
3,
3,
3,
3,
3,
3,
2,
2,
2,
2,
1,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
.211
.151
.107
.070
.035
.354
.569
. 690
.753
. 696
. 655
.478
.341
.099
.013
.824
.516
.362
.049
.735
.500
.288
.115
.876
.702
.784
.350
.108
.073
.343
.346
.469
.325
.295
.243
.183
.100
.060
.030
.005
.121
.186
.326
.336
.320
.275
.248
.277
.228
.210
.146
.114
.122
.075
.026
.183
.323
.313
.541
.435
.173
.010
.245
.276
.335
27.
27.
28.
28.
28.
28.
28.
29.
29.
29.
29.
29.
30.
30.
30.
30.
30.
30.
31.
31.
31.
31.
31.
31.
32.
33.
33.
33.
33.
34.
34.
34.
34.
34.
35.
35.
35.
35.
35.
35.
36.
36.
36.
37.
37.
37.
38.
38.
38.
38.
39.
39.
39.
39.
40.
40.
40.
40.
41.
41.
41.
41.
41.
42.
42.
,590
,930
,090
,340
,500
, 680
,940
,180
,370
,490
, 670
,800
,040
,280
,410
,530
,780
,960
,140
,270
,510
, 640
,700
,820
,990
,300
,480
, 660
,910
,090
,210
,390
,580
,700
,010
,130
,320
,560
, 690
,940
,000
,270
, 620
,070
,430
,780
,050
,230
,500
,760
,120
,390
, 660
,920
,090
,370
, 640
,910
,250
,430
,790
,790
,960
,220
,490
23.
23.
23.
23.
23.
24.
24.
25.
25.
26.
26.
27.
27.
28.
28.
29.
29.
30.
30.
31.
31.
31.
31.
32.
32.
32.
32.
33.
33.
33.
34.
34.
34.
34.
34.
34.
34.
34.
34.
34.
34.
34.
34.
35.
35.
36.
36.
36.
36.
37.
37.
38.
38.
38.
39.
39.
39.
40.
40.
41.
41.
41.
41.
41.
41.
,934
,896
,868
,844
,822
,326
,793
, 609
,992
,351
, 692
,308
,870
,379
, 621
,063
, 646
,159
, 611
,010
,362
, 674
,865
,121
,200
,343
,938
,435
,761
,947
,098
,224
,255
,242
,235
,210
,195
,178
,164
,149
,137
,556
,890
,405
,839
,213
,541
, 694
,971
,335
,843
,270
,499
,825
,106
,480
,811
,192
,800
,081
,165
,164
,159
,150
,133
3,
4,
4,
4,
4,
4,
4,
3,
3,
3,
2,
2,
2,
1,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
1,
1,
. 656
.034
.222
.496
. 678
.354
.147
.571
.378
.139
.978
.492
.170
.901
.789
.467
.134
.801
.529
.260
.148
.034
.165
.301
.790
.957
.542
.225
.149
.143
.112
.166
.325
.458
.775
.920
.125
.382
.526
.791
.863
.714
.730
. 665
.591
.567
.509
.536
.529
.425
.277
.120
.161
.095
.984
.890
.829
.718
.450
.349
. 625
. 626
.801
.070
.357
60
-------
92.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
93.
94.
94.
94.
95.
95.
96.
97.
101.
102.
102.
104.
117.
122.
122.
122.
123.
123.
124.
124.
125.
126.
132.
139.
146.
150.
151.
151.
152.
153.
155.
157.
162.
167.
168.
16.
,967
,000
,017
,033
,067
,100
,133
,167
,233
,317
,383
,467
,550
, 633
,733
,850
,967
,100
,300
, 600
,000
, 617
,733
,200
,517
,367
,750
,733
,117
,517
, 667
,800
,083
,483
,267
,700
,450
,417
,550
,233
,083
,900
,100
,250
, 600
,550
,550
,783
,883
,550
,917
,000
99.
99.
100.
100.
101.
102.
103.
104.
105.
107.
108.
110.
111.
112.
114.
115.
117.
118.
120.
122.
125.
129.
134.
135.
143.
144.
144.
145.
148.
148.
148.
149.
151.
154.
159.
161.
164.
168.
181.
188.
192.
193.
198.
199.
205.
210.
219.
227.
244.
257.
261.
0.
,270
,430
,100
,500
,390
,330
,210
,020
, 620
,440
,730
,200
,530
,730
,070
, 620
,000
,340
,340
,960
,990
,590
,260
,770
,380
,130
,490
,740
,000
,400
, 610
,060
,280
,300
,150
,420
,890
,580
, 610
,420
,200
,890
,950
,380
, 640
,080
,070
,790
,790
, 690
,200
,308
98.
99.
99.
100.
101.
102.
103.
104.
106.
108.
109.
111.
112.
114.
115.
117.
119.
120.
123.
126.
129.
133.
138.
140.
146.
147.
147.
148.
148.
148.
148.
149.
152.
155.
161.
163.
167.
171.
186.
193.
196.
197.
197.
197.
204.
209.
217.
225.
239.
250.
253.
0.
,910
,259
,586
,065
,093
,268
,367
,379
,219
,221
, 659
,314
,843
,273
,871
, 601
,209
,914
,233
,273
, 672
,805
,965
,510
,930
,372
,527
,046
,552
,559
, 658
,368
,075
,718
,105
,549
,254
,322
,739
,726
,766
,799
,441
,777
,953
,578
,566
,113
, 653
,848
,840
,313
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
1.
1.
1.
1.
1.
2.
2.
2.
3.
3.
4.
4.
4.
3.
3.
3.
2.
0.
0.
0.
0.
0.
1.
1.
2.
2.
2.
5.
5.
4.
3.
1.
1.
0.
0.
1.
2.
5.
6.
7.
0.
,360
,171
,514
,435
,297
,062
,157
,359
,599
,781
,929
,114
,313
,543
,801
,981
,209
,574
,893
,313
, 682
,215
,705
,740
,550
,242
,037
,306
,552
,159
,048
,308
,795
,418
,955
,129
,364
,742
,129
,306
,566
,909
,509
, 603
, 687
,502
,504
, 677
,137
,842
,360
,005
42.
42.
43.
43.
44.
44.
44.
44.
45.
45.
45.
45.
46.
46.
46.
47.
47.
47.
47.
47.
48.
48.
48.
48.
49.
49.
49.
50.
50.
50.
50.
50.
50.
50.
51.
51.
51.
51.
51.
51.
52.
52.
52.
52.
52.
53.
53.
53.
53.
53.
54.
,500
,990
,750
,990
,330
,570
,830
,990
,320
,490
,820
,990
,240
,490
,740
,080
,240
,490
,740
,990
,240
,490
,820
,990
,160
,490
, 660
,000
,000
,000
,200
,390
,590
,780
,040
,160
,290
,550
, 680
,930
,200
,200
,430
,520
,830
,080
,320
,550
,800
,950
,200
41.
41.
42.
42.
42.
43.
43.
43.
43.
44.
44.
45.
45.
45.
45.
46.
46.
46.
47.
47.
48.
48.
49.
49.
50.
50.
50.
50.
50.
50.
50.
50.
51.
51.
51.
52.
52.
52.
53.
53.
53.
54.
54.
54.
54.
54.
55.
55.
55.
56.
56.
,126
,911
,199
,439
,789
,092
,356
,589
,990
,416
,718
,051
,351
, 622
,914
,222
,498
,780
,146
, 603
,081
, 622
,242
,416
,099
,143
,158
,210
,260
,260
,541
,789
,110
,419
,844
,029
,298
,576
,490
,844
,989
,037
,290
,385
,802
,987
,287
,548
,988
,280
,351
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
2,
2,
2,
.374
.079
.551
.551
.541
.478
.474
.401
.330
.074
.102
.939
.889
.868
.826
.858
.742
.710
.594
.387
.159
.132
.422
.426
.939
. 653
.498
.210
.260
.260
.341
.399
.520
. 639
.804
.869
.008
.026
.810
.914
.789
.837
.860
.865
.972
.907
.967
.998
.188
.330
.151
61
-------
Listing of TF-OPT
c
C. PROGRAM
C. PURPOSE
C.
C.
C. DEVELOPED BY
C. BASED ON
C.
C.
C.
C.
C.
C
: TF-OPT. FOR
: INVERSE PARAMETER ESTIMATION OF TWO-PHASE FLOW
CAPILLARY PRESSURE AND PERMEABILITY FUNCTION IN
MULTI-STEP OUTFLOW EXPERIMENT OF SOIL COLUMN
: JIAYU CHEN (HYDROLOGIC SCIENCE, LAWR, UCD)
: - TPH1D FROM DR. JOHN NIEBER (UNIV. MINNESOTA)
FOR TWO-PHASE UNSATURATED FLOW MODEL
- MLSTPM (LAWR PAPER NO. 100021 OF UNIVERSITY
OF CALIFORNIA AT DAVIS) FOR LEVENBERG-MARQUARDT
OPTIMIZATION
: OCT., 1997
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
DIMENSION FC(MNOB),FC1(MNOB),FC2(MNOB),R(MNOB)
DIMENSION FREEPAR(MPAR),X1(MPAR),X2(MPAR)
DIMENSION DD(MPAR,MPAR),A(MPAR,MPAR),AS(MPAR,MPAR)
DIMENSION E(MPAR),EWJAC(MPAR),C(MPAR),CHI(MPAR)
CHARACTER*60 INFILE,OUTFILE
INCLUDE 'tfcombk.dat'
DATA ZERO/0./
C... initiation
WRITE(*,*) 'Enter the input file name:'
READ(*,'(A60)') INFILE
OPEN(UNIT=11,FILE=INFILE,STATUS='OLD')
WRITE(*,*) 'Enter the output file name:'
READ(*,'(A60)') OUTFILE
OPEN(UNIT=12,FILE=OUTFILE,STATUS='UNKNOWN')
CALL INPUT(FREEPAR)
DO 4 I=1,NPAR
XI(I)=FREEPAR(I)
X2 (I)=X1(I)
4 CONTINUE
C... first model call
NIT = 0
WRITE(*,*) 'NIT=',NIT
CALL MODEL(FC)
IF (ABS(MODE) .eq. 1) THEN
CALL SIMULATION_REPORT(FC)
GO TO 8888
ENDIF
CALL OPTIMIZATION_INITIAL_REPORT(FC)
C.. check data fitness
SSQ=0.
DO 10 I = 1,NOB
R(I) = WGHT(I)*(FO(I)-FC(I))
SSQ = SSQ+R(I)*R(I)
10 CONTINUE
sssq=SQRT(SSQ/NOB)
IF (NPAR .NE. 0) THEN
62
-------
WRITE(12,1040) (PNAMA(I),I = 1,NPAR)
WRITE(*,1040) (PNAMA(I),I = 1,NPAR)
WRITE(12,1042) NIT,sssq,(XI(I),I = 1,NPAR)
WRITE(*,1042) NIT,sssq,(XI(I),I = 1,NPAR)
ELSE
WRITE(12,1040)
WRITE(12,1042) NIT,sssq
ENDIF
IF (MIT .EQ. 0 .OR. NPAR .EQ. 0) GO TO 8000
1040 FORMAT(/T2,'ITERATION',3X,'SSQ',2X,A16,IX,7(A12,IX)
1042 FORMAT(4X,I3,2X,E12.4,3X,F8.4,1X,F12.10,6(F8.4,1X))
C... optimization loop
GA = 0.02
STOPCR =0.01
DO 18 I = 1,NPAR
18 E(I) = 0.
20 NIT = NIT+1
NET = 0
GA = 0.1*GA
C. . . evaluate weighted j acobian j(i,j) = dr(i)/dxl(j) and ewj ac(j)
WRITE(*,*) 'NIT=',NIT
DO 38 J = 1,NPAR
PVALA(TRFPAR(J)) = 1.01*X1(J)
IF (XI(J) .EQ. 0.) STOP 'X1(J)=0, DEVIDED BY ZERO ERROR.'
EWJAC(J) = 0.
CALL MODEL(FC1)
DO 36 I = 1,NOB
QJAC(I,J) = WGHT(I)*(FC1(I)-FC(I))
EWJAC(J) = EWJAC(J)+QJAC(I,J)*R(I)
36 CONTINUE
EWJAC(J) = 100.*EWJAC(J)/Xl(J)
PVALA(TRFPAR(J)) = XI(J)
38 CONTINUE
DO 44 I = 1,NPAR
DO 42 J = 1,1
SUM = ZERO
DO 40 K = 1,NOB
SUM = SUM+QJAC(K,I)*QJAC(K,J)
40 CONTINUE
DD(I,J) = 10000.*SUM/(X1(I)*X1(J))
DD(J,I) = DD(I,J)
42 CONTINUE
SCAL=DD(1,1)
IF (SCAL .LT. l.OE-30) SCAL=1.0E-30
SCAL=DSQRT(SCAL)
IF (E(I) .LT. SCAL) E(I)=SCAL
44 CONTINUE
50 DO 53 I = 1,NPAR
DO 52 J = 1,NPAR
A(I, J) = DD(I,J)/(E(I)*E(J) )
52 CONTINUE
53 CONTINUE
DO 54 I = 1,NPAR
C(I) = EWJAC(I)/E(I)
63
-------
CHI(I) = C(I)
A (I,I) = A(I,I)+GA
54 CONTINUE
CALL QRSOLV(A,NPAR,C)
STEP = 1.
56 NET = NET+1
IF ( NET .GE. MAXTRY ) THEN
WRITE(*,*) 'MEET MAXIMUM NET, NO FURTHER REDUCTION IN SSQ.'
WRITE(12,*) 'MEET MAXIMUM NET, NO FURTHER REDUCTION IN SSQ.
GO TO 96
ENDIF
DO 58 I = 1,NPAR
X2(I) = C(I)*STEP/E(I)+X1(I)
IF (X2(I) .LT. PMINA(I)) X2(I)=PMINA(I)
IF (X2(I) .GT. PMAXA(I)) X2(I)=PMAXA(I)
C(I) = (X2(I)-XI(I))*E(I)/STEP
PVALA(TRFPAR(I)) = X2(I)
58 CONTINUE
SUM1 = ZERO
SUM2 = ZERO
SUMS = ZERO
DO 62 I = 1,NPAR
SUM1 = SUM1+C(I)*CHI(I)
SUM2=SUM2+C(I)*C(I)
SUM3=SUM3+CHI(I)*CHI(I)
62 CONTINUE
DUM=SUM2*SUM3
IF (DUM .EQ. 0.) DUM=0.000000001
ARG=SUM1/DSQRT(SUM2*SUM3)
ANGLE=57.29578*DATAN2(DSQRT(1.-ARG*ARG),ARG)
DO 64 I=1,NPAR
IF ( X1(I)*X2(I) .LE. 0.) GO TO 70
64 CONTINUE
SUMB = ZERO
CALL MODEL(FC2)
DO 66 I = 1,NOB
R(I) = WGHT(I)*(FO(I)-FC2(I))
66 SUMB = SUMB+R(I)*R(I)
C SUMB=SQRT(SUMB/NOB)
IF (NET .GE. MAXTRY) THEN
WRITE(*,*) 'NO FURTHER REDUCTION IN SSQ.'
WRITE(12,*) 'NO FURTHER REDUCTION IN SSQ.'
GO TO 96
ENDIF
IF ( SUMB/SSQ-1. GT. 0.) GO TO 70
IF ( SUMB/SSQ-1. LE. 0.) GO TO 80
70 DO 75 I = 1,NPAR
PVALA(TRFPAR(I)) =X1(I)
75 CONTINUE
IF ( ANGLE-30.0 .GT. 0.) THEN
GA=10.*GA
IF (GA .GT. 100.) GA=100.
GO TO 50
64
-------
ELSE
STEP = 0.5*STEP
GO TO 56
ENDIF
80 ssum=sqrt(sumb/nob)
WRITE(*,1042) NIT,ssum,(X2(J),J=1,NPAR)
WRITE(12,1042) NIT,ssum,(X2(J),J=1,NPAR)
DO 90 I = 1,NPAR
DUM = ABS(C(I)*STEP/E(I))/(1.OE-20+ABS(X2(I)))-STOPCR
IF (DUM .GT. 0.) THEN
DO 82 J = 1,NPAR
82 XI(J) = X2(J)
DO 83 J = 1,NOB
83 FC(J) = FC2(J)
SSQ=SUMB
IF (NIT .LT. MIT) THEN
GO TO 20
ELSE
WRITE(*,*) 'MEET MAXIMUM NIT, NO FURTHER REDUCE IN SSQ.
GO TO 96
ENDIF
ENDIF
90 CONTINUE
C END OF ITERATION LOOP
96 CONTINUE
WRITE(*,*) 'CONVERGENCE.'
WRITE(12,*) 'CONVERGENCE.'
CALL MATINV(DD,NPAR)
C WRITE RSQUARE, CORRELATION MATRIX
SUMS=SUMB
ssum=sqrt(SUMB/NOB)
SUMS1=0.0
SUMS2=0.0
DO 98 1=1,NOB
FOS=FO(I)
SUMS1=SUMS1+FOS
SUMS2=SUMS2+FOS*FOS
98 CONTINUE
RSQ= l.-SUMS/(SUMS2-SUMS1*SUMS1/NOB)
WRITE(*,*) 'RSQ=',RSQ
WRITE(12,*) 'RSQ=',RSQ
write(*,*) 'ssq=',ssum
write(12,*) 'ssq=',ssum
DO 100 I=1,NPAR
IF (E(I) .LT. l.OE-30) E(I)=1.0E-30
E(I)=DSQRT(E(I))
100 CONTINUE
WRITE(*,*) 'Correlation Matrix:'
WRITE(12,*) 'Correlation Matrix:'
WRITE(*,*) (I,I=1,NPAR)
WRITE(12,*) (I,I=1,NPAR)
DO 104 I=1,NPAR
DO 102 J=1,I
AS (J,I)=DD(J,I)/(E(I)*E(J))
65
-------
102 CONTINUE
WRITE(*,105) I,(AS(J,I),J=l,I)
WRITE(12,105) I,(AS(J,I),J=1,I)
104 CONTINUE
105 FORMAT(IX,16,IX,4(F14.3,IX))
C CALCULATE 95% CONFIDENCE INTERVAL
106 ZZ = 1./FLOAT(NOB-NPAR)
SDEV = DSQRT(ZZ*SUMS)
WRITE(*,1052)
WRITE(12,1052)
1052 FORMAT(//11X,'NON-LINEAR LEAST-SQUARES ANALYSIS: FINAL RESULTS'/
111X,48(1H=)/53X,'95% CONFIDENCE LIMITS'/11X,'VARIABLE',8X,'VALUE',
27X,'S.E.COEFF.',4X,'LOWER',8X,'UPPER')
TVAR=1.96+ZZ*(2.3779+ZZ*(2.7135+ZZ*(3.187936+2.466666*ZZ**2)))
DO 108 I=1,NPAR
SECOEF=SNGL(E(I))*SDEV
TSEC=TVAR*SECOEF
TMCOE=X2(I)-TSEC
TPCOE=X2(I)+TSEC
WRITE(*,109) PNAMA(I),X2(I),SECOEF,TMCOE,TPCOE
WRITE(12,109) PNAMA(I),X2(I),SECOEF,TMCOE,TPCOE
108 CONTINUE
109 FORMAT(1X,A12,IX,4(F14.3,IX))
CALL OPTIMIZATION_FINAL_REPORT(FC2)
8000 CLOSE(UNIT=11)
CLOSE(UNIT=12)
8888 STOP 'OKAY'
END
c
SUBROUTINE OPTIMIZATION_INITIAL_REPORT(FC)
C
C PURPOSE: TO REPORT THE MATCH OF OBSERVATION AND INITIAL-GUESS
C CALCULATION
C
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
DIMENSION FC(MNOB)
INCLUDE 'tfcombk.dat'
WRITE(12,*) 'INITIAL OBS-CAL FITTING:'
WRITE(12,*) 'TIME OBS-HC CAL-HC DFHC OBS-Q CAL-Q DFQ'
DO 6 I=1,NTOB
NH=2*I-1
NQ=2*I
DFHC=abs(FO(NH)-FC(NH))
DFQ=abs(FO(NQ)-FC(NQ))
WRITE(12,7) FTIME(I),FO(NH),FC(NH),DFHC,FO(NQ),FC(NQ),DFQ
6 CONTINUE
I=NTOB+1
NTH=2*I-1
DFTH=ABS(FO(NTH)-FC(NTH))
WRITE(12,8) FTIME(I),FO(NTH),FC(NTH),DFTH
7 FORMAT(IX,7(F9.3,IX))
8 FORMAT(1X,4(F9.3,1X))
WRITE(12,*)
66
-------
RETURN
END
SUBROUTINE OPTIMIZATION_FINAL_REPORT(FC)
C
C PURPOSE: TO REPORT THE OPTIMIZATION RESULTS
C
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
DIMENSION FC(MNOB)
INCLUDE 'tfcombk.dat'
C PREPARE FINAL OUTPUT
IF (IEQ .EQ. 1) PVALA(6) = 1-1/PVALA(5)
DO 2 I=1,MPAR
WRITE(*,*) PNAMI(I),'=',PVALA(I)
WRITE(12,*) PNAMI(I),'=',PVALA(I)
2 CONTINUE
CKSW=(3600*980)*PVALA(3)/CMUW
CKSNW=CKSW*RATIOK
WRITE(12,*) 'Ksw (cm/hr) =',CKSW
WRITE(12,*) 'Ksnw (cm/hr)=',CKSNW
WRITE(12,*)
WRITE(12,*) 'FINAL OBS-CAL FITTING:1
WRITE(12,*) 'TIME OBS-HC CAL-HC DFHC OBS-Q CAL-Q DFQ'
DO 6 I=1,NTOB
NH=2*I-1
NQ=2*I
DFHC=abs(FO(NH)-FC(NH))
DFQ=abs(FO(NQ)-FC(NQ))
WRITE(12,7) FTIME(I),FO(NH),FC(NH),DFHC,FO(NQ),FC(NQ),DFQ
6 CONTINUE
I=NTOB+1
NTH=2*I-1
DFTH=ABS(FO(NTH)-FC(NTH))
WRITE(12,8) FTIME(I),FO(NTH),FC(NTH),DFTH
7 FORMAT(IX,7(F9.3,IX))
8 FORMAT(1X,4(F9.3,1X))
RETURN
END
c
SUBROUTINE SIMULATION_REPORT(FC)
C
C PURPOSE: TO REPORT THE SIMULATION RESULTS
C
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
DIMENSION FC(MNOB)
INCLUDE 'tfcombk.dat'
WRITE(12,*) 'Observation depth Z=',Z(iobsnode)
WRITE(12,*) ' Time he ohc Qw OQw Sw Snw'
WRITE(*,*) 'Observation depth Z=',Z(iobsnode)
WRITE(*,*) ' Time he ohc Qw OQw Sw Snw'
DO 5 I=1,NTOB
NH=2*I-1
67
-------
NQ=2*I
HA=FC (NH)
HW=0.
CALL COEFT (HW, HA, CKW, CKA, CWW, CAA, CWA, SW, SA)
WRITE (12, 7) FTIME (I) , FC (NH) , FO (NH) , FC (NQ) , FO (NQ) , SW, SA
WRITE (*, 7) FTIME ( I ) , FC (NH) , FO (NH) , FC (NQ) , FO (NQ) , SW, SA
5 CONTINUE
7 FORMAT (IX, 7 (F9. 3, IX) )
RETURN
END
c ----------------------------------------------------------------------
SUBROUTINE QRSOLV (A, NP, B)
C
C PURPOSE: TO SOLVE LINEAR SYSTEM A*X=B BY QR-DECOMPOSITION
C WHERE A IS J'*J , B IS J'*R, AND ' DENOTES TRANSPOSE.
C THE SOLUTION X IS THE PARAMETER CORRECTION
C
PARAMETER (MNOB=5 0 0 , MPAR=8 , MTYP=5 , NDIM=1 0 0 )
IMPLICIT REAL*8 (A-H,0-Z)
DIMENSION A(MPAR,MPAR) , B (MPAR) ,A1 (MPAR) ,A2 (MPAR)
C
C REDUCE A TO UPPER TRIANGULAR FORM BY HOUSEHOLDER TRANSFORMATIONS
C ----------
IF(NP.EQ.l) THEN
B (NP) =B (NP) /A(NP,NP)
GO TO 300
ENDIF
NR=NP-1
DO 200 K=1,NR
IF (A(K,K) .EQ. 0.) THEN
A1(K) = 0.
GO TO 200
ENDIF
SS = 0.
DO 20 I = K,NP
SS = SS+A(I,K) *A(I,K)
20 CONTINUE
SIGM = DSQRT (SS)
IF (A(K,K) . LT. 0.) SIGM = -SIGM
A(K,K) = A(K,K)+SIGM
TERM = SIGM*A(K,K)
Al (K) = TERM
A2 (K) = -SIGM
DO 100 J = K+1,NP
SS = 0.0
DO 80 I = K,NP
SS = SS+A(I,K) *A(I, J)
80 CONTINUE
SS = SS/TERM
DO 90 I = K,NP
A(I,J) = A(I, J)-SS*A(I,K)
90 CONTINUE
100 CONTINUE
200 CONTINUE
A2 (NP) = A(NP,NP)
C
68
-------
C APPLY TRANSFORMATIONS TO B
DO 230 J = 1,NR
SS = 0.
DO 210 I = J,NP
SS = SS+A(I,J)*B(I)
210 CONTINUE
SS = SS/A1(J)
DO 220 I = J,NP
B(I) = B(I)-SS*A(I,J)
220 CONTINUE
230 CONTINUE
C SOLVE TRIANGULAR SYSTEM
B(NP) = B(NP)/A2(NP)
DO 260 I = NR,1,-1
SS = 0.
DO 250 J=I+1,NP
SS = SS+A(I,J)*B(J)
250 CONTINUE
B(I) = (B(I)-SS)/A2(I)
260 CONTINUE
C DONE, SOLUTION IS RETURNED IN B
300 RETURN
END
c
SUBROUTINE MATINV(A,NP)
C
C PURPOSE : TO INVERT J'*J
C
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
DIMENSION A(MPAR,MPAR),TINDX(6,2)
DO 2 J=l,6
2 TINDX(J,1)=0
1 = 0
4 AMAX=-1. ODO
DO 12 J=1,NP
IF(TINDX(J,1).NE.0.0) GO TO 12
6 DO 10 K=1,NP
IF(TINDX(K,1).NE.0.0) GO TO 10
8 P=DABS(A(J,K))
IF(P.LE.AMAX) GO TO 10
IR=J
IC=K
AMAX=P
10 CONTINUE
12 CONTINUE
IF(AMAX) 30,30,14
14 TINDX(IC,1)=IR
IF(IR.EQ.IC) GO TO 18
DO 16 L=1,NP
P=A(IR,L)
A(IR,L)=A(IC,L)
16 A(IC,L)=P
1 = 1 + 1
TINDX(1,2)=IC
69
-------
18 P=1./A(IC,IC)
A(IC,1C)=1.
DO 20 L=1,NP
20 A(IC,L)=A(IC,L)*P
DO 24 K=1,NP
IF(K.EQ.IC) GO TO 24
P=A(K,IC)
A(K,1C)=0.0
DO 22 L=1,NP
22 A(K,L)=A(K,L)-A(IC,L) *P
24 CONTINUE
GO TO 4
26 IC=TINDX(I,2)
IR=TINDX(1C, 1)
DO 28 K=1,NP
P=A(K, IR)
A(K,IR)=A(K,IC)
28 A(K,1C)=P
1 = 1-1
30 IF(I) 26,32,26
32 RETURN
END
SUBROUTINE INPUT(FREEPAR)
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=100)
IMPLICIT REAL*8 (A-H,0-Z)
INCLUDE 'tfcombk.dat'
DIMENSION FREEPAR(MPAR),NOBS(MNOB),SUMFO(MNOB),AVFO(MTYP)
DIMENSION SETWT(MTYP)
READ(11,*)
PART I: MODELLING DATA
READ(11,*)
READ i
'11 *
* -L -1 r
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*
READ(11,*)
DO 2 I=1,ITIM
READ(11,*) TIM(I),PB(I
CONTINUE
RATI0K=CMUW/CMUNW
NNP,IFNODE,IOBSNODE,SLENTH,PLTLENTH,DIAM
(PLTPROP(I),1=1,2)
RHOW,RHONW,CMUW,CMUNW
TMAX,DTMAX,IDT,DTI, ERR
IEQ
MODE
hnwO,HBO
ITIM
70
-------
AREA=3.14156*DIAM**2/4
nl=IFNODE-l
n2=NNP-l
dzl=PLTLENTH/nl
dz2=SLENTH/(NNP-IFNODE)
do 10 i=l,nl
10 z(i)=(i-l)*dzl
do 20 i=nl,n2
20 z(i+l)=PLTLENTH+(i-IFNODE+l)*dz2
C.. PART II: OPTIMIZATION DATA
READ(11,*)
READ(11,*)
READ(11,*) MAXTRY,MIT
READ(11,*)
READ(11,*) NTOB,IPAR,ITYP
NOB=2*NTOB+1
NOA=NOB-1
READ(11,*)
READ(11,*)
READ(11,*) (SETWT(I),1=1,ITYP)
READ(11,*)
DO 29 1=1,IPAR
READ(11,'(A16)') PNAMI(I)
29 CONTINUE
READ(11,*)
READ(11,*)
DO 30 1=1,IPAR
READ(11,*) PVALI(I),PMINI(I),PMAXI(I),IOPT(I)
30 CONTINUE
READ(11,*)
READ(11,*)
DO 40 1=1,NTOB
READ(11,*) DUMTIME,HC,QW,HB(I)
FTIME(I)=DUMTIME
FO(2*I-1)=HC
FO(2 *I)=QW
IDATTYP(2*1-1)=1
WGHT(2*1-1)=SETWT(1)
IDATTYP(2*1)=2
WGHT(2*1)=SETWT(2)
40 CONTINUE
READ(11,*) FTIME(NTOB+1),FO(NTOB*2+1)
IDATTYP(NTOB*2+1)=3
WGHT(NTOB*2+1)=SETWT(3)
C... wls : adjustment weights according to IDATTYP
DO 50 I = 1,ITYP
SUMFO(I) = 0.
NOBS (I) = 0
50 CONTINUE
DO 60 I = 1,NOB
SUMFO(IDATTYP(I)) = SUMFO(IDATTYP(I))+FO(I)
NOBS(IDATTYP(I)) = NOBS(IDATTYP(I))+1
60 CONTINUE
DO 70 I = 1,ITYP
IF (NOBS(I) .GT. 0) AVFO(I) = SUMFO(I)/REAL(NOBS(I)
71
-------
70 CONTINUE
DO 72 I = 1,ITYP
SETWT(I) = SETWT(I) * DABS( DBLE(AVFO(2) / AVFO(I)) )
72 CONTINUE
DO 80 I = 1,NOB
WGHT(I) =WGHT(I) * DABS(DBLE(AVFO(2) /AVFO(IDATTYP(I))))
80 CONTINUE
C... rearrange parameter array
NPAR = 0
DO 90 I = 1,1PAR
PVALA(I) = PVALI(I)
IF (IOPT(I) .EQ. 1) THEN
NPAR = NPAR+1
PNAMA(NPAR) = PNAMI(I)
FREEPAR(NPAR) = PVALI(I)
PMINA(NPAR) = PMINI(I)
PMAXA(NPAR) = PMAXI(I)
TRFPAR(NPAR) = I
ENDIF
90 CONTINUE
RETURN
END
c
SUBROUTINE MODEL(FC)
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5, NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
INCLUDE 'tfcombk.dat'
DIMENSION FC(MNOB)
DIMENSION A(2*NDIM,2*NDIM),GG(2*NDIM,2*NDIM)
DIMENSION G(2*NDIM),PHI(2*NDIM),PHINEW(2*NDIM)
DIMENSION SW(NDIM),SA(NDIM)
DATA CUMBW,CUMSA,cumq/0.,0.,0./
DATA IDT,DTI /I, O.I/
C... INITIALIZATION
WRITE(*,*) 'START SIMULATION NOW.'
CALL INITIATION(PHI,PHINEW)
IF (MODE .EQ. 1) THEN
CALL VOLUME(PHI,VSW1,VSA1)
VW=VSW1*AREA
WRITE(12,*) 'INI. SOIL WATER VOLUME (POROSITY=',pvala(1),'):',VW
ENDIF
CALL VOLUME(PHI,VSW1,VSA1)
VW=VSW1*AREA
TIME=0.
NUM=0
NTRY=0
cumq=0.
NP2=NNP*2
C. . . . START A TIME STEP ,
50 CONTINUE
C... SET UP THE GLOBAL MATRIX
CALL EQ(A,G,GG,PHINEW,PHi;
72
-------
C... DIRECT EQ SOLVER
CALL BDEQSOL(A,G,NP2)
C... CHECK CONVERGENCE OF SOLUTION: PICARD ITERATION
NUM=NUM+1
CALL CONV(G,PHINEW,DMAXW,DMAXA)
IF (DMAXW.GT.ERR .OR. DMAXA.GT.ERR) THEN
IF (NUM.EQ.200) GO TO 63
IF (NUM.LT.200) GO TO 50
63 DT=0.75*DT
DT1=DT
NUM=0
GO TO 50
ENDIF
C.... A SUCCESS TIME STEP DT
TIME=TIME+DT
C... DETERMINE THE FLUX OF WATER AND AIR AT THE BOTTOM AND TOP BOUNDARIES
CALL FLUX(G,GG)
C... CALCULATE SATURATIONS
NP1=2*NNP-1
N=0
DO 100 I=1,NP1,2
N=N+1
CPW=G(I)-Z(N)*RHOW
CPA=G(I+1)-RHONW*Z(N)
CALL COEFT(CPW,CPA,CKW,CKA,CWW,CAA,CWA,SW(N),SA(N))
100 CONTINUE
C... STORAGE CHANGE & MASS BALANCE ERROR
CALL VOLUME(G,VSW,VSA)
DQW=QBW
DQA=QSA
DVW=VSW-VSW1
DVA=VSA-VSA1
cumq=cumq-DVW*AREA
ERRW=ABS(ABS(DQW)-ABS(DVW))
ERRA=ABS(ABS(DQA)-ABS(DVA))
C... OUTPUT RESULTS FOR THIS TIME STEP
CALL OUTPUT(G, FC)
C... UPDATE FOR NEXT TIME STEP
99 IF (TIME.GE.TMAX) GO TO 1000
VSW1=VSW
VSA1=VSA
DO 156 1=1,NNP
PHI(2*I)=PHINEW(2*I)
PHI(2*1-1)=PHINEW(2*1-1)
156 CONTINUE
C... CONTROLABLE TIME STEP DT
NUM=0
DTOLD=DT
CALL NEWDT
73
-------
IF (DT.EQ.0.0) THEN
WRITE(*,*) 'TIME STEP IS ZERO'
GO TO 1000
ENDIF
GO TO 50
1000 RETURN
END
c
SUBROUTINE INITIATION(PHI,PHINEW)
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
INCLUDE 'tfcombk.dat'
DIMENSION PHINEW(2*NDIM),PHI(2*NDIM)
CUMBW=0.
CUMSA=0.
NTIM=1
NTOUT=1
TIMOUT=FTIME(NTOUT)
DT=DT1
DTC=DT
HBB=(HBO+HB(1))/2
IF (MODE .GT. 0) THEN
DO 5 1=1,NNP
PHI(2*1-1)=HBB
PHI(2*I)=hnwO
HW=PHI(2*1-1)-RHOW*Z(I)
HA=PHI(2*I)-RHONW*Z(I)
IF(HA .LE. HW) stop 'Pc<=0 Error —> Check the I.e.!'
5 CONTINUE
ENDIF
IF (MODE .LT. 0) THEN
HW_TOP=hnwO-RHONW*Z(NNP)
DO 6 1=1,NNP
PHI(2*1-1)=HBB-HW_TOP
PHI(2*I)=RHONW*Z(NNP)
HW=PHI(2*1-1)-RHOW*Z(I)
HA=PHI(2*I)-RHONW*Z(I)
IF(HA .LE. HW) stop 'Pc<=o Error —> Check the I.e.!'
6 CONTINUE
ENDIF
NP2=NNP*2
DO 10 1=1,NP2
10 PHINEW(I)=PHI(I)
RETURN
END
c
SUBROUTINE OUTPUT(G,FC)
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
INCLUDE 'tfcombk.dat1
DIMENSION G(2*NDIM)
DIMENSION HW(NDIM),HA(NDIM),HTW(NDIM),HTA(NDIM)
DIMENSION FC(MNOB)
74
-------
J=0
NP1=2*NNP-1
DO 90 I=1,NP1,2
J=J+1
HW(J)=G(I)-Z(J)*RHOW
HTW(J)=G(I)
HA(J)=G(I+1)-RHONW*Z(J)
HTA(J)=G(I+1)
90 CONTINUE
IF (NTOUT .GT. NTOB) GO TO 99
IF ((TIME-0.0001).LT.TIMOUT.AND.(TIME+0.0001).GT.TIMOUT) THEN
NH=2*NTOUT-1
NQ=2*NTOUT
FC(NH)=HA(IOBSNODE)-HW(IOBSNODE)
FC(NQ)=cumq
c DFH=ABS(FC(NH)-abs(FO(NH)))
c DFQ=ABS(FC(NQ)-abs(FO(NQ)))
c WRITE(*,95) TIMOUT,FC(NH),abs(FO(NH)),DFH,IDATTYP(NH)
c WRITE(12,95) TIMOUT,FC(NH),abs(FO(NH)),DFH,IDATTYP(NH)
c WRITE(*,95) TIMOUT,FC(NQ),abs(FO(NQ)),DFQ,IDATTYP(NQ)
c WRITE(12,95) TIMOUT,FC(NQ),abs(FO(NQ)),DFQ,IDATTYP(NQ)
NTOUT=NTOUT+1
IF (NTOUT .GT. NTOB) THEN
N=IOBSNODE
HAA=FTIME(NTOUT)
HWW=0.
CALL COEFT(HWW,HAA,CKW,CKA,CWW,CAA,CWA,SWW,SAA)
NTH=2*NTOUT-1
FC(NTH) = SWW*PVALA(1)
c DF=ABS(FC(NTH)-FO(NTH))
c WRITE(*,95) FTIME(NTH),FC(NTH),FO(NTH),DF,IDATTYP(NTH)
c WRITE(12,95) FTIME(NTH),FC(NTH),FO(NTH),DF,IDATTYP(NTH)
TIMOUT = FTIME(NTOUT-1)
GO TO 99
ENDIF
TIMOUT=FTIME(NTOUT)
HBB=(HB(NTOUT-1)+HB(NTOUT))/2
ENDIF
95 FORMAT(IX,4(F9.3,IX),15)
99 RETURN
END
c
C... SUBROUTINE TO BUILD THE TWO-PHASE FLOW COEFFICIE MATRICES
C. .
SUBROUTINE EQ(A,G,GG,PHINEW,PHI)
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
INCLUDE 'tfcombk.dat'
DIMENSION GG(2*NDIM,2*NDIM),DD(2*NDIM,2*NDIM),A(2*NDIM,2*NDIM)
DIMENSION G(2*NDIM),PHINEW(2*NDIM),PHI(2*NDIM)
DIMENSION COND(4,4),CAP(4,4),DTH(2*NDIM)
C... DETERMINE HALF-TIME AND ESTIMATED NEW-TIME HYDRAULIC HEAD
NP2=NNP*2
DO 291 1=1,NP2
DO 290 J=1,NP2
GG(I,J)=0.0
75
-------
DD(I,J)=0.0
290 CONTINUE
291 CONTINUE
DO 330 1=1,4
DO 330 J=l,4
CAP(I,J)=0.0
330 COND(I,J)=0.0
NELEM=NNP-1
DO 320 N=1,NELEM
ZLENTH=Z(N+l)-Z(N)
HW1=PHINEW(2*N-1)-Z(N)*RHOW
HA1=PHINEW(2*N)-RHONW*Z(N)
HW2=PHINEW(2*(N+1)-1)-Z(N+l)*RHOW
HA2=PHINEW(2*(N+l))-RHONW*Z(N+l)
CALL COEFT(HW1,HA1,CKW1,CKA1,CWW1,CAA1,CWA1,SW1,SA1)
CALL COEFT(HW2,HA2,CKW2,CKA2,CWW2,CAA2,CWA2,SW2,SA2)
C... EVALUATE THE CAPACITANCE INTEGRAL FOR ELEME N
CAP(1,1)=ZLENTH*CWWl/2.
CAP(1,2)=ZLENTH*CWAl/2.
CAP(2,2)=ZLENTH*CAAl/2.
CAP(3,3)=ZLENTH*CWW2/2.
CAP (3,4)=ZLENTH*CWA2/2.
CAP(4,4)=ZLENTH*CAA2/2.
CAP (2,1)=CAP(1,2)
CAP (4,3)=CAP(3,4)
C... EVALUATE THE CONDUCTIVITY INTEGRAL FOR ELEME N
TKW=0.5*(CKW1+CKW2)/ZLENTH
TKA=0.5*(CKA1+CKA2)/ZLENTH
COND(1,1)=TKW
COND(1,3)=-TKW
COND(3,1)=-TKW
COND(3,3)=TKW
COND(2,2)=TKA
COND(2,4)=-TKA
COND(4,2)=-TKA
COND(4,4)=TKA
CONSTRUCT THE GLOBAL STIFFNESS MATRICES GG AND DD
DO 360 1=1,4
NROW=2*N-2+I
DO 350 J=l,4
NCOL=2*N-2+J
GG(NROW,NCOL)=GG(NROW,NCOL)+COND(I, J)
DD(NROW,NCOL)=DD(NROW,NCOL)+CAP(I,J)
350 CONTINUE
360 CONTINUE
320 CONTINUE
DO 61 1=1,NP2
DO 60 J=1,NP2
A(I,J)=GG(I,J)*DT+DD(I,J)
60 CONTINUE
61 CONTINUE
76
-------
CALL DTHET(PHINEW,PHI,DTH)
DO 80 1=1, NP2
G(I)=0.
DO 70 J=1,NP2
G(I) =G(I) +DD (I, J) * PHINEW (J)
70 CONTINUE
G(I) =G(I) -DTK (I)
80 CONTINUE
C... TOP & BOTTOM B . C . TREATMENT
DO 90 J=1,NP2
A(l, J)=0.
A(NP2, J) =0.
90 CONTINUE
A(NP2,NP2)=1.
IF (MODE . GT. 0) THEN
PHINEW ( 1 ) =HBB+Z ( 1 ) *RHOW
PHINEW(NP2) =PB (NTIM) + Z (NNP) *RHONW
ENDIF
IF (MODE . LT. 0) THEN
PHINEW ( 1 ) =HBB+Z ( 1 ) *RHOW-PB (NTIM)
PHINEW (NP2) =RHONW*Z (NNP)
DO 95 1=1, NNP
DO 93 J=1,NP2
A(2*I, J)=0.
93 CONTINUE
A(2*I,2*I)=1.
PHINEW (2*1) =RHONW*Z (NNP)
G(2*I)=PHINEW(2*I)
95 CONTINUE
ENDIF
G(l) =PHINEW(1)
G(NP2)=PHINEW(NP2)
RETURN
END
SUBROUTINE BDEQSOL(AA,BB,NN)
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
DIMENSION AA(2*NDIM,2*NDIM),BB(2*NDIM),X(2*NDIM),Y(2*NDIM)
DIMENSION ALPHA(2*NDIM),GAMA(2*NDIM),BETA(2*NDIM)
C... TRANSFOR QUINDIAGONAL MATRIX TO UPPER DIAGONAL MATRIX
1 = 1
ALPHA(I) =AA(I, I)
IF (ALPHA(I) .EQ. 0.) GO TO 400
BETA(I)=AA(I,1+1)/ALPHA(I)
GAMA(I)=AA(1,1+2)/ALPHA(I)
Y(I)=BB(I)/ALPHA(I)
1=2
ALPHA(I)=AA(I,I)-AA(I,I-1)*BETA(I-1)
IF (ALPHA(I) .EQ. 0.) GO TO 400
BETA(I)=(AA(I,I+1)-AA(I,I-1)*GAMA(I-1))/ALPHA(I)
GAMA(I)=AA(1,1+2)/ALPHA (I)
77
-------
Y(I)=(BB(I)-AA(I,I-1)*Y(I-1))/ALPHA(I)
DO 100 I=3,NN-2
A=AA(I,I-2)
B=AA(I,1-1)
C=AA(1,1)
D=AA(I,1+1)
E=AA(I,I+2)
F=BB (I)
DUM=B-A*BETA(I-2)
ALPHA(I)=C-A*GAMA(I-2)-DUM*BETA(1-1)
IF (ALPHA(I) .EQ. 0.) GO TO 400
BETA(I)=(D-DUM*GAMA(I-1))/ALPHA(I)
GAMA(I)=E/ALPHA (I)
Y(I)=(F-A*Y(I-2)-DUM*Y(I-1))/ALPHA(I)
100 CONTINUE
I=NN-1
DUM=AA(I,1-1)-AA(I,1-2)*BETA(I-2)
ALPHA(I)=AA(I,I)-AA(I,I-2)*GAMA(1-2)-DUM*BETA(1-1;
IF (ALPHA(I) .EQ. 0.) GO TO 400
BETA(I)=(AA(I,I+1)-DUM*GAMA(I-1))/ALPHA(I)
GAMA(I)=0.
Y(I)=(BB(I)-AA(I,I-2)*Y(I-2)-DUM*Y(I-l))/ALPHA(I)
I=NN
DUM=AA(I,1-1)-AA(I, 1-2)*BETA(I-2)
ALPHA(I)=AA(I,I)-AA(I,I-2)*GAMA(1-2)-DUM*BETA(1-1;
IF (ALPHA(I) .EQ. 0.) GO TO 400
BETA(I)=0.
GAMA(I)=0.
Y(I)=(BB(I)-AA(I,I-2)*Y(I-2)-DUM*Y(I-l))/ALPHA(I)
C... BACKWARD SUBSTITUTION FROM LAST ROW
X(NN)=Y(NN)
X(NN-1)=Y(NN-1)-BETA(NN-1)*X(NN)
DO 200 J=l,NN-2
I=NN-1-J
X(I)=Y(I)-BETA(I)*X(1+1)-GAMA(I)*X(1+2)
200 CONTINUE
DO 300 1=1,NN
BB(I)=X(I)
300 CONTINUE
GO TO 500
400
500
r1
WRITE (*, *)
WRITE (*,*)
WRITE (*, *)
WRITE (*,*)
STOP
RETURN
END
'Sigular EQ : '
'- In 1st iteration: Check your
'- After some iterations: Check
1 Stop here. '
initial parameter set.'
your parameter limits . '
SUBROUTINE NEWDT
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
INCLUDE 'tfcombk.dat'
78
-------
IF (IDT.EQ.O) THEN
IF (DT.LT.DT1) THEN
DT=DT1
ELSE
IF (NUM.LT.4) DT=1.04*DT
IF (NUM.GE.12) DT=DT/1.04
IF (DT.GT.DTMAX) DT=DTMAX
DT1=DT
ENDIF
ELSE
DT=DTC
ENDIF
IF (NTOUT.GT.NTOB) GO TO 101
IF ((TIME+DT) .GT. TIMOUT) DT=TIMOUT-TIME
101 IF ((TIME+DT) .GT. TMAX) DT=TMAX-TIME
IF (NTIM .LT. ITIM) THEN
I=NTIM+1
IF (TIME+DT .GT. TIM(I)) DT=TIM(I)-TIME
IF (TIME.GE.(TIM(I)-O.0001).AND.TIME.LE.(TIM(I)+0.0001)) THEN
NTIM=NTIM+1
DT=0.001
ENDIF
ENDIF
RETURN
END
c
SUBROUTINE CONV(G,PHINEW,DMAXW,DMAXA)
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=100)
IMPLICIT REAL*8(A-H,0-Z)
INCLUDE 'tfcombk.dat'
DIMENSION G(2*NDIM),PHINEW(2*NDIM)
DIMENSION SP(2),WP(2),WPS(2),WP1(2)
EP1=20.
EP2=20.
NP2=2*NNP
IIT=1
PMAX1=0.0
PMAX2=0.0
DO 90 1=1,NNP
PMAX1=MAX(PMAX1,ABS(G(2*I-1)))
90 PMAX2=MAX(PMAX2,ABS(G(2*I)))
DMAXW=0.0
DMAXA=0.0
DO 100 1=1,NNP
pl=g(2*I-l)
p2=g(2*I)
if (pl.eq.0.0) pl=l.
if (p2.eq.0.0) p2=l.
RERRW=DABS( (PHINEW(2*1-1)-G(2*1-1) ) /pi)
RERRA=DABS((PHINEW(2*1)-G(2*I))/p2)
AERRW=DABS(PHINEW(2*1-1)-G(2*1-1) )
AERRA=DABS(PHINEW(2*1)-G(2*I))
79
-------
IF (RERRW.GT.DMAXW.AND.AERRW.GT.0.001) THEN
DMAXW=RERRW
DPHIW=G(2*1-1)-PHINEW(2*1-1)
NODE1=I
ENDIF
IF (RERRA.GT.DMAXA.AND.AERRA.GT.0.001) THEN
DMAXA=RERRA
DPHIA=G(2*I)-PHINEW(2*I)
NODE2=I
ENDIF
100 CONTINUE
IF (IIT.EQ.O) THEN
DO 115 1=1,NNP
PHINEW(2*1-1)=G(2*1-1)
PHINEW(2*I)=G(2*I)
HW=G(2*1-1)-Z(I)*RHOW
HA=G(2*I)-RHONW*Z(I)
IF ((HA .LE. HW) .AND. (I .GE. IFNODE)) THEN
PHINEW(2*1-1)=HA+RHOW*Z(I)-0.01
ENDIF
115 CONTINUE
ELSE
IF (NUM.EQ.l) THEN
SP(1)=0.00
SP(2)=0.00
ELSE
IF (DPHIWO.EQ.0.00) DPHIWO=1.0
IF (DPHIAO.EQ.0.00) DPHIAO=1.0
SP(1)=DPHIW/(WP(1)*DPHIWO)
SP(2)=DPHIA/(WP(2)*DPHIAO)
ENDIF
IF (SP(1).GE.-l.OO) THEN
WPS(1)=(3.+SP(1))/(3.+DABS(SP(1)))
ELSE
WPS(!)=!/(2.*DABS(SP(1)))
ENDIF
IF (SP(2).GE.-l.OO) THEN
WPS(2)=(3.+SP(2))/(3.+DABS(SP(2)))
ELSE
WPS(2)=!/(2.*DABS(SP(2)))
ENDIF
IF (WPS(1)*DABS(DPHIW).LE.EP1) THEN
WP1(1)=WPS(1)
ELSE
WP1(1)=EP1/DABS(DPHIW)
ENDIF
IF (WPS(2)*DABS(DPHIA).LE.EP2) THEN
WP1(2)=WPS(2)
ELSE
WP1(2)=EP2/DABS(DPHIA)
80
-------
ENDIF
TMP1=G(1)
TMP2=G(NP2)
DO 103 1=1, NNP
G( 2*1-1 )=WP1 (1) * (G (2*1-1 )-PHINEW (2*1-1) ) +PHINEW (2*1-1 )
PHINEW ( 2*1-1 )=G( 2*1-1)
G(2*I) =WP1 (2) * (G (2*1) -PHINEW (2*1) ) +PHINEW(2*I)
PHINEW(2*I)=G(2*I)
HW=G( 2*1-1) -Z (I) *RHOW
HA=G(2*I)-RHONW*Z (I)
IF ((HA . LE. HW) .AND. (I . GE . IFNODE) ) THEN
PHINEW ( 2*1-1 )=HA+RHOW*Z (I) -0.01
ENDIF
103 CONTINUE
G(l) =TMP1
G(NP2)=TMP2
PHINEW (1) =G(1)
PHINEW (NP2)=G(NP2)
WP(1)=WP1 (1)
WP(2)=WP1 (2)
DPHIWO=DPHIW
DPHIAO=DPHIA
ENDIF
RETURN
END
c ----------------------------------------------------------------
SUBROUTINE COEFT (HW, HA, CKW, CKA, CWW, CAA, CWA, SW, SA)
PARAMETER (MNOB=5 0 0 , MPAR=8 , MTYP=5 , NDIM=1 0 0 )
IMPLICIT REAL*8 (A-H,0-Z)
INCLUDE 'tfcombk.dat'
C. .
C... VAN GENUCHTEN RETENSION FUNCTION
C. .
PA=3. 1415927
TS=PVALA(1)
TR=PVALA(2)
SR=TR/TS
CKSW= (3600*980) *PVALA(3) /CMUW
CKSA=CKSW*RATIOK
HAW=HA-HW
IF (N.LE.NNP .AND. N . GE . IFNODE) THEN
A=PVALA(4)
B=PVALA(5)
C=PVALA(6)
DDD1=PVALA(7)
DDD2=PVALA(8)
IF (HAW .GE. 0.) THEN
C ......... VGM_model
IF (IEQ . EQ. 1) THEN
SWE= (!./(!.+ (A*HAW) **B) ) **AM
SW=(1.-SR) *SWE+SR
SA=1. 0-SW
CWW=TS*A* (B-l. ) * (l.-SR) * (SWE** (l./AM) )
CWW=CWW* (l.-SWE** (1. /AM) ) **AM
81
-------
CAA=CWW
CWA=-CWW
CKW=CKSW*(SWE**DDD1)*(1.-(l.-SWE**(l./AM))**AM)**2.
CKA=CKSA*((1.0-SWE)**DDD1)*(l.-SWE**(1./AM))**(2.*AM)
ENDIF
VGB_model
IF (IEQ .EQ. 2) THEN
AM=(1.-2./B)
SWE=(1./(l.+(A*HAW)**B))**AM
SW=(l.-SR)*SWE+SR
SA=1.0-SW
CWW=TS*A*(B-1.)*(1.-SR)*(SWE**(1./AM))
CWW=CWW*(l.-SWE**(l./AM))**AM
CAA=CWW
CWA=-CWW
CKW=CKSW*SWE**(2.)*(1.-(l.-SWE**(1./AM))**AM)
CKA=CKSA*(1.0-SWE)**(2.)*(l.-SWE**(l./AM))**(AM)
ENDIF
BCM_model
IF (IEQ .EQ. 3) THEN
IF (HAW .LE. A) SWE=1.
IF (HAW .GT. A) SWE=(A/HAW)**B
SW=(l.-SR)*SWE+SR
SA=1.0-SW
IF (HAW .LE. A) CWW=0
IF (HAW .GT. A) CWW=(TS-TR)*B*(A**B)/HAW**(B+l)
CAA=CWW
CWA=-CWW
CKW=CKSW*SWE**(C+2+2/B)
CKA=CKSA*(1.0-SWE)**C*(l.-SWE**(1+1/B))**2
ENDIF
BCB_model
IF (IEQ .EQ. 4) THEN
IF (HAW .LE. A) SWE=1.
IF (HAW .GT. A) SWE=(A/HAW)**B
SW=(l.-SR)*SWE+SR
SA=1.0-SW
IF (HAW .LE. A) CWW=0
IF (HAW .GT. A) CWW=(TS-TR)*B*(A**B)/HAW**(B+l)
CAA=CWW
CWA=-CWW
CKW=CKSW*SWE**(3+2/B)
CKA=CKSA*(1.0-SWE)**(2)*(l.-SWE**(1+2/B))
ENDIF
BRB_model
IF (IEQ .EQ. 5) THEN
SWE=1/(1+A*HAW**B)
SW=(l.-SR)*SWE+SR
SA=1.0-SW
CWW=(TS-TR)*A*B*HAW**(B-l)*SWE**2
CAA=CWW
CWA=-CWW
CKW=CKSW*SWE**2*(1-(1-SWE)**(1-2/B) )
CKA=CKSA*(1-SWE)**(3-2/B)
ENDIF
GDM_model
IF (IEQ .EQ. 6) THEN
DUM=-0.5*A*HAW
82
-------
SWE=EXP(DUM)*(1-DUM)
SW=(l.-SR)*SWE+SR
SA=1.0-SW
CWW=(TS-TR)*1/4*A**2*HAW*EXP(DUM)
CAA=CWW
CWA=-CWW
CKW=CKSW*EXP(-A*HAW)
CKA=CKSA*(1-EXP(DUM))**2
ENDIF
C LNM_model
IF (IEQ .EQ. 7) THEN
DUM=DLOG(HAW/A)/B
X=DUM/(2**0.5)
SWE=0.5*ERFCC(X)
SW=(l.-SR)*SWE+SR
SA=1.0-SW
CWW=(TS-TR)/((2*PA)**0.5*B*HAW)
CWW=CWW*EXP(-(DLOG(HAW/A))**2/(2*B**2))
CAA=CWW
CWA=-CWW
X=(DUM+B)/(2**0.5)
CKRW=0.5*ERFCC(X)
CKW=CKSW*SWE**C*CKRW**2
CKA=CKSA*(1.0-SWE)**C*(1-CKRW)**2
ENDIF
ELSE
WRITE(*,*) 'Error! (Pc<=0) —> Stop!1
WRITE(12,*) 'Error! (Pc<=0) —> Stop!'
WRITE(*,*) 'n=',n,' ha=',ha,' hw=',hw
WRITE(*,*)'Check:ini. cond.; ini. para, guess; or para.limits'
STOP
ENDIF
ELSE IF (N.LT.IFNODE.AND.N.GE.l) THEN
CKW=PLTPROP(2)
CKA=1.0E-30
CWW=PLTPROP(3)
CAA=CWW
CWA=-CWW
SWE=1.00
SAE=0.0
SW=1.00
SA=0.00
ENDIF
200 CONTINUE
RETURN
END
c
SUBROUTINE DTHET(PHIN,PHIO,DTH)
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
INCLUDE 'tfcombk.dat'
DIMENSION PHIN(2*NDIM),PHIO(2*NDIM),DTH(2*NDIM)
DIMENSION HW2 (2) ,HW1(2) ,HA2(2) ,HA1(2) ,SW2 (2) ,SW1 (2) ,SA2(2) , SA1(2)
NELEM=NNP-1
NP2=2*NNP
DO 10 1=1,NP2
DTH(I)=0.0
83
-------
10 CONTINUE
por=PVALA(l)
DO 100 N=IFNODE,NELEM
DO 50 J=l,2
K=N+J-1
HW2(J)=PHIN(2*K-1)-Z(K)*RHOW
HW1(J)=PHIO(2*K-1)-Z(K)*RHOW
HA2(J)=PHIN(2*K)-RHONW*Z(K)
HA1(J)=PHIO(2*K)-RHONW*Z(K)
CALL COEFT(HW2(J),HA2(J),CKW,CKA,CWW,CAA,CWA,SW2(J),SA2(J))
CALL COEFT(HW1(J),HA1(J),CKW,CKA,CWW,CAA,CWA,SW1(J),SA1(J))
50 CONTINUE
ZLENTH=Z(N+l)-Z(N)
J=N
J1=J+1
DTH(2*J-1)=DTH(2*J-1)+POR*ZLENTH*(SW2(1)-SW1(1))/2.
DTK(2*J1-1)=DTH(2*J1-1)+POR*ZLENTH*(SW2(2)-SW1(2))/2.
DTK(2*J)=DTH(2*J)+POR*ZLENTH*(SA2(1)-SA1(1))/2.
DTK(2*J1)=DTH(2*J1)+POR*ZLENTH*(SA2(2)-SA1(2))/2.
100 CONTINUE
DO 200 N=1,IFNODE-1
DTK(2*N-1)=0.
DTH(2*N)=0.
200 CONTINUE
RETURN
END
c
SUBROUTINE FLUX(G,GG)
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=100)
IMPLICIT REAL*8 (A-H,0-Z)
INCLUDE 'tfcombk.dat'
DIMENSION G(2*NDIM),GG(2*NDIM,2*NDIM)
NP2=NNP*2
QBW=GG(1,1)*(G(1)-G(3))
QSA=-GG(NP2,NP2)*(G(NP2-2)-G(NP2))
CUMBW=CUMBW+QBW*DT*AREA
CUMSA=CUMSA+QSA*DT*AREA
RETURN
END
c
SUBROUTINE VOLUME(PHI,VSW,VSA)
PARAMETER (MNOB=5 0 0,MPAR=8,MTYP=5,NDIM=10 0)
IMPLICIT REAL*8 (A-H,0-Z)
INCLUDE 'tfcombk.dat1
DIMENSION PHI(2*NDIM),HW(2),HA(2),SW(2),SA(2)
POR=PVALA(1)
VSW=0.0
VSA=0.0
NELEM=NNP-1
DO 200 N=IFNODE,NELEM
ZLENTH=Z(N+1)-Z(N)
L=N-1
84
-------
DO 100 1=1,2
L=L+1
K=2*L-1
HW(I)=PHI(K)-Z(L)*RHOW
HA(I)=PHI(K+1)-RHONW*Z(L)
CALL COEFT(HW(I),HA(I),CKW,CKA,CWW,CAA,CWA,SW(I),SA(I))
VSA=VSA+POR*ZLENTH*SA(I)/2.
VSW=VSW+POR*ZLENTH*SW(I)/2.
100 CONTINUE
200 CONTINUE
RETURN
END
C. . .
C.... COMPLEMENTARY ERROR FUNCTION
FUNCTION ERFCC(X)
IMPLICIT REAL*8 (A-H,0-Z)
Z=ABS(X)
T = I./(1+0.5*Z)
ERFCC = T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(0.37409196+
* T*(.09678418+T*(-. 18628806+T*(.27886807+T*(-1.13520398+
* T*(1.48851587+T*(-.82215223+T*.17087277) ))))))))
IF (X .LT. 0.) ERFCC = 2.-ERFCC
RETURN
END
85
-------
References
Adamson, A.W. 1990. Physical chemistry of surfaces. John Wiley & Sons., Inc.
Bard, Y. 1994. Nonlinear Parameter Estimation, pp. 341. Academic Press. Orlando. FL.
Bentsen, Ramon G. 1994. Effect of hydrodynamic forces on capillary pressure and relative
permeability. Transport in Porous Media 17:121-132.
Brooke, A, D. Kendrick and A. Meeraus. 1992. General Algorithm Modeling System (GAMS),
Release 2.25; A user's guide. Page 155. Boyd & Fraser Publishing Co., MA.
Brooks, R.H. and A.T. Corey. 1964. Hydraulic properties of porous media. Colorado
Univiversity Hydrology Paper No. 3. Colorado State University. Fort Collins. CO.
Burdine, N.T. 1953. Relative permeability calculations from pore size distribution data.
Petroleum Trans., Am. Inst. Mining Eng. 198:71-77.
Busby, R.D., RJ. Lenhard, and D.E. Rolston. 1995. An investigation of saturation-capillary
pressure relations in two- and three-fluid systems for several NAPLS in different porous media.
Ground Water 33:570-578.
Carrera , J., and S.P. Neuman. 1986a. Estimation of aquifer parameters under transient and
steady state conditions, 1. Maximum likelihood incorporating prior information. Water Resour.
Research 22:199-210.
Carrera , J., and S.P. Neuman, 1986b. Estimation of aquifer parameters under transient and
steady state conditions, 2. Uniqueness, stability, and solution algorithms. Water Resour. Research
22:211-227.
Celia, M. A., E. T. Bouloutas, and R. L. Zarba. 1990, A general mass-conservative numerical
solution for the unsaturated flow equation. Water Resour. Research 26: 1483-1496.
Celia, M.A. and P. Binning. 1992. A mass conservative numerical solution for two-phase flow in
porous media with application to unsaturated flow. Water Resour. Research 28:2819-2828.
Demond, A.H. and P.V. Roberts. 1993. Estimation of two-phase relative permeability
relationships for organic liquid contaminants. Water Resour. Research 29:1081-1090.
Demond, A.H. and P.V. Roberts. 1991. Effect of interfacial forces on two-phase capillary
pressure-saturation relationships. Water Resour. Research 27:423-437.
86
-------
Eching, S.O. and J.W. Hopmans. 1993. Optimization of hydraulic functions from transient
outflow and soil water pressure data. Soil Sci. Soc. Amer. J. 57:1167-1175.
Eching, S.O., J.W. Hopmans, and O. Wendroth. 1994. Unsaturated hydraulic permeability from
transient multi-step outflow and soil water pressure data. Soil Sci. Soc. Amer. J. 58:687-695.
Finsterle, S., and Karsten Pruess. 1995. Solving the estimation-identification problem in two-
phase flow modeling. Water Resour. Research 31:913-924.
Gardner, W.R. 1956. Calculation of capillary permeability from pressure outflow data. Soil Sci.
Soc. Amer. Proc. 20:317-320.
Hopmans, J.W., T. Vogel, and P.D. Koblik. 1992. X-ray tomography of soil water distribution in
one-step outflow experiments. Soil Sci. Soc. Amer. J. 56:355-362.
Hornung, U. 1983. Identification of nonlinear soil physical parameters from an output-input
experiment, in numerical treatment of inverse problems in differential and integral equations,
edited by P. Deuflhard and E. Hairer, pp. 227-237. Birkhauser. Boston.
Kool, J.B., J.C. Parker and M.Th. van Genuchten. 1985a. ONESTEP: A nonlinear parameter
estimation program for evaluating soil hydraulic properties from one-step outflow experiments.
Virginia Agricultural Experiment Station Bulletin 85-3.
Kool, J.B., J.C. Parker, and M. Th. van Genuchten. 1985b. Determining soil hydraulic properties
from one-step outflow experiments by parameter estimation: I theory and numerical studies. Soil
Sci. Soc. Amer. J. 49:1348-1354.
Kool, J.B. and J.C. Parker. 1988. Analysis of the inverse problem for transient unsaturated flow.
Water Resour. Research 24: 817-830.
Lenhard, R.J. and J.C. Parker. 1987. Measurement and prediction of saturation-pressure
relationships in three-phase porous media systems. J. Contam. Hydrol. 1:407-424.
Leverett, M.C. 1941. Capillary behavior in porous solids. Trans. AIME 142:152-169.
Luckner, L. M.Th. van Genuchten, and D.R.Nielsen. 1989. A consistence set of parametric
models for the two-phase flow of immiscible fluids in the subsurface. Water Resour. Research 25:
2187-2193.
Marquardt, D.W. 1963. An algorithm for least-squares estimation of non-linear parameters. J.
Soc. Ind. Appl. Math. 11:431-441.
Miller, E.E. and R.D. Miller. 1956. Physical Theory for Capillary Flow Phenomena. Journal of
Applied Physics 27: 324-332.
87
-------
More, JJ. 1977. The Levenberg-Marquardt algorithm: Implementation and theory.Lecture Notes
in Mathematics 630. Edited by G.A. Watson. Springer-Verlag. New York.
Mualem, Y. 1976. A new model for predicting the hydraulic permeability of unsaturated porous
media. Water Resour. Research 12:513-522.
Parker, J.C., RJ. Lenhard, and T. Kuppusany. 1987. A parametric model for constitutive
properties governing multiphase flow in porous media. Water Resour. Research 23:618:624.
Press, W.H., S.A. Teukolsky, W.T. Vetterling and B.P. Flannery. 1992. Numerical Recipes in
C, The Art of Scientific Computing. Second Edition, Cambridge Univ. Press.
Richards, L.A. 1931. Capillary conduction of Liquids in porous medium. Physics 1:318-333.
Russo, David, E. Bresler, U. Shani, and J.C. Parker, 1991. Analysis of infiltration events in
relation to determining soil hydraulic properties by inverse problem methodology. Water Resour.
Research 27:1361-1373.
Schroth, M.H., J.D. Istok, SJ.Ahearn, and J.S. Selker. 1995. Geometry and position of light
nonaqueous-phase liquid lenses in water-wetted porous media. J. Contam. Hydrol. 19:269-287.
Toorman, A. F., PJ. Wierenga, and R.G. Flills. 1992. Parameter estimation of hydraulic
properties from one-step outflow data. Water Resour. Research 28: 3021-3028.
van Dam, J.C., J.N.M. Strieker and P. Droogers. 1994. Inverse method to determine soil
hydraulic functions from multistep outflow experiments. Soil Sci. Soc. Amer. J. 58:647-652.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of
unsaturated soils. Soil Sci. Soc. Amer. J. 44:892-898.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of
unsaturated soils. Soil Sci. Soc. Amer. J. 44:892-898.
Vemuri, V. and WJ. Korplus. 1981. Digital computer treatment of partial differential equations,
Prentice-Hall Inc. Englewood Cliffs, New Jersey 07632.
Whitaker, S.A. 1986. Flow in porous media II: the governing equation for immiscible, two-phase
flow. Transport in Porous Media 1: 105-125.
Whitaker, S.A. 1994. The closure problem for two-phase flow in homogeneous porous media.
Chemical Engineering Science 49: 765-780.
-------
Wraith, J.M., and D. Or. 1997. Nonlinear parameter estimation using spreadsheet software. J.
Natural Resource and Life Education. In Press.
Zachmann, D.W., P.C.Duchateau, and A.Klute. 1981. The calibration of the Richards' flow
equation for a draining column by parameter identification. Soil Sci. Soc. Amer. J. 45:1012-
1015.
Zachmann, D.W., P.C.Duchateau, and A.Klute. 1982. Simultaneous approximation of water
capacity and soil hydraulic conductivity by parameter identification. Soil Sci. 134:157-163.
89
------- |