Simulation
of
~onedimensional
Water
flow',
including
temperature
and
. hysteresis
effects
Agronomy and Soils D)epartmental
Series No. 10o2
Alabama Agricultural
Experiment Station
Auburn LI niversity Auburn
L'niversity. Alabama
Gale A. Buchanan,
D)irector September
1985
4,0
We
Simulation of Onedimensional Water Flow,
Including Temperature and Hysteresis Effects
J.W. Hopmans and J.H. Dane
Graduate Research Assistant and Associate Professor
of Agronomy and Soils
TABLE OF CONTENTS
Page
SUMMARY .. 4
INTRODUCTION 0**....... ...... 5
THEORET ICAL............................7
Temperature dependent hydraulic properties......7
Hysteresis.......................9
I. Universal independent domain model.......9
IIe A dependent domain model...................16
III. Extended similarity hypothesis.........19
IV. Modified dependent domain theory........20
METHODS 0 ......................
Hysteresis model.................
Water flow model.. . .. .
RESULTS AND DISCUSSION 00.............
Hysteresis model.0.0.......0.............
Water flow model..........................
LITERATURE CITED..........................
Appendix A...
Appendix B .. . . . . . . . .
22
22
23
26
26
29
32
34
40
Information contained herein is available to all regardless
of race, color, sex, or national origin.
0 0 00 00 0 0 0
0 00 00
LaritY
LIST OF FIGURES AND TABLES
Page
Figure 1. Schematic diagram of a pore.10
Figure 2. The filled pore diagrams in the Fp
plane (shadowed domain) for (a) the main wetting
process, (b) the main drying process, and (c) for
the primary drying scanning curve.13
Figure 3. Hysteretic curves at 20
0
C for the sandy
soil used in the computer simulations............24
Figure 4. Wetting scanning curves for the sandy soil
predicted by Mualem's modified dependent domain
theory.27
Figure 5. Primary loops for the sandy soil predicted by
Mualem's modified dependent domain theory.28
Figure 6. Volumetric water content (0) profiles at
temperatures of 20 (solid lines) and 40
0
C
(dashed lines) after 1 h of infiltration with a
pressure head (h) top boundary condition ........ 31
Appendix Figure Bl. Flow chart of predictorcorrector
model....................................42
Appendix Table Al. Listing of hysteresis model.36
Appendix Table Bl. Definition of main program variables
and subroutines.43
Appendix Table B2. Required input data of model.47
Appendix Table B3. Listing of simulation model.49
SUMMARY
The pressure head form of the general water flow equa
tion was numerically solved using the predictorcorrector
method (2). The model accounts for temperature effects on
the soil hydraulic properties. Hysteresis was considered by
implementing Mualems modified dependent domain theory. Scan
ning curves predicted by the model were compared with those
predicted by Mualem (8). The effects of temperature and
hysteresis on soil water movement were investigated for var
ious boundary conditions. Results indicated that these
effects depend on the type of surface boundary condition
applied and that hysteresis tends to dominate temperature
effects. A description and listing of the model is included
in the Appendix.
5
INTRODUCTION
A previous Departmental Series (2), reported investiga
tion of the effect of temperature dependent hydraulic prop
erties on soil water movement for a variety of initial and
boundary conditions. Although not considered, it was antic
ipated that hysteresis in the soil water content  pressure
head relationship (8(h)) would influence the water content
distributions, especially when wetting and drying cycles
occur. As a result of hysteresis, the 0(h) function deter
mined during drainage is not the same as the one found upon
rewetting. Furthermore, the relation between 0 and h will
depend on the volumetric water content at which the reversal
from drainage to wetting (or wetting to drainage) occurs.
Therefore, instead of a single valued isothermal 0(h) func
tion, we actually deal with a multivalued hysteresis func
tional (4).
A first attempt
in modeling
soilwater hysteresis,
using the independent domain concept, was carried out by
Poulovassilis (10). Inherent to the independent domain con
cept are two assumptions: (1) the pore space is made up of
pores or domains, with each pore size defined by two pres
sure head values, one at which the pore drains and one at
which the pore fills, i.e., the draining and filling of each
pore takes place independent of the state of the remaining
pores in the system, and (2) the water volume difference
between the drained and the filled status of each pore is
6
independent of the pressure
head.
The independent domain concept was
tested by Topp and
Miller (15), Talsma (12) and Topp
(13). Their conclusion was
that, in general, the theory predicted
scanning curves mod
erately well, except possibly at
the high water content
ranges. This failure was associated
with airentry values of
different pore sizes on the main
drying curve. Poulovassilis
and Childs (11) and Topp (14) resolved
this inadequacy by
including a domain dependence factor
so that the draining
and wetting of each pore were assumed to be dependent
on the
state of neighboring pores.
A computational scheme based
on the independent domain
model was introduced by Mualem (5,6).
Mualem and Dagan (9)
and Mualem (7) introduced a domain dependence
factor, Pd(e),
to account for the impedance of air
entry by water, which is
especially important for soils having
high airentry values.
Mualem (8) defined Pd(e)
as the ratio between the volume
of
pores actually emptied and the
volume which could have been
emptied had all pores guaranteed access
to air from neigh
boring pores.
After treating the temperature
aspects of soil hydrau
lic properties and a review of Mualem's hysteresis theory,
the hysteresis model of Mualem (8) was incorporated
in the
soil water flow model introduced by Hopmans
and Dane (2).
The objective was to investigate the combined
effects of
temperature and hysteresis on soil water movement,
which may
7
be especially important in predicting actual field
situations.
THEORETICAL
Temperature Dependent Hydraulic Properties
Assuming onedimensional flow, the general water flow
equation in its pressure head form can be written as (4):
ah 3h
C(h,T)=  {K(h,T)[ + 1]} , [1]
at az az
Where C is the specific water capacity (slope of water
retention curve), T is temperature, z is distance (0 at ref
erence level and > 0 above it), t is time, h is soilwater
pressure head, and K denotes the hydraulic
conductivity.
Hopmans and Dane (2) derived the following expression to
compute the soilwater pressure head (hT) at any temperature
(T), assuming knowledge of a reference soilwater pressure
head value (href) at a reference temperature
(Tref):
hT = a(T)href , [2]
8
where a(T) is defined as
a(T) = 1 + (TTref)Y(T)
and y(T) denotes the temperature coefficient of surface
ten
sion of soil water. Furthermore, the water capacity C(h) was
found to be dependent on temperature by
C(hT) = (1/a(T))C(href) , [3]
while the hydraulic conductivity, KT(e), at any temperature,
T, was calculated from
KT = (Pref/1T)Kref , [4]
where kref and PT denote the viscosity of water (Ns/m
2
) at
the reference temperature and the soil temperature in ques
tion, respectively, and Kref is the hydraulic conductivity
value at the reference temperature. The effect of tempera
ture dependent hydraulic properties on soil water movement
was investigated by Hopmans
and Dane (2), using Eq. [2]
through [4]. Hopmans and Dane (3) employed Eq. [2] through
[4] to obtain solutions to Eq. [1] using the DouglasJones
9
approximation (implicit method with implicit linearization).
Hysteresis
The following analyses pertain to various hysteresis
theories, which were initially introduced by Poulovassilis
(see references) and further developed by Mualem. A main
drying curve refers to the relationship between soil water
pressure head and volumetric water content when the soil
water pressure head in an initially "saturated" soil is
decreased until a limiting low water content is reached. A
main wetting curve is defined as the relationship between
soil water pressure head and water content when the soil
water pressure head is increased until saturation, starting
at the limiting water content value. When drying starts at
some soil water pressure head along the main wetting curve,
one refers to a primary drying curve.
I. Universal independent domain model
Mualem (6) distinguishes between two parameters that
maximum
pore diameter
pore opening
FIG. 1. Schematic diagram of a pore.
10
11
pore fills or drains is characterized by the parameters p
and r, respectively. A bivariate pore water distribution
function f(r,p) can be defined that describes the relative
volume of pores of radii p p+dp having openings of radii
r + r+dr:
f(r,p)drdp = dVp(rr+dr,pp+dp)/V , [5]
V being the total volume of the sample and dVp the change in
waterfilled pore volume. In normalizing r and p we obtain:
rRmin
pRmin
r= and p= , [6]
RmaxRmin
RmaxRmin
where R denotes specific values for p or r. By the capillary
law, all size measures, r, p, and R, are proportional to
1/h, where h is the corresponding soil water pressure head.
Rmax and Rmin correspond to hmax (at maximum water content,
Ou) and hmin (at residual water content, Omin), respec
tively, figure 3. The radii F and p change in the range
from zero to one, assuming that both r and p vary between
Rmin and Rmax. In addition, we define
0
=66min
,
where e and
0
min are the actual and residual water content. After any
number of wetting and drying cycles, O(R) can then be
obtained from integration of f(r,p) over the domain of
12
waterfilled pores. R is defined by Eq. [6],
where p or r is
replaced by R, and O(R) conveniently
replaces O(F,p).
However, so far no indication has been given as
to how
f(F,p) can be determined. Assuming that the
probability den
sity functions of F and p are independent, the
bivariate
distribution function can be written as the product
of the
two marginal distribution functions, i.e.,
f(r,p) = b(F)l(p) . [7]
By definition, f(F,p), b(F), and 1(p) are stricktly
posi
tive. Equation [7] states that the pores of any
group F have
the same distribution function 1(p). Similarly,
any p has a
pore distribution defined by b(F). The pore
water distribu
tion function, f(F,p), is mapped in figure 2. The area
of
the rectangles represents the total pore volume probability
space (0 F ( 1, 0 ( p 1). In figure 2a,
it is assumed
that when h(A) changes to h(R+dR), as during wetting, all
pores having radii R p R4+dR become waterfilled. In
a
drainage process, when h(R) diminishes to h(RdR), the
pores
of the group with radii of openings F in the range RdR F
4 1 having pore radii RdR 4 p 4 R are drained
(Figure 2b).
Since the domains of F and p are positive definite,
L(R) and B(R) can be defined as:
13
.......... ..... 0. .
/.. ...***. :.*..: .:. :
.......... /
........... 7
AdA RdR
wetting
1 1
(a)(b
1 0
O
..*..............
t.......... ...
... ........ X
A......
...
...../
.... ....
.S. .....
/ I
. ..... ..
... ....
./... ..
.. .... X
FI. . hefild.or dagas.n he.
lae shdoe
domain) for (a the main.wettng.process,
(b the main.dry
ing poces, and(C) or te priary ryin.scaning.urve
14
L(f) = f l(p)dp and
0
B(f) = f b(r)dr
0
so that L(0) = B(0) = 0. The effective water content after
wetting along the main wetting curve (Gw(R)) can be calcu
lated from (Fig. 2a):
R1
Ow(R) = f l(p)dp f b(?)d? = L(R)B(1) . [9]
0 0
Assuming B(1) = 1, and since h(R) is uniquely defined
(whether R is p or r), we can therefore calculate L(h):
L(h) = Ow(h)
, [10]
from which it can be shown that L(1) = L(hmax) = Ou
Following the main drying curve, the water content
Gd(R) is obtained from figure 2b:
2 1
Od(R
)
= f l(p)dp J b(r)dF
+
0 0
1 R
f1(p)dp fb(F)d , [11]
S0
[8]
15
which, after integrating and rearranging,
becomes
Od() 
B(R) = 
Gu Ow(R)
or as a function of h
Od(h) 
Ow(h)
B(h) =
O
u
 w(h)
[12]
[13]
Given the functions B(h) and L(h), we can
now derive
h(0)curves for any arbitrary scanning process.
For example,
the primary drying scanning curve for the drainage
depicted
in figure 2c can be calculated from:
( R
10 R
R1
= f 1(p)dp f b(?)d? +
0 0
S
(p)dp
b()d?
, [14]
0
which becomes
0 hminh= Ow(h) +
( hhi
n
h
[Ow(hl)
0
w(h)] [Od(h)0w(h)]
, [15]
[GuOw(h)]
16
Shl
where (h h indicates a wetting process from
h = hmin to hl, followed by a drying process where h
decreases from h = h
i
to h.
With this model, the 0(h)relationship in
the hystere
sis domain can be explicitly predicted from the two main
boundary curves alone. In comparing his hysteresis model
with experimental results, Mualem noted that the proposed
model deviated substantially from measured scanning curves
if hysteresis takes place at hvalues larger
than the air
entry value. This behavior is characteristic of independent
domain models.
II. A dependent domain model
Pores can only drain when, in addition to a continuous
water phase, there is a continuous air phase. Therefore,
drainage of pores in the domain of airentry values will
depend on the status of neighboring pores (9).
Pore water blockage against air entry is defined by the
function Pd(
0
) (bar indicates average domain dependence
factor Pd(G) for main wetting and drying curve):
17
Pd(=) = Ao / AG
O
, [16]
where AO is the actual decrease in water content on drying
from h
i
to h while AGOo is the value predicted for a drying
process where pore blockage does not occur (i.e., indepen
dent domain). For sufficiently low water contents, airentry
restrictions are negligible,
so that Pd(
0
) = 1. Pd(G) is
assumed to be a function of the water content only and not
of the history of the drying process (14).
Upon wetting at low water contents, one might expect
pore air blockage against water entry. It was verified from
experiments that air blockage was only of minor importance
(9) and will therefore be neglected.
In using the dependent domain model we need to deter
mine, in addition to L(h) and B(h), a third functional rela
tionship, Pd(G). This function can be calculated from the
two main curves, and a primary drying curve.
L(h) can still be calculated from Eq. (10). From the
main drying curve we can write
AG = und(h) = Pd(G)Ano = Pd(i)[Gufdo(h)] , [17]
where GuOd(h) is the volume of pores actually drained in
18
the drying process, while OuOdo(h) pertains
to the drained
water volume as calculated with the independent
domain
model. Equation [17] can be transformed,
using Eq. [11],
to
OuOd(h) = Pd(e)[lB(h)].[OuL(h)]
[18]
Since Pd(
0
) was assumed to be a function of the final water
content only, a third expression can be found
from the pri
mary drying curve (Eq. [14])
AG = w,(hl)
= Pd(G)AGo
(hmin
= Pd(G)[lB(h)][L(hl)L(h)]
. [19]
After solving for L(h) as in the independent
domain model,
Eq. [18] and [19] can be used to determine simultaneously
B(h) and Pd(
0
) with the aid of the measured main drying
curve and one measured primary drying curve. Once L(h),
B(h) and Pd(
0
) have been obtained from the main wetting,
main drying, and one primary drying curve,
any hysteretic
path can be predicted.
19
III. Extended similarity hypothesis.
Less experimental data are required if one of the dis
tribution functions b(F) or 1(p) is known. Mualem (7)
assumed
f(r,p) = b(f)b(p) . [20]
In such a case we have only one unknown function b, and the
required data for determination of f(?,p) are therefore
reduced to only one of the two main curves. The assumption,
l(p) = b(p), is valid if the areal and volumetric porosity
are equal, as in homogeneous porous media (7).
Combination of Eq. [9] and [20] results in
Ow(R
)
B(R)B(1) , [21]
which yields upon substitution of R = 1 (Ow(
1
) = Ou)
B(1) = (0u)
1 /2 ,
[22]
and therefore
B(h) = Ow(h)/(0u)
1
/
2
. [23]
20
Substitution of 1(p) = b(p) into Eq. [11] yields, for the
main drying curve:
Od(h) = [2euew(h)].
[w(h)/Ou] *
[24]
Equation [24] represents a relationship between the main
wetting and the main drying curve.
IV. Modified dependent domain theory.
This last modified analysis combines the results
of the
dependent domain model (II) with the extended similarity
hypothesis (III). Mualem's (8) expression for the domain
dependence factor is obtained by substitution of Gd(h) in
Eq. [24] for Odo(h) in Eq. [17]:
Gu(
0
uOd(h))
Pd(0) 2.
[25]
(uew(h))2
Similarly, Mualem (8) derived expressions for the primary
and secondary wetting and drying scanning curves. For exam
ple, the primary drying scanning curve can be obtained from:
21
OW~h  Pd(O)A~o
.(\hmin
hl h
(hm in
where A0
0
o is determined from integration over
the rectangle
ABCD in Fig. 2c. By substitution
of B(h), as defined by Eq.
[23], 0 can be expressed solely in
terms of Ou and the main
wetting curve:
hl h)
Ow (hl)Pd(
0
)x
F [26a]
ou
Similarly, the primary wetting and secondary
wetting and
drying curves can be derived:.
hl h)
 ()d(hl) + Pd(l) x
(
0
uOw(hl)) (OW(h)Ow(hl))
ou
hl h)
hmin h
=
?)(hmi
n
hl )
+ d(02)x
(hma
x
I [ 26b]I
20 (hmi
22
(OuOw(h2)) (O,(h)Gw(h2))
, [26c]
eu
and
hmax h2)
(hmax h2
0 = e hPd(O) x
(Ouew(h)) (O,(h2)Sw(h))
[26d]
eu
METHODS
Hysteresis Model
A FORTRAN computer program, Appendix A, based on the
modified dependent domain theory developed by
Mualem (8) was
written to simulate hysteresis. Mualem compared
predicted
scanning curves with experimental data for
three porous
media: glass beads, sand, and a sandy loam
soil. Hysteresis
simulations derived using this computer
program were com
pared with those presented by Mualem (8)
for the sandy soil
to check for correspondence with Mualem's
computer simula
tions. Mualem's main wetting and drying curves
were fitted
using van Genuchten's (16) closedform analytical
model. The
domain dependence factor, Pd(
0
), for this sandy soil was
presented, as a function of 0, in Mualem (8).
23
Water Flow Model
The general water flow equation in its pressure head
form (Eq. [1]) can be solved, provided the specific water
capacity function is known. Because there are infinitely
many hysteretic curves, the water capacity function is not
uniquely defined. However, the water capacity at any point
along a scanning curve can be computed using the soil water
pressure head derivatives of Mualem's scanning curves (i.e.,
Eq. [26a] to [26d]). Differentiating Eq. [26b] with respect
to h, for example, yields
dO dOw(h) Ow(hl) dOw(h)
= : d(
0
1).( . ) , [27]
dh w dh
0
u dh
where the derivatives on the right hand side denote the
slope of the main wetting curve at soil water pressure head
h.
The volumetric water content at any soilwater pressure
head value can be calculated from the soilwater pressure
head value at the last reversal point, provided the main
wetting and drying curves, figure 3, and Pd(
0
) are defined.
All hysteresis calculations (Eq. [25] and Eq. [26]) were
done at the reference temperature, using Eq. [2]
for
24
1.25
hmin"
1.00 
sandy soil, T=201C
0.75 
E
S  primary wetting scanning curve
0.50
Amain drying curve
0.25
main wetting curve
hmax
0 1 I 1
0.05 Omin 0.10 0.15 0.20 0.25 u 0.30
O, m
3
m
3
FIG. 3. Hysteretic curves at 20
oC for the sandy soil used
in the computer simulations.
25
conversion to the reference temperature.
As in Eq. [3], the
temperature dependent water capacity
function can be written
as
1 dO
C(hT) = . [28]
a(T) dhref
w
For the water flow simulations
it was assumed that the
K(6)function was dependent only
on temperature (i.e., it
was not subject to hysteresis). The hydraulic
properties of
the sandy soil, for which water movement
was simulated, were
experimentally determined by Haverkamp
et al. (1) and are
identical to those employed in simulation
no.2 in Hopmans
and Dane (3).
The main drying and wetting curves need
to be defined
to include hysteresis in the 6(h) relationship.
The water
retention curve determined by by Haverkamp
et al. (1) was
assumed to be the main drying curve. The
main wetting curve
was calculated from an analytical relationship
between the
main wetting and drying curve (Eq. [24]).
Using the same dependency of the soil's hydraulic prop
erties on temperature as in Hopmans
and Dane (3), infiltra
tion was simulated at soil temperatures
of 20 and 40
0
C for
both a pressure head and flux boundary
condition at the soil
surface. The initial and boundary
conditions were:
26
h(z,0,T) = 0.615 m ,
0.8 z 0 m
h(0,t,T) = 0.3 m or q(0,t,T)
= 0.1369 m h

, t > 0
h(0.8,t,T) = 0.615 m , t > 0
, [29]
where h, t, and T were previously
defined and q is the flux
density of water (q < 0 for downward
flow).
RESULTS AND DISCUSSION
Hysteresis model
Our predicted primary wetting scanning curves,
figure
4, and the scanning loops, figure 5, for the sandy
soil were
identical from those presented by Mualem (8). These
results
indicated correspondence between Mualem's
and the present
computer simulation output. Therefore, the source code
was
added to an already existing water
flow model (2). Computer
simulations of temperature and hysteresis effects
on onedi
mensional, unsaturated soil water flow
were performed using
the expanded Hopmans and Dane model.
27
64.0
56.0
0O 1
E 40.0
C
2 32.0
CO
24.0
16.0
8.0
0
0.00 0.05 0.10 0.15 0.20
0.25 0.30 035
Volumetric water content
FIG. 4. Wetting scanning curves
for the sandy soil predicted
by Mualem's modified dependent
domain theory.
28
48.0
S40.0
E
U
"0
. 32.0
I 24.0
8.0
0
0.00
Volumetric water content
FIG. 5. Primary loops for the sandy soil predicted by Mua
lem's modified dependent domain theory.
29
Water Flow Model
Isothermal infiltration into the sandy profile,
with an
initial pressure head condition of 0.615
m, followed the
primary wetting scanning curve (Fig. 3).
Changing the soil
temperature from 20 to 40
0
C shifts all 0(h)curves with
equal proportion.
Constant flux infiltrations were simulated with
soil
temperatures of either 20 OC or 40 OC. The
increase in temp
erature resulted in lower water content values in the
trans
mission zone, figure 8, (3). Inclusion of hysteresis
calcu
lations had no effect on the water content distribution
for
either of the two temperature regimes. The
water content in
the wetted transmission zone attained a value that
sustained
a hydraulic conductivity near, or equal to, the
applied flux
density at the surface. This flux density was
unaltered when
hysteresis calculations were included.
The 0(h) relation
ships were altered, however, producing
differences in the
soil water pressure head profiles.
Water content profiles resulting from a constant
pres
sure head at the top boundary are shown in figure 6. Water
storage after one hour of infiltration, ignoring
hysteresis,
had increased 42.4 mm in the 20 OC profile
(solid line, no.
1) and 52.3 mm in the 40 C profile (dashed line, no. 1).
The increase in amount of infiltrated water with the
increase in soil temperature was primarily caused by a cor
30
responding increase in hydraulic conductivity.
The amount of infiltrated water was significantly less
when hysteresis was considered (19.7 and 22.5 mm for soil
temperatures of 20 OC and 40 OC, respectively), Following
the primary wetting scanning curve, the water content corre
sponding to a surface pressure head boundary condition of_
0.30 m was only .175 at 20
0
C. If hysteresis was ignored
(figure 6, profiles numbered as 1) the corresponding water
content at a pressure head of 0.30 m was .223 (main drying
curve in figure 5). The hydraulic conductivity and, there
fore, the flux density sharply decreased when hysteresis was
considered. The water content profiles in figure 6 indicate
that the temperature effect is less pronounced when hystere
sis is considered.
A program description and listing of the water flow
model is given in Appendix B. With this model, actual field
conditions may be simulated. However, the present model
assumes a homogeneous soil profile and must be altered to
include various 'Soil layers with different hydraulic proper
ties if it is to be used for field situations.
31
9, m
3
m
3
0 0.05 0.10 0.15 0.20 0.25
pressure head
boundary condition
0.2 I
IE 2 I
I I
I I
I I
II
/I
I
initialc t on IsI
0.2 I
I
/ I
/ /
/ I
/' /
E 2 /2 1 11
.c 0.4 
3~I / /
I /
I
0.8 I
Fig. 6. Volumetric water content
(6) profiles at tempera
tures of 20 (solid lines) and
40
0
C (dashed lines) after 1
hour of infiltration with a
pressure head (h) top boundary
condition.
32
LITERATURE CITED
(1) Haverkamp, R., M. Vauclin, J.
Touma, P.J. Wierenga, and
G. Vachaud. 1977. A Comparison of Numerical
Simulation
Models for Onedimensional
Infiltration. Soil Sci. Soc.
Am. J. 41:285293.
(2) Hopmans, J.W. and J.H. Dane.
1984. Numerical Solution of
the Onedimensional Water Flow Equation with
and without
Temperature Dependent Hydraulic Properties.
Dept. of
Agronomy and Soils. Departmental Series
No. 94. Alabama
Agricultural Experiment Station. Auburn
University.
(3) . 1985.
Effect of Tempera
ture Dependent Hydraulic Properties on Soil
Water Move
ment. Soil Sci. Soc. Am. J. 49:5158.
(4) Klute, A. 1969. The Movement of Water
in Unsaturated
Soils. The Progress of Hydrology. IN: Proc.
of the First
Int. Seminar for Hydrology Prof., Urbana, Ill.,
Vol.
II, 821886.
(5) Mualem, Y. 1973. Modified Approach to Capillary
Hystere
sis Based on a Similarity Hypothesis. Water
Resour. Res.
9(5):13241331.
(6) . 1974. A Conceptual Model
of Hysteresis. Water
Resour. Res. 10(3):514520.
(7) . 1977. Extension of
the Similarity Hypothesis
Used for Modeling the Soil Water Characteristics.
Water
Resour. Res. 13(4) :773780.
(8) . 1984. A Modified Dependentdomain
Theory of
Hysteresis. Soil Sci. 137(5):283291.
(9) and G. Dagan.
1975. A Dependent Domain
Model
of Capillary Hysteresis. Water Resour. Res.
11(3):452460.
33
(10) Poulovassilis, A. 1962. Hysteresis
of Pore Water, an
Application of the Concept of Independent
Domains. Soil
Sci. 93:405412.
(11) and E.E.
Childs. 1971. The Hystere
sis of Pore Water: The Nonindependence of
Domain. Soil
Sci. 112:301312.
(12) Talsma, T. 1970. Hysteresis in Two
Sands and the Inde
pendent Domain Model. Water Resour. Res. 6(3):964970.
(13) Topp, G.C. 1971a. Soil Water
Hysteresis in Silt Loam
and Clay Loam Soils. Water Resour. Res. 7(4):914920.
(14) . 1971b. Soilwater Hysteresis:
the Domain
Model Theory Extended to Pore Interaction
Conditions.
Soil Sci. Soc. Am. Proc. 35:219225.
(15) and E.E. Miller.
1966. Hysteretic Moisture
Characteristics and Hydraulic Conductivities
for Glass
bead Media. Soil Sci. Soc. Am. Proc. 30:156162.
(16) Van Genuchten, M.Th. 1978. Calculating the Unsaturated
Conductivity with a New Closedform
Analytical Model.
Water Resour. Program. Dept. of Civil
Eng., Princeton
Univ., 78WR08.
34
APPENDIX A
Fortran Program for Simulation of
Hysteresis
Program description
Hysteresis is simulated thereby using
Eq. [25] and
[26]. After reading the input data
(series of pressure head
values) in DAINl, the parameter IPA is set
to either 0 or 1.
IPA defines whether initially the soil is following
the main
wetting (IPA = 0) or main drying curve (IPA
= 1). After
each pressure head value, the program checks if
a reversal
point occurs. In that case:
UU = [H(I)H(I2)]*[H(I)H(I)] > 0
where the numbers in parentheses indicate successive
pres
sure head values with time. For UU > 0, water content
values
are then determined following statement
number 100 in the
program listing through statement
number 500 (Appendix Table
Al). IPA is updated (soil is wetting
or drying) and N
(number of times a reversal point occurs)
is increased with
1. A(N+I) takes the value of the volumetric water content at
this last reversal point.
Two different procedures are followed to determine the
domain dependance factor Pd(e). If wetting occurs (e.g. Eq.
35
[26a]), Pd is calculated for the water content at the last
reversal point (Pd(eN) in function P). However, when drying
occurs (e.g. Eq. [26b]), Pd can only be implicitly deter
mined. For a given starting value of e (last known water
content), an iteration process starts, until the two latest
calculated water content values differ by less than 0.0001.
The main wetting and drying curves, necessary to compute the
water content at any scanning curve, are defined in the
function FTH. Finally, a plot is generated that shows the
scanning curves for the given pressure head sequence (Fig. 3
and 4).
36
Appendix Table Ai. Listing
of Hysteresis Program.
//HYSTER JOB (AYL59,124), 'JAN
HOPMANS ,NOTIFY=AYL59JHMSGCLASS=P
/*ROUTE PRINT RMT4
/*JOBPARM LINES=4K,TIMEO09
/7 EXEC FTGVCLG
//SYSIN DD *
C GOPTIONS DEVICE=CALCOMP;
CTHIS PROGRAM PLOTS THE MAIN
WETTING AND DRYING CURVE FOR
THE SANDY
C SOIL, DESCRIBED IN MUALEM(1984).
DIMENSION A(1000),TH(1000),U(1000)rUR(1000)
COMMON THS
CALL PLOTS(0,0,0)
C FIRST WE PLOT THE MAIN
DRYING AND WETTING CURVE
C
IPAR=O
J=O
DO 40 1=1,61,5
J=J+1
U(J)=1.0*I
TH(J)=FTH(IPARU(J))
WRITE(6,1) JTH(J),U(J)
40 CONTINUE
1 FORMAT(I5,2X,2E12.5)
IPAR=1
DO 55 1=1,61,5
J=J+1
U(J)=1.0*I
TH(J)=FTH( IPARU(J))
WRITE(6,1) J,TH(J),U(J)
55 CONTINUE
J=J+1
JJ=J
J8=J+7
J9=J+8
J16=J+15
J17=J+16
J24=J+23
J25=J+24
37
10
c
c 
c
c
WRITE(6,10) (U(I),1=J17,J24)
READ(2,10) (U(I),I=J25,J32)
WRITE(6,10) (U(I),I=J25,J32)
FORMAT(8F8.1)
INITIALLY THE SOIL IS FOLLOWING THE MAIN WETTING OR DRYING CURVE
IPA=O WETTING , IPA1l  DRYING
IPA =1
IPAR=IABS (IPA1)
N = 0
1=1
A(JJ) =FTH(IPAU(JJ))
TH(JJ) =A(JJ)
WRITE(6,1000) JJ,IPAIPARU(JJ),TH(JJ)
JJJ=J32
JJ1=JJ+1
DO 600 1=JJ1,JJJ
c IF(I.EQ.2.AND.U(2).LT.U(1).AND.IPA.EQ.0)
c IF(I.EQ.2.AND.U(2).GT.U(1).AND.IPA.EQ.1)
IF(I.EQ.JJ+1) GO TO 50
UU = (U I1 ( )*U I1 ()
WRITE(6,11) UU
11 FORMAT(3X'UU ',E12.5)
IF(UU.GT.0.0) GO TO 100
50 IF(N.EQ.0) TH(I)=FTH(IPAU(I))
IF(N.EQ.0) GO TO 500
GO TO 300
100 IPAR=IPA
IPA=lABS (IPA1)
N=N+1
tR(N) =U( I1)
GO TO 100
GO TO 100
IF(N.EQ.1) A(N+1)=FTH(IPARtJR(N))
IF(N.GT.1) A(N+1)=TH(I1)
WRITE(6,101) NUR(N),A(N+1)
101 FORMAT(# N UR A(N+1) 'j15j2El2.5)
300 IF(IPA.EQ.0) TH(I) = A(N+1) +P(IPArUR(N),U(I),A(N+1))
IF(IPA.EQ.0) GO TO 500
XXX=TH( Ii)
350 TH(I)= A(N+1) +P(IPAUR(N),U(I),XXX)
IF(ABS(TH(I)XXX).LE.0.0001) GO TO 500
XXX=TH( I)
WRITE(6,1000) IIPArIPARU(I) ,TH(I)
GO TO 350
500 IF(TH(I).GT.THS) TH(I)=THS
IF(IPA.EQ.1.AND.TH(I).GT.FTH(1,U(I))) TH(I)=FTH(1,U(I))
IF(IPA.EQ.0.AND.TH(I).LT.FTH(0,U(I))) TH(I)=FTH(0,U(I))
600 WRITE(6,1000) IIPAIPARU(I),TH(I)
CONTINUE
1000 FORMAT(3X,316,5X,2E12.5)
CALL PLOT(0.5,0.5,3)
CALL AXIS(0.0,0.0, 'VOLUMETRIC WATERCONTENT' ,+23,8.0,0.0,0.0,0.05)
38
CALL AXIS(0.0,0.0, 'PRESSURE HEAD' ,+13,8.0,90.0,0.0,8.0)
DO 1500 K=1,JJJ
U (K) =U (K)
TH(K)=TH(K)+0 .08
WRITE(6,1100) U(K),TH(K)
1100 FORMAT(2X,2(E12.5,2X))
1500 CONTINUE
TH(JJJ+1)0O.0
TH(JJJ+2 )0.05
U(JJJ+1)0 .0
U(JJJ+2 )8.0
CALL LINE(TH,UJJJ, 1 +1, 3)
CALL PLOT(0.0,0.0,+999)
STOP
END
FUNCTION FTH(IPAU)
C THIS SUBROUTINE WILL CALCULATE THETA FROM PRESSURE HEAD AT
EITHER
C THE MAIN WETTING (PARO0.0) OR MAIN DRYING CURVE
(PAR=1.0)
C
REAL N1,M1
COMMON THS
UU= U
TF (TPA.EQO0) GO TO 100
C NEXT 7 STATEMENTS FOR DRYING CURVE
C
C N1=6.97894
C A1=0.00809
C R1=0.0340
C S1=0.365
N1=12.5284
A1=0.02 422
R1=0.06268
S1=0.310
THS=S10 .08
M11l.  (1./Ni)
TE = (1.0/(1.0+(A1*UU)**N1))**M1
FTH=R1+(SiRi ) *Eo .08
C FTH=(S1R1)*TE
C Q= (1.0+(A1*UU)**N1)**(M11.0)
C CC= (S1R1)*M1*Q*N1*(A1**N1)*UU**(N11)
C WRITE(6,50) UUCC
50 FORMAT(' PRESSURE',E12.5,'WATER CAP',E12.5)
GO TO 300
C NEXT 7 STATEMENTS FOR WETTING CURVE
100 CON t1 TINUEtT Tr
39
S1=0. 310
C THS=SiR1
THS=Si0.08
M11l.  (1./Ni)
TE = (1.0/(1.0+(A1*tUt)**N1))**M1
FTH=R1+ (SiRi) *TE 0.08
C FTH=(S1R1)*TE
C Q =(i.0+(Ai*UU)**Ni)**(Mii.0)
C CC= (S1Ri)*Mi*Q*N1*(Ai**N1)*IJU**(N11)
C WRITE(6,50) UUCC
300 RETURN
END
FUNCTION P(IPAUU1,UtJ2VTT)
C
C  TO DETERMINE THE DOMAIN DEPENDENCE FACTOR.....
C P1=1.0
C TUO0.365
C TU=TU.0329
C P1=1.0
C P1=1.06962+0. 64649*TT22 .875297*TT**2+31 .2176*TT**3
TtJO.3i0.08
TT=TT+0 .08
P1=0.881787+4. 46257*TT45 .8645*TT**2+71 .920154*TT**3
TT=TT0 .08
IF(P1.LT.0.0) P1=0.0
IF(P1.GT.i.0) P1=1.0
IF(IPA.EQ.0) GO TO 100
C FIRST FOR DRYING
A1=(FTH(0,UUi)FTH(0,UU2) )/TU
A2 =(TUFTH(0,UU2))*A1
P = P1*A2
WRITE(6,23) A1,A2,TTP1,P
23 FORMAT(' Al A2 TH P1 P',5E12.5)
GO TO 200
C NEXT FOR WETTING
100 A1=(FTH(0,UU2)FTH(0,UUi))/TU
A2 =(TUFTH(0,UU1))*A1
P =+P1*A2
WRITE(6,23) A1,A2 ,TTP1,P
200 RETURN
END
40
APPENDIX B
Water Flow
Simulation Model
Execution of
the program
The simulation
model consists
of a main program,
6 sub
routines (INIPLO,
TEMP, CORTEM,
PLO, MASSBA, and
CONDI), and
7 functions
(FTH, FTE, PP,
FC, FCC, FK, and
UIN). Appendix
Figure 1 shows
a flow chart of
the model. The
input data are
read from the
data file DATIN.
Scanning curves
are deter
mined in FTH, while
FTE defines
the main wetting
and drying
curves. Calculation
of the water
capacity function
occurs in
FC and FCC, and
the hydraulic
conductivity function
is
defined in FK.
The function UIN
provides the initial
pres
sure head condition
versus depth.
Upon execution
of the model, a
listing is printed
of
the soil's hydraulic
properties and
the initial conditions.
INIPLO generates
a plot of the
initial pressure
head and
water content distribution,
while TEMP sets
the initial
temperature distribution.
CORTEM determines
the temperature
coefficients of
pressure head and
hydraulic conductivity
as
a function of
temperature.
If the solution
does not satisfy
the criterion of
the
mass balance equation (calculations done in MASSBA), the
time step is decreased
and the solution
process is repeated.
The simulation,
on the other hand,
proceeds if the
mass bal
41
ance criterion (EPS) is met. The time step size will
increase for subsequent calculations if the mass balance is
less than 1/10 of the imposed criterion. The simulation
stops when the maximum simulation time (TEND) is reached.
The subroutine PLO genererates a plot of the water content
and soil water pressure head distributions at the predefined
times (0). CONDI allows for transient top and bottom bound
ary conditions as well as for a variable temperature distri
bution in both time and space.
42
yes
APP. FIG. Bi. Flow chart of predictorcorrector
model.
43
Appendix Table BI. Definition of Main Program Variables
and Subroutines
HYS
HYS
HYS
HYS
HYS
HYS
HYS
F(I) temperature correction factor to water retention
curve for grid point i
HO(I) pressure head value at time TI and
grid point i
H1(I) pressure head value at time TIII and grid point i
H2(I) pressure head value at time TII and grid point i
l(I,l) value of IPA at grid point i
L(I,2) nr. of reversal points at grid point
i and at time I
1(I,3) =1, if grid point i follows one of the two main
curves
=0, if grid point i follows any scanning curve
2(I,l) pressure head at last reversal point for grid
point i
2(I,2) H2(I)
2(I,3) TH(I)
2(I,4) water content at last reversal point for grid
point i
O(10) array containing the times (s) that output is
desired
TE(I) temperature at grid point i
44
TH(I) water content at grid point i
V(I) temperature correction for hydraulic
conductivity
at grid point i
V0(I) water flux at grid point i and time TI (cm s
1
V2(I) water flux at grid point i and time TII (cm s
1
)
WAT(2) amount of water stored
in profile, cm
Z(I) depth of grid point i ( < 0),
cm
ALP specifies whether the top boundary
condition is
a pressure head or flux:
= 0, pressure head
= 1, flux
DELMO change in stored water over current
time step
(cm), calculated from 2 consecutive
water
content profiles
DELFLU change in stored water over current
time step
(cm), calculated from fluxes at boundaries
DT time step (s)
DZ space step (cm)
EMB absolute mass balance at current
time step (cm)
EPS criterion for mass balance
45
IPA specifies whether grid point
is drying or wetting
= 1, drying
= 0, wetting
NO number of times output is desired
NZ number of space steps
NZl number of grid points, NZ+1
NIZ NZl  1
N2Z NZl  2
OVERAL relative mass balance since start of
simulation (%)
REMB relative mass balance at current time step (%)
SCALF residual water content of specific soil
TEND end of simulation (s)
THS saturated water content of specific soil
TI time since start of simulation (time level j)
TII time since start of simulation (time level j+l)
TIII time since start of simulation (time level j+1/2)
TIHR TII (hrs)
UBOT bottom boundary condition (cm)
UTOP top boundary condition, pressure head (cm) or
flux (cm
s
1
)
46
ZBOT depth of profile (cm)
CONDI provides boundary conditions
and temperature
distribution with depth and
time
CORTEM calculates temperature
correction for water
retention curve and hydraulic
conductivity
function
FCFCC calculates water capacity
from pressure head
and temperature
FK computes hydraulic conductivity
from pressure
head and temperature
FTHFTE computes water content
from pressure head and
temperature
INIPLO plots initial pressure
head and temperature
distribution
MASSBA determines mass balance
PLO at specified times output is
printed and plotted
PP computes domaindependent
factor
TEMP provides initial temperature
distribution
UIN provides initial pressure
head distribution
47
Appendix Table B2.
Required Input Data
1. Input data file
Column Format
Variable
Description
15 15
NZ number of space
steps
1120 FlO.4
ZBOT profile
depth (cm)
2130 F1O.4
UTOP top boundary
condition,
pressure head (cm)
or
flux
(cm
h
1
)
3140 FI0.4
UBOT bottom
boundary condition
(cm)
4150 F10.4
DT(l) initial time
step (s)
5160 F10.2
TEND simulation
time (s)
6170 FlO.4
EPS error criterion
mass
balance
15 F5.1 ALP
see Appendix Table
1
610 15
NO number of
times that output
must be printed (TEND
included)
1115 15
IPA see Appendix
Table 1
180 8F10.1 0(8)
array containing times
(s)
that output must be
printed (TEND included)
48
2. Initial conditions
The initial conditions are listed in the function
UIN(Z),
which accepts only pressure head values, and in the subrou
tine TEMP(Z,TE) which contains initial temperature distribu
tion.
3. Soil properties
Analytical expressions for e(h), K(h), and C(h)
are defined
in the functions FTE, FK, and FCC.
4. Transient boundary conditions and temperature
distribu
tions can be defined in the subroutine CONDI.
49
Appendix Table B3. Listing of Simulation Model
//HYSMODE JOB (AYL59 ,124),r'JAN HOPMANS' ,NOTIFY=AYL59JH,MSGCLASS=P
/1* *MSGCLASS=A
/*ROUTE PRINT RMT4
//*JOBPARM LINES=20K,TIMEQO0900,F0RMS=6201
/*JOBP1\RM LINES=20K,TIMEOO09
//*OUTPUT L64
/1* EXEC LISTPARM='BEG=13'
//*SYSPRINT DD SYSOUT=(A,,L64)
//*SYSIN DD DSNAME=AYL59JH.INFTEM.LIB(VARTE) ,DISP=SHR
// EXEC ZAP
//SYSIN DD*
DSNAME=AYL59JH .OUTDAT .CNTL
//*EXEC FORTGCLG ,PARM='XREF ,MAP'
//*EXEC FTGVCLG, PARM='LIST, MAP'
//*EXEC FTGVCLG
//*EXEC FTGCCLG
//*******EXEC FTGVCLG
1EXEC FORTHCLGPARM'IXREFMAP'
//*EXEC WATFIV
//*FTO1FOO1 DD DSNAME=AYL59JH.OUTMOD.DATA,UNIT=DISK,
/*DISP=(NEW,CATLG) ,SPACE=(TRK, (5,5) ,RLSE) ,LABELRETPD=3,
/*DCB=(RECFM=FBFLRECL=80,BLKSIZE=6160)
//*FTO3FQO1 DD DSN=AYL59JH.MODEL.LIB(DATIN) ,DISP=SHRLABEL=(,, ,IN)
//*WATFIV.SYSIN DD *
//*JOB DUMMYPAGES=1000,TIME=1900
//FORT.SYSIN DD*
C
C DEBUG UNIT(9),INIT(SCALFHYS1)
C* A ONE DIMENSIONAL SIMULATION MODEL*
C* USING THE PREDICTORCORRECTOR METHOD*
C* TIME STEP : VARIABLE*
C* SPACE STEP: FIXED*
VERSION DECEMBER 1984C JAN HOPMANS
50
C ****************************************************************
C
C
C
C
C *** THE DATA ARE READ FROM DISK *****************************
INTEGER TT,HYS1(220,3)
REAL HO(220),H1(220),H2(220),TH(20),VO(220),V2(220),DT(2)
REAL Z(220),A(220),B(220),C(220),CC(220),D(220),WAT(2),0(10)
REAL TE(220),F(220),V(220),CO(220),GR(220),HYS2(220,5),SCALF
REAL SLO(220),SLOP(220)
COMMON AAAJJJNZ1,ZBOT,ALP,UTOP,EMBREMBDELMOTTDELFLU,
1 TEND,HYS1,HYS2,SCALF,THS,IFLAG
READ(3,25) NZ,ZBOT,UTOPUBOT,DT(1),TENDEPSALPNOIPA
READ(3,26) (O(I),I=1,NO)
26 FORMAT(8F10.1)
C
25 FORMAT(I5,5X,4F10.4,F10.2,F10.4,/,F5.1,215)
C
C CONVERT FLUX TOP BOUNDARY TO CM/SEC
IF(ALP.EQ.1.0) UTOP = UTOP/3600
C
C * THESE DATA ARE WRITTEN TO UNIT 6 ************************
WRITE(6,45) NZZBOTUTOPALPUBOTDT(1),TENDEPSIPA,
1 (O(I),I=1,NO)
45 FORMAT(' INITIALIZATIONS AND BOUNDARY CONDITIONS ',2X,/,
1' NR. OF SPACE STEPS .........
IS,/,
2' DEPTH OF PROFILE (CM)......',Fl0.5,/,
3' TOP BOUNDARY CONDITION ....',F10.6,' ALPHA =lF5.1,/,
4' BOTTOM BOUNDARY CONDITION .',FlO.5,/,
5' INITIAL TIME STEP (SECON)',Fl0.5,/,
6' MODEL STOPS AT ............ ',F10.2,' SECON',/,
7' ERROR CRITERION MASS BALANCE',Fl0.5,3(/),
8' WETTING OR DRYING...........'r15r/f
9' OUTPUT IS PRINTED AT ',2(/),
92X,8Fl0.1)
C
C
C C
TOP BOUNDARY CONDITION:
FLUX: ALP 1.0
PRESSURE HEAD: ALP = 0.0
BOTTOM BOUNDARY CONDITION:
PRESSURE HEAD ONLY
C
C
C REMEMBER THE DARCY CONVENTION
C
C POSITIVE FLUX  > UPWARD FLOW
C NEGATIVE FLUX  > DOWNWARD FLOW
C
C
C
C
C
C
C
C
C
C
C
C
C
C
51
C ********************************
C
C
C
C
C SOME INTIALIZATIONS
C
55 DELM 0.0
DELF =0.0
EPS =.001
IFLAG=0
NO = 1
TT =1
DZ = ZBOT/NZ
NZ1=NZ+1
N1Z=NZ1
N2Z=NZ2
T1=0.0
AAA =0.30
JJJ = 0
KKK =1
CALL PLOTS(0,0,0)
SCALFO .075
DO 56 11,fNZ1
HYS1(I,2)0O
HYS1(I,3)=1
HYS2(I,1)0O.0
HYS2(I,4)=0.0
HYS1(I,1)=IPA
56 CONTINUE
DO 57 1=1,NZ1
WRITE(6,5 8) HYS1(I,1)
57 CONTINUE
58 FORMAT(' IPAl,15)
C
C INITIAL VALUES OF ZUV AND TH AT TIME ZERO
PMAX=O.
DO 60 1=1,NZ1
Z(I) =FLOAT(I~1)*DZ
HO(I)=UIN(Z(I))
PMAX = AMIN1(PMAXHO(I))
60 CONTINUE
C
C THIS SUBROUTINE GIVES AND PLOTS THE TEMP. DISTRIBUTION IN PROFILE.
C ATm TriME ZEROI (ONLY M rb .TTtlTO BrUSED iF " TEMPERlATURE TIS CO l " ll~r 4rNSTANT WITH TIME):
52
C
C FOR TEMP. DISTRIBUTION AND/OR BOUNDARY CONDITIONSIF TRANSIENT:
IF(ALP.EQ.1.0) CALL CONDI(ZTEUTOPUBOTTIFrVH0(NZ))
WRITE(6,8)
8 FORMAT( iHi,' DEPTH, TEMPERATURE AND TEMPERATURE
CORRECTION ',
1 'FACTORS FOR'!' PRESSURE HEAD AND HYDRAULIC CONDUCITIVITY RESP')
DO 61 I=1,NZ1
WRITE(6,9) Z(I) ,TE(I) ,F(I) ,V(I)
61 CONTINUE
9 FORMAT(2X,4(2XF9.4))
C
DO 62 1=1,NZ1
TH(I)=FTE(IPAH0(I),F(I)) + SCALF
C(I)=FCC(IPAHO(I) ,F(I))
SLO(I)=C(I)
SLOP(I )=C( I)
HYS2(I,2 )=H0(I)
HYS2(I,3)=TH(I)
62 CONTINUE
C
C LISTING OF THE SOIL'S PHYSICAL PROPERTIES...
WRITE(6,10)
10 FORMAT( iHi,' THE FOLLOWING TABLE GIVES THE
HYDRAULIC PROPERTIES',
1' OF THE SOIL CONSIDERED'/7' SOIL TEMPERATURE
IS REFERENCE TEMP',
22(1),' PRESSURE WATER CONTENT '
3 'CONDUCTIVITY WATER CAPACITY' '7)
FF= 1.
Vv= 1.
DO 20 L=10,100,5
U=1 .0*1
THET = FTE(IPAUFF) + SCALF
COND = FK(UVVFFKKK)
CAP = FCC(IPAUFF)
WRITE(6,50) UTHETCONDCAP
C WRITE(1,50) UTHETCONDCAP
20 CONTINUE
DO 30 I=100,1400,100
U=1.0*I
THET = FTE(IPAUFF) + SCALF
COND = FK(UVVFFKKK)
CAP = FCC(IPAUFF)
WRITE(6,50) UTHETCONDCAP
C WRITE(1,50) UTHETCONDCAP
30 CONTINUE
DO 40 I=1500,15I500, 1 000
53
40 CONTINUE
C
C A PLOT OF INITIAL
CONDITIONS:
C
CALL INIPLO(ZHOTHPMAX)
DO 65 I=20,NZ1
CON,.= FK( .5*(HO(I)+HO(I1)),
5*CV(I)+V(I1)), 5*(F(I)+F(I1
VO(I)= CON*((HO(I)HO(Ii))/DZ)
+ CON
65CONTINUE
VO(1) =FK(HQ(1) ,V(1)
F(1) ,1)
C
C **LIST THE INITIAL VALUES
OF DEPTHTHETArPRESSURE
HEAD AND FLUX RESP.
C
+ TEMPERATURE.
WRITE(6t,6)
66 FORMAT(1H1,/,'
INITIAL CONDITIONS ARE:
',2(/),
1' NODE DEPTH THETA PRESSURE HEAD '
2' FLUX','
TEMPERATURE',' WATER
CAPACITY',2(/))
DO 70 I=1,NZ1
WRITE(6,75) I,Z(I) ,TH(I)
,H0(I) ,V0(I) ,TE(I) ,C(I)
70 CONTINUE
75 FORMAT(2XI14X,4(3XE12,5)
,6XF6.2,3XE12.5)
DELMO = 0. 0
DO 80 I=1,NZ
DELMO=DELMO(TH(I)+TH(I+1))
80 CONTINUE
WAT(1) = DZ*DELMO/2.
WRITE(6,81) TIWAT(1)
81 FORMAT(1H1,' AT
TIME ',F10.5,' WATER
IN PROFILE IS ',E12.5,'
CM')
C
C
C IF CONTSTANT FLUX AT
TOP THEN:
IF(ALP.EQ.1.0) VO(1)=UTOP
IF(ALP.EQ.1.0) GO TO
300
C
C FOR CONSTANT PRESSURE
HEAD TOP AND BOTTOM BOUNDARY
CONDITION
TII= TI + DT(TT)
C
CCCCCCCCCCCCCCCC P R
E D I C T 0 R CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
85 DO 90 I=2,NZ
A(I)=((2*DZ**2)/DT(TT))*FC(TH(I),HO(I),F(I),IE)/.
1 FK(HO(I) ,V(I) ,F(I) ,I)
B(I)=(FK(HO(I+1),V(I+1),F(I+)I+)FK(HO(11),V(I1),F(I1l),
1 I1))/(4*FK(H0(I),V(I),F(I),I))
CC(I)=2*DZ*B(I)
90 CONTINUE
C WRITE(6,91)
TII
91 FORMAT(2(/),'
PREDICTOR AT TIME ',F1O.51'
SEC '
H1(1)=UTOP
H1(NZ1 )UBOT
54
D(2)=A(2)*H0(2)  B(2)*HO(3)  CC(2)
 H1(1) + B(2)*HO(1)
D(NZ )= B(NZ)*HO(N1Z )A(NZ)*H0 (NZ ) CC(NZ )
H1(NZ1)B(NZ ) *HO~(1)
DO 100 I=3,N1Z
D(I)= B(I)*H0(I1)  A(I)*HO(I)
 B(I)*HO(I+1) CC(I)
100 CONTINUE
DO 105 I=2,NZ
C WRITE(6,104) Z(I) ,A(I)
,B(I) ,CC(I) ,D(I)
104 FORMAT(' Z A B CC Dl,2X,5E15.5)
105 CONTINUE
C
C SOLVE FOR PRESSURE HEAD BY THOMAS
ALGORITHM.
C
C(2)= 1.0/(2.0 + A(2))
D(2)= D(2)/(2.0 + A(2))
DO 120 1=3,NZ
Y = 2.0  A(I)  C(I1)
C(I) =1.0/y
D(I) =(D(I) D(I1))/Y
120 CONTINUE
C
DO 130 1=2,NZ
C WRITE(6,125) IC(I),D(I)
125 FORMAT(2X,'C D 'f13f2El5.5)
130 CONTINUE
C(NZ)0O.0
H1(NZ) = D(NZ)
N2Z = N1Z  1
DO 140 11,fN2Z
J = NZ I
H1(J)= D(J)  C(J)*H1(J+1)
140 CONTINUE
DO 160 1=1,NZ1
C WRITE(6,150) Z(I),HO(I),H1(I)
150 FORMAT(' Z UO U1',2X,3E15.5,/)
160 CONTINUE
C
CCCCCCCCCCCCCCCCCCC C 0 R R E C T 0 R CCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C WRITE(6,161) TII
161 FORMAT(2(/),' CORRECTOR AT TIME
',F1O.5,' SEC '
DO 180 1=2,NZ
D(I)=FTH(H1(I),F(I),I)
A(I)=((2.0*DZ**2)/DT(TT))*FC(D(I),H1(I),F(I),IE)/
B(I)FK(H1~(I) ,V(I+),F(I1 I1
F(1I1 V(),()
CC If =2 IvIf 0*DZ*B( a1% 1rI % T YITT r I T 7
55
D(NZ)=H0(N1Z)+(2.0A.(NZ))*H0(NZ)+2.0*B(NZ)*H1(N1Z)
1 2.0*B(NZ)*H2(NZ1)  2.0*CC(NZ)  H2(NZ1)
 H0(NZ1)
DO 200 I=3,N1Z
D(I) = 110(11) + (2.0A(I))*H0(I) HO(I+1) + 2.*B(I)*H1(I1)
1  2.0*B(I)*H1(I+1)  2.0*CC(I)
200 CONTINUE
DO 205 I=2,NZ
C WRITE(6,104) Z(I) ,A(I) ,B(I) ,CC(I) ,D(I)
205 CONTINUE
C
C SOLVE FOR PRESSURE HEAD BY THOMAS ALGORITHM.
C
C(2)= 1.0/(2.0 + A(2))
D(2)= D(2)/(2.0 + A(2))
DO 220 I=3,NZ
Y = 2.0  A(I)  C(I1)
C(I) =1.0/y
D(I) =(D(I) D(I1))/Y
220 CONTINUE
C
DO 230 1=2,NZ
C WRITE(6,125) IIC(I),D(I)
230 CONTINUE
C(NZ)0O.0
H2(NZ) D(NZ)
DO 240 1=1,N2Z
J = NZ I
H2(J) =D(J) C(J)*H2(J+1)
240 CONTINUE
DO 260 1=1,21
C WRITE(6,250) Z(I) ,HO(I) ,H1(I) ,H2(I)
250 FORMAT(# Z H TH V TEMP C',3X,6E13.4)
260 CONTINUE
C
C SETUP MASS BALANCE:
GO TO 335
C
C FOR VARIABLE FLUX TOP BOUNDARY CONDITION
C
300 TII= TI + DT(TT)
C
CCCCCCCCCCCCCCCCCCCCCCCC PREDICTOR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
301 CONTINUE
yrTI~ITr75DTTT
56
PPPP=FC(TH(I) ,HO(I) ,F(I) ,I,SLO(I))
A(I)CC(2*DZ**2)/DT(TT))*PPPP/
1 FK(HO(I) ,V(I) ,F(I) ,I)
SLOC I)=PPPP
IF(I.EQ.1) GO TO 302
B(I)=(FK(HO(I+1) ,V(I+1) ,F(I+1) ,I+1)  FK(HO(I1) ,V(I1) ,F(I1),
1 Ii) )/C4*FKCHO(I) ,V(I) ,F(I) 'I))
302 IF(I.EQ.1) B(I)= (FK(HO(I+1),V(I+1),F(I+1),I+1)
1 FK(FU1,V(I) ,F(I) ,I) )/(4*FK(IiOCI) ,V(I) ,F(I) ,I))
CC(I)=2*DZ*BCI)
303 CONTINUE
C WRITE(6,91) TIT
H1(NZ1) =UBOT2
D(1)=A(1)*H0C1)  CC(1) + BC1)*D1 D2
D(NZ)= B(NZ)*HO(N1Z )A(NZ)*H0(NZ ) CC(NZ) H1CNZ1)B(NZ)*HOCNZ1)
IF(TII.GT.95000.0) WRITE(6,308) TIIHO(1) ,HO(2)
308 FORMAT(2X,3E12.5)
DO 304 I=2,N1Z
D(I) BCI)*H0(I1)  ACI)*HO(I)  BCI)*HO(I+1) CCCI)
304 CONTINUE
DO 305 11,fNZ
C WRITE(6,104) Z(I) ,A(I) ,B(I) ,CC(I) ,D(I)
305 CONTINUE
c
C SOLVE FOR PRESSURE HEAD BY THOMAS ALGORITHM:
C(1)= 2.0/(2.0 + A~l))
D(1)= D(1)/(2.0 + A~l))
DO 306 1=2,NZ
Y = 2.0  A(I)  C(I1)
C(I) =1.0/y
D(I) =(D(I) D(I1))/Y
306 CONTINUE
C
DO 307 11,fNZ
C WRITE(6,125) IC(I),D(I)
307 CONTINUE
C(NZ)0O.0
H1(NZ) = DCNZ)
DO 310 11,fN1Z
J = NZ I
H1CJ) =DCJ)  CCJ)*H1(J+1)
310 CONTINUE
DO 311 11,INZ1
C WRITEC6,150) ZCI),HOCI),H1CI)
031 T CONf%'TINUET7
57
D2 = 2*DZ*(UTOP2 + FK(H1(1),V(1),F(1),1))/FK(H1(Il),V(1),F(1),1)
D3 = 2*DZ*(UTOP+ FK(H1(1) ,V(1) ,F(1) ,1))/FK(H1(1) ,V(1) ,F(1) ,1)
Fiji H1(2) + 02
C WRITE(6,312) D2,FU1
312 FORMAT(' D2,FU1',2E12.5)
DO 315 11,fNZ
D(I)= FTH(H1(I) ,F(I) ,I)
PPPP=FC(D(I) ,H1(I) ,F(I),I ,SLOP(I'))
A(I)=((2.0*DZ**2)/DT(TT))*PPPP/
SLOP( I)=PPPP
IF(I.EQ.1) GO TO 313
B(I)=(FK(H1(I+1) ,V(I+1) ,F(I+1) ,I+1)  FK(H1(I1) ,V(I1) ,F(I1),
1 I1))/(4*FK(H1(I),V(I),F(I),I))
313 IF(I.EQ.1) B(I) = (FK(H1(I+1),V(I+1),F(I+1),I+1)
1 FK(FU1,V(I) ,F(I) ,I) )/(4*FK(H1(I) ,V(I) ,F(I) ,I))
C WRITE(6,314) IB(I)
314 FORMAT(' I B(I)',13,E12.5)
CC(I)=2.0*DZ*B(I)
315 CONTINUE
H2 (NZ1 )UBOT
D(1=(.01~) )*HO (1) 2*HO (2 )+2.0*B(1) *D2D1D32*CC(1)
D(NZ)=H0(N1Z)+(2.0A(NZ))*H0(NZ)+2.O*B(NZ)*H1(N1Z)
1 2.0*B(NZ)*H2(NZ1)  2.0*CC(NZ)  H2(NZ1)  HO(NZ1)
DO 320 I=2,N1Z
D(I) = HO(I1) + (2.0A(I))*HQ(I)  H0(I+1) + 2.*B(I)*H1(I1)
1  2.0*B(I)*H1(I+1)  2.0*CC(I)
320 CONTINUE
DO 321 11,FNZ
C WRITE(6,104) Z(I) ,A(I) ,B(I) ,CC(I) ,D(I)
321 CONTINUE
C
C SOLVE FOR PRESSURE HEAD BY THOMAS ALGORITHM:
C(1)= 2.0/(2.0 + A(l))
D(1)= D(1)/(2.0 + A(l))
DO 325 I=2,NZ
Y = 2.0  A(I)  C(I1)
C(I) =1.0/y
325 CONTINUE
C
DO 326 1=1,NZ
C WRITE(6,125) I,C(I),D(I)
326 CONTINUE
58
331 CONTINUE
C
335 CALL MASSBA(THH2 , Z VO ,V2 DTDZVFCOGR)
C
C DO 3 I=1,NZ1
C WRITE(6,2)HYS1(I,1) ,HYS1(I,2) ,HYS1(I,3)
C2 FORMAT(3I7)
C3 CONTINUE
C DO 5 1=1,NZ1
C WRITE(6,4)HYS2(I,1) ,HYS2(I,2) ,HYS2(I,3) ,HYS2(I,4) ,HYS2(I,5)
C4 FORMAT(5E12.5)
C5 CONTINUE
C IF(TII.GT.1.O0) EPS= 0.001
WRITE(6,501) TIIDELMODELFLUfEMBREMB
501 FORMAT(' TIME DELMO DELFLU EMD REMB l,5E15.5)
TT = 1
IF(EMB.GT.EPS) GO TO 510
IF(EMB.LT.0.1*EPS) DT(TT)= 1.5*DT(TT)
IF(TII.LT.7200.AND.DT(TT) .GT.100.0) DT(TT)=100.0
C
C TIME STEP IS DECREASED IF THE REL. MASS BALANCE IS TOO LARGE:
C
C IF(TII.GT.50.0.AND.REMB.GT.0.5) GO TO 510
GO TO 520
510 DT(TT) = 0.5* DT(TT)
TII = TI + DT(TT)
IF(TII.GT.O(NO)) GO TO 900
IF(TI.EQ.0.0) GO TO 520
IF(TI.GT.1000.0.AND.DT(TT).LT.0.5) GO TO 520
C DO 515 I=1,NZ1
C TH(I)=FTH(HO(I) ,F(I) ,I)
C515 CONTINUE
IF(ALP.EQ.1.0) GO TO 301
GO TO 85
520 TI = TII
WAT(1)= WAT(1) + DELMO
DELF = DELF + DELFLU
DELM = DELM + DELMO
OVERAL =((ABS(DELMDELF))/DELF)*100
WRITE(6,502) TIIWAT(1),EMBREMBOVERALUTOPUBOT
502 FORMAT(' TI WAT EMB RE OVER TOP BOT',F8.1,6(Ell.3))
C WRITE(6,451) TII, WAT(1)
451 FORMAT(/,# WATER IN PROFILE AT TIME ',F1O.2,'SEC:',E12.5,'CM')
IF(TI.EQ.TEND) GO TO 1000
IF(TT .EQ7200.0) DmTTT)5.0
59
DO 700 1=1,NZ1
HYS2(I,3)=TH(I)
TH(I) = FTH(H2(I),F(I),I)
IF(TI.EQ.O(NO)) WRITE(6,250)
Z(I) ,H2(I) ,TH(I) ,V2(I) ,TE(I)
,SLOP(I)
HO(I)=H2 (I)
HYS2(I,2)=H0(I)
VO(I) =V2(I)
700 CONTINUE
IF(TI.EQ.O(NO)) GO TO 800.
GO TO 750
800 AAA = AAA  0.3
JJJ = JJJ + 1
WRITE(6,451) TI, WAT(1)
CALL PLO(ZTHTIOPMAXfHO)
NO = NO + 1
IF(TII.GT.O(NO)) GO TO
900
750 CONTINUE
WRITE(6,751) H0(1) ,TH(1)
,HO(40) ,TH(40)
751 FORMAT(' H AND
THETA1,4E12.5)
IF(ALP.EQ.1.0) GO
TO 301
GO TO 85
900 TT=2
DT(TT) =O(NO) TI
T11= TI + DT(TT)
DO 940 1=1,NZ1
HYS2(I,3)=TH(I)
TH(I) = FTH(HO(I),F(I),I)
VO(I) =V2(I)
HYS2(I,2)=H0(I)
940 CONTINUE
IF(ALP.EQ.1.0) GO TO
301
GO TO 85
1000 TIHR=TI/3600
WRITE(6,531) TIHR
DO 950 11,fNZ1
TH(I) = FTH(H2(I),F(I),I)
WRITE(6,250) Z(I) ,H2(I) ,TH(I)
,V2(I) ,TE(I)
950 CONTINUE
AAA= AAA  0.3
JJJ JJJ+ 1
WRITE(6,451) TI, WAT(1)
CALL PLO(ZTHTIOPMAXfH2)
NO = NO + 1
CvALL PLOT(0.0,0.0,n na999)
60
C DEBUG UNIT(9),INIT(SCALFIHYS1)
REAL UFHY2(220,5),AAURSCALF
INTEGER HY1(220r3),MC(220)
COMMON A1,J1,J2 ,A2 ,A3,A4,A5,A6,A7,J3,A8,A9,EIY1,HY2 SCALFTHS
IPA=HY1(IJ,1)
IPAR=IABS( IPA1)
N =HY1(ILJ,2)
IF(N.GT.0) UR= HY2(IJ,1)
AA=HY2(IJ,4)
MC(IJ)=HY1(IJ,3)
C WRITE(6,1) IJNIPAURAAISCALF
1 FORMAT(316,3E12.5)
IF(IPA.EQ.0.AND.(UHY2(IJ,2)).LT..O1) GO TO 100
IF(IPA.EQ.1.AND.(UHY2(IJ,2)).GT.0.001) GO TO 100
IF(HY1(IJ,2) .EQ.0) FTEH=FTE(IPAUF)
IF(HY1(IJ,2) .EQ.0) MC(IJ)=1
IF(HY1(IJr2).EQ.0) GO TO 500
GO TO 300
100 IPAR=IPA
IPA= lABS (IPA1)
HY1(IJ, 1)=IPA
N=N+1
tR=HY2 (IJ,2)
MC (I J)= 0
C WRITE(6,1) IJNIPAURAASCALF
IF(N.EQ. 1) AA=FTE( IPARURF)
IF(N.GT.1) AA=HY2(IJ,3) SCALF
HY2(IJ,4)=AA
300 IF( IPA.EQ. 0) FTH=AA+PP( IPAURU,AAF)
C WRITE(6,2) AAFTH
IF(IPA.EQ.0) GO TO 500
XXX=HY2(IJ,3)SCALF
350 FTH=AA+PP(IPAURUX.XXF)
XYZ = ABS(FTHXXX)
IF(XYZ.LE.0.0001) GO TO 500
XXX=FTH
C WRITE(6,351) XXXXYZ
351 FORMAT(' ITERl,2E12.5)
GO TO 350
500 IF(FTH.GT.THS) FTH=THS
IF(IPA.EQ.1.AND.FTH.GE.FTE(1,UF))
IF(IPA.EQ.1.AND.FTH.GE.FTE(1,UF))
IF(IPA.EQ.0.AND.FTH.LE.FTE(0,U,F))
IF(IPA.EQ.0.AND.FTH.LE.FTE(0,UrF))
IF(N.GT.0) HY2(IJ,1)=UR
HY1(IJ, 1)=IPA
HY1(IJ 2 )=N
HY1(IJ 3 )=MC(IJ)
C WRITE(612) FTHSCALF
2 FORMAT(2X,2E12.5)
FTH=FTH+SCALF
C WRITE(6,2) FTHSCALF
FTH=FTE C1,U, F)
MC (I J)= 1
FTIBFTE(0 ,U,F)
MC (I J)= 1
61
RETURN
END
FUNCTION FTE(IPAUF)
C DEBUG UNIT(9),INIT(SCALFHYS1)
C * DEFINITION OF MAIN WETTING AND DRYING CURVE....
REAL UFUUTEHY2 (220,5)
REAL M1,Ni
INTEGER HY1(220r3)
COMMON B1,J1,J2 ,B2 ,B3,B4,B5,B6,B7,J3,A8,A9,HY1,HY2 ,SCALFTHS
SCALF=0.075
UU=(l1./F)*U
C IF(HH.GE.1.0) GO TO 10
C GO TO 30
C WRITE(6,i) UUUIPA
1 FORMAT(2E12,5,14)
IF(IPA.EQ.0) GO TO 100
UU =UU
FTE =1 .611E+06* .212/( 1. 611E+06+ABS(UJ) **3 .96 )+0.075
THS=. 287SCALF
C WRITE(6,2) N1,Mi
2 FORMAT(2E12.5)
FTE=FTESCALF
C WRITE(6,2) TEFTE
GO TO 300
100 N1=5.19408
A1=0.0371
R1=0.07 44
S1=0.287
THS=Si SCALF
M11l.  (1./Ni)
TE= (1.0/(1,0+(A1*UU)**N1))**Mi
FTE=R1+(SiRi) *TESCALF
C WRITE(6,2) FTESCALF
300 RETURN
END
FUNCTION PP(IPAUU1,UU2,T2,F)
C DEBUG UNIT(9),INIT(SCALFHYS1)
C** COMPUTATION OF DOMAIN DEPENDENCE FACTOR
REAL SCALFUU1 ,UU2 ,T2 ,HY2 (220,5)
INTEGER HYi(220r3)
COMMON Bi,JiJ2 ,B2 ,B3,B4,B5,B6,B7,J3,A8,A9 ,HYiHY2 ,SCALFTHS
SCALF=0 .075
TU=0 .283SCALF
T2=T2+SCALF
C WRTIDTE(6,V3)4:1 NC'A'LFrTUT2 M
62
A1=(FTE(0,UU1,F)FTE(0,UU2,F))/TU
A2=(TUFTE(0,UU2,F) )*A1
PP=P1 *A2
GO TO 200
C NEXT FOR WETTING
100 A1=(FTE(0,UU2,F)FTE(0,UU1,F) )/TU
A2=(TUFTE( 0 UU1,F) )*Al
PP=+P1 *A2
c WRITE(6,2) A1,A2,PPTT,SCALF
2 FORMAT(5E12.5)
200 RETURN
END
FUNCTION FC(T1,UFIJSL)
C DEBUG UNIT(9),INIT(SCALFHYS1)
C * CALCULATION OF WATER CAPACITY AT ANY SCANNING CURVE FROM
C * WATER CAPACTITY AT2 MAIN CURVES (DEFINED IN FCC) ...
REAL UFUU,N1,M1,HY2(220,5) ,T1,TTTSCALFSL
INTEGER HY1(220f3)
COMMON AlJ1 J2 ,A2 ,A3 ,A4 A5,A6 ,A7,J3 ,A8,A9 ,HY1 ,HY2 ,SCALFTHS
IPA=HY1(IJ,1)
TU=0. 28 30 .075
SCALFO .075
IF(HY1(IJ,3).EQ.1) GO TO 300
IF(IPA.EQ.0) GO TO 200
100 B1=FCC(0,UF)
B2=((FTE(0,HY2(IJ,1),F))/TU)*FCC(0,UF)
B3=(2/TU)*(FTE(0,U,F) )*FCC(0,U,F)
B4=2 .0622824. 54188*T1+168. 526*T1**2380 .0273*T1**3
IF(B4.LT.0.0) B4=0.0
IF(B4.GT.1.0) B4=1.0
FC =B4*(B1 +B2+B3)
C B5=24.54188 + 337.052*T1  1140.0819*T1**2
C B6=SL
C B7= FTE(0,HY2(IJ,1),F)  FTE(0,UF)  FTE(0,UF)*
C 1(FTE(0,HY2(IJ,1),F))/TU + ((FTE(0,UF))**2)/TU
C FC = FC (B5*B6*B7)
GO TO 400
200 Bi = FCC(0,UF)
B2=( (FTE(0,HY2(IJ,1) ,F) )/TU)*FCC(0,UF)
TTT=HY2(IJ,4)+0.075
B3=2 .0622824. 54188*TTT+168 .526*TTT**2380 .0273*TTT**3
C WRITE(6,12) TTTB3
12 FORMAT(# TTT B3',2E12.5)
IF(B3.LT.0.0) B3=0.0
TIVf(3 GT.1. 0)N B3=1I.0
63
FUNCTION FCC(IPUF)
C ** WATER CAPACITY VALUES FROM PRESSURE
HEAD DATA AT 2 MAIN CURVES.
C DEBUG UNIT(9),I.NIT(SCALFHYS1)
REAL N1,MiUUU
COMMON B1,J1,J2,B2,B3,B4,B5,A6,A7,J3,A8,A9,HY1,HY2
,SCALFTHS
C
UU=U/F
IF(U.LT.130.O) FCCO0.0
IF(U.LT.130.0) GO TO 100
C GO TO 30
IF(IP.EQ.0) GO TO 50
UU = UU
FCC=1. 611E+06* .212*3. 96*ABS(UU) **2
.96
FCC=FCC/( 1. 611E+06+ABS(UU) **3 .96) **2
IF(UU.GT.1.O) FCCO0.0
GO TO 100
50 N1=5.19408
A1=0.0371
Ri=0 .0744
S1=0.287
M11l.  (1./Ni)
Q=(1.0+(A1*UU)**Ni) ** (Mii .0)
FCC=(SiR ) *Mi*Q*Ni* (Ai**N ) *UU**(Nii)
100 FCC=FCC/F
C WRITE(6,ii) FUUFCC
11 FORMAT(' UU FCC',3Ei2.5)
RETURN
END
FUNCTION FK(HVFFJJ)
C * HYDRAULIC CONDUCTIVITY VALUES FROM PRESSURE
HEAD DATA.
C DEBUG UNIT(9),INIT(SCALFHYSi)
REAL HVFHHWC
INTEGER JJ
HH= (i.0)*H
IF(HH.GT.0.0) HH=0.0
WC=FTH(HH,FJJ)
GO TO 30
C GOTO 10
FK=34.*i.i75E+06/(3600*(i.i75E+06+ABS(HH)**4.74))
C FK=4.428E02*124.6/(3600*(i24.6+ABS(HH)**i.77))
FK = V*FK
GO TO 30
10 HH=ABS(H)
F=0.58420234  0.09268778*HH + 0.0005i873*HH**2
FMTiO **F
64
C ** FUNCTION OF DEPTH.
C DEBUG UNIT(9),INIT(SCALF,HYS1)
UIN = 61.5
RETURN
END
SUBROUTINE INIPLO (Z ,HOrTH, PMAX)
C DEBUG UNIT(9),INIT(SCAL.FHYS1)
COMMON AAAJJJNZ1, ZBOT
REAL Z(220),TH(220),HO(220)
CALL PLOT(1.O,9.0,3)
K=NZ1+1
L =K + 1
Z(K) = 0'.0
Z(L) =ZBOT/8.0
HO(K)= 0.
HO(L)=PMAX/4.
TH (K)= 0.
TH (L) =0 .05
CALL AXIS(0.0,0.0,'PRESSURE HEAD THETA',+19,8.O,0.O,0.0,HO(L))
CALL AXIS(0.0,0.0, 'DEPTH CM' ,8,8.0,270.0,0.0,Z(L))
CALL LINE(HO, Z NZ1, 1 +1, 1)
CALL AXIS(0.0,O.3, ' ',+1,8.0,0.0,0.0,O.05)
CALL LINE(TH, Z NZ1, 1 ,+1 2)
CALL SYMBOL(2.0,8.2,0.1O,'THETA AND PRESSURE
HEAD',0.0,+23)
CALL PLOT(0.0,0.O,999)
RETURN
END
SUBROUTINE TEMP(CZ, TE)
C SETS AND PLOTS INITIAL TEMPERATUE
DISTRIBUTION IN PROFILE
C
C DEBUG UNIT(9),INIT(SCALFHYS1)
COMMON AAAJJJNZ1, ZBOT
REAL Z(220),TE(220)
DO 100 11,fNZ1
C TE(I) = ((+25.*Z(I))/ZBOT)
+ 40.
TE(I) =40.0
100 CONTINUE
K= NZ1 + 1
L = K+1
Z(K) = 0.0
Z(L) = ZBOT/8.0
TE(K) = 10.0
TE(L) = 5.0
CALL PLOT(1.0,9.0,3)
CALL ,AXIS.0,0.0, 'TEMPE.RATUJRE' ,+11,8.0,0.0,10.0,5.0)
65
C DEBUG UNIT(9),INIT(SCALFHYS1)
C
C DETERMINES THE TEMP. COEFFICIENT OF PRESSURE HEAD AND HYDRAULIC
C CONDUCTIVITY:
C
COMMON AAAJJJNZ1
REAL TE(220) ,Z(220) ,FACT(220.),VIS(220) ,T(220)
DO 200 1=1,NZ1
SUM=0.0
T(I )10 .0*TE(I)
IT=INT(T(I))
IF(IT.LE.200) GO TO 150
DO 100 J=210,IT
E = J/10.
SIG 75.594 0.1328*E0.000537*E**2+2.2719E06*E**3
DSIG= .1328  0.001074*E + 6.8157E06*E**2
GAM =(2.0/SIG)*DSIG*.1
SUM = SUM + GAM
100 CONTINUE
FACT(I)= 1.0 + SUM
GO TO 200
150 IF(IT.EQ.200) GO TO 195
DO 190 J=IT,200
E = J/10.
SIG = 75.594 0.1328*E.0.000537*E**2+2.2719E06*E**3
DSIG= .1328  0.001074*E + 6.8157E06*E**2
GAM = (2.0/SIG)*DSIG*.1
SUM = SUM + GAM
190 CONTINUE
195 FACT(I)= 1.0  SUM
200 CONTINUE
DO 400 11,FNZ1
IF(TE(I).LT.20.) GO TO 300
A = 1.3272*(20.TE(I)) 0.001053*(TIE(I)20.)**2
B= TE(I) + 105.
C = 10**(A/B)
VI =0.01002*C
IF(TE(I).EQ.20.0) VI=.01002
GO TO 350
300 A = 998.333+8.1855*(TE(I)20.)+0.00585*(TE(I)20.)**2
B= (1301./A)  3.30233
VI =10**B
350 VIS(I) = 0.01002/VI
400 CONTINUE
RTRNlT
66
COMMON AAArJJJNZI,ZBOTALPUTOPEMBREMBDELMOTTDELFLUI
1 TEND,HY1,HY2 ,SCALFrTHS
DO 3 I=1,NZ1
WRITE(6F2)HY1(II),HY1(I,2) ,HY1(I,3)
2 FORMAT(317)
3 CONTINUE
DO 5 11,fNZ1
WRITE(6,4)HY2(I,1) ,HY2(I,2) ,HY2(I,3) ,HY2(I,4)
4 FORMAT(4E12.5)
5 CONTINUE
K= NZ1 + 1
L= K + 1
Z(K) = 0.0
Z(L) = ZBOT/8.O
H(K)=0.O
H(L) =PMAX/4.
P = PMAX/4.
TH(K) =0.0
TH(L) = 0.05
IF(TIME.GT.O(1)) GO TO 20
CALL PLOT(10.0,9.0,3)
CALL AXIS(0,.0, 'OfVOLUMETRIC WATERCONTENT' ,+23,8.0,0.O,0.0,O.05)
CALL AXIS(0.0,0.0, 'DEPTH CM' ,8,8.0,270.0,O.0,Z(L))
CALL AXIS(0,0,8.0,'PRESSURE HEAD',+13,8.0,180.0,O.0,P)
CALL SYMBOL(1.0,0.3,0.10, 'SEC' ,0.0,+3)
CALL SYMBOL(6.0,6.0,0.10,'JAN HOPMANS',0.0,+11)
20 CALL LJINE(THrZNZ1,1,+1,JJJ)
CALL LINE(H, ZNZ1,1,+1,JJJ)
CALL SYMBOL( 0. 3,AAA,0 .10,JJJ,0. 0,i)
CALL SYMBOL(0.5,AAA,0.10, 'TIME' ,0.0,+4)
CALL NUMBER( 1.0 AAA,0 .10 TIME,0 .0,1)
CALL PLOT(0.0,0.0,+3)
RETURN
END
SUBROUTINE PLOCO(CG, Z)
C DEBUG UNIT(9),INIT(SCALFHYS1)
REAL Z(100) ,C(100) ,G(100)
COMMON AAA, JJJ ,NZ1, ZBOT, ALP ,UTOPrEMB, REMB, DELMO, TT, DELFLU, TEND
NZ=NZ1 1
K=NZ1
L =K + 1
DO 10 Jl1,NZ
Z(J)=Z(J+1)
10 CONTINUE
DO 20 11,NZM
67
Z(K)=0 .0
Z(L)=ZBOT/8.0
CALL PLOT(2.0,9.0,3)
CALL AXIS(0.O,0.O, 'DEPTH CM' ,8,8.0,270.0,0.0,Z(L))
CALL AXIS(0O.O, '01PRESSURE GRADIENT' ,17,8.0,0.0,0.0,1.00)
CALL AXIS(O.0, .43, 'CONDUCITIVITY' ,12,10.0,0.OO.0,O.0005)
CALL LINE(GZNZ,1,+1,1)
CALL LINE(CZNZ,1,+1,2)
CALL PLOT(O.OO.0,999)
RETURN
END
SUBROUTINE MASSBA(THH2rZV0,V2 ,DTDZ ,VFCOGR)
C
C CALCULATES MASS BALANCE OVER EACH TIME PERIOD.
C
C DEBUG UNIT(9),INIT(SCALFHYS1)
REAL TH(220) ,H2(220) ,Z(220) ,V0(220) ,V2(220) ,DT(2) ,V(220) ,F(220)
REAL CO(220),GR(220),P(220)
INTEGER TTJJJ
COMMON AAA, JJJ ,NZ1, ZBOT, ALP, UTOP, EMB ,REMB, DELMO, TT, DELFLU, TEND
NZ = NZ1 1
DO 10 I=1,NZ1
P(I)= FTH(H2(I),F(I),I)  TH(I)
C WRITE(6,349) ITH(I),H2(I),P(I)
349 FORMAT(' I TH H DELTH',14f3El5.5)
10 CONTINUE
DELMO=0.0
DO 20 I=1,NZ
DELMO = DELMO  (P(I) + P(I+1))
399 FORMAT(' DELTH ',E12.5)
20 CONTINUE
C WRITE(6,399) DELMO
DELMO = DELMO *DZ / 2.0
C
IF(ALP.EQ.1.0) V0(1)=UTOP
DO 50 I=2,NZ1
GR(I1)=(H2(I)H2(I1))/DZ
C01=FK(H2(I),V(I),F(I),T)
C02FK(H2(I1),V(I1),F(I1),I1)
CON=(CO1+C02 )/2.
C CON = FK( .5*(H2(I)+H2(I1)), .5*(V(I)+V(I~1) ),.5*(F(I)+iF(I~1)))
V2(I) = CON * ((H2(I)  H2(I1))/DZ) +CON
CO(I1) =CON
50 CONTINUE
Ti)I(ALEQ1.01) V( )1O
68
EMB ABS(DELMO  DELFLU)
REMB (EMB/ABS(DELFLU))*100
RETURN
END
SUBROUTINE CONDI (Z
TUTUBTIMF,VU)
C
C TEMPERATURE DISTRIBUTION
AND TRANSIENT BOUNDARY
CONDITIONS
C AS A FUNCTION OF TIME:
C
C DEBUG UNIT(9),INIT(SCALFHYS1)
REAL Z(220),T(220),UTUBTIMF(220)fV(220)rXX
COMMON AAAJJJNZlZBOTB3,B4,B5,A6,A7,J3,A8,A9,H1,H2
,SCTSIFLAG
C GO TO 600
xx=100000.
IF(TIM.GT.(7200.+XX))
GO TO 100
IF(TIM.LE.7200.) GO TO
200
IF(TIM.GT.XX) GO TO
300
100 Al =(2*3.141*(TIM7200))/86400
UT=0 .0000025+0.000002
5*SIN(Al)
C UT=0.000004
C IFLAG~l
GO TO 400
200 UT= 3./3600
C IF(IFLAG.EQ.1)
UTO0.0000025?SIN((2*3.14*(TIM4000))/86400)
GO TO 400
300 UT =(.../3600)*(TIMXX)/3600.
IF(TIM.GT.(XX+3600))
UT=2./3600
400 IF(TIM.GT.7200.)
UB =U
IF(TIM.LE.7200.) UB
61.5
WW =7.272E05
TA =.30.0
AO = 15.0
DD = 22.6
DO 500 1=1,NZ1
A2 = Z(I)/DD
A3 =(EXP(A2))*SIN(WW*TIM
+A2)
T(I) = TA + AO*A3
C T(I) =20.0
500 CONTINUE
GO TO 700
C600 UT = 13.69/3600.
600 UT=30.0
UB =61.5
C UT=5.0/360O
650 CTINUE1 (
69
//GO.FTO1FOO1 DD DSNAME=AYL59JH.0U
9
?DAT.CNTLUNIT=DISK,
//DISP=(NEW,CATLG),SPACE=(TRK, (5,5) ,RLSE)LABELRETPD=3,
/DCB=(RECFM=FB,LRECL8O,BLKSTZE=616O)
//*Go.FTO9FOO1 DD DSN=AYL59JH.DUMMY.DATA,UNIT=DISK,
/*DISP=(NE~,CATLG),SPACE=(TRK,(5,5),RLSE),LABELRETPD=3,
/*DCB=(RECFM=FBLRECL=132,BLKSIZE=132O)
//GO.FTO3FOO1 DD DSN=AYL59JH.WAFLOW.LIB(DATIN)
,DISP=SHRLABEL=(,, ,IN)
//GO.SYSIN DD*
S
4