Simulation 
of
~one-dimensional
-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 One-dimensional 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 F-p
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 B-l. Flow chart of predictor-corrector
model....................................42
Appendix Table A-l. Listing of hysteresis model.36
Appendix Table B-l. Definition of main program variables
and subroutines.43
Appendix Table B-2. Required input data of model.47
Appendix Table B-3. Listing of simulation model.49
SUMMARY
The pressure head form of the general water flow equa-
tion was numerically solved using the predictor-corrector
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 
soil-water 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 air-entry 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 air-entry 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 one-dimensional 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 soil-water
pressure head, and K denotes the hydraulic 
conductivity.
Hopmans and Dane (2) derived the following expression to
compute the soil-water pressure head (hT) at any temperature
(T), assuming knowledge of a reference soil-water pressure
head value (href) at a reference temperature 
(Tref):
hT = a(T)href , [2]
8
where a(T) is defined as
a(T) = 1 + (T-Tref)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 Douglas-Jones
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(r-r+dr,pp+dp)/V , [5]
V being the total volume of the sample and dVp the change in
water-filled pore volume. In normalizing r and p we obtain:
r-Rmin 
p-Rmin
r= and p= , [6]
Rmax-Rmin 
Rmax-Rmin
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
=6-6min
, 
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
water-filled 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 water-filled. In 
a
drainage process, when h(R) diminishes to h(R-dR), the 
pores
of the group with radii of openings F in the range R-dR F
4 1 having pore radii R-dR 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)d-r
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]
[Gu-Ow(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 h-values 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 air-entry 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, air-entry
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 = u-nd(h) = Pd(G)Ano = Pd(i)[Gu-fdo(h)] , [17]
where Gu-Od(h) is the volume of pores actually drained in
18
the drying process, while Ou-Odo(h) pertains 
to the drained
water volume as calculated with the independent 
domain
model. Equation [17] can be transformed, 
using Eq. [11],
to
Ou-Od(h) = Pd(e)[l-B(h)].[Ou-L(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)[l-B(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) = [2eu-ew(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
u-Od(h))
Pd(0) 2. 
[25]
(u-ew(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
u-Ow(hl)) (OW(h)-Ow(hl))
ou
hl h)
hmin h
= 
?)(hmi 
n
hl )
+ d(02)x
(hma 
x
I [ 26b]I
20 (hmi
22
(Ou-Ow(h2)) (O,(h)-Gw(h2))
, [26c]
eu
and
hmax h2) 
(hmax h2
0 = e hPd(O) x
(Ou-ew(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) closed-form 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 soil-water pressure
head value can be calculated from the soil-water 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-
A-main 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 one-di-
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 One-dimensional 
Infiltration. Soil Sci. Soc.
Am. J. 41:285-293.
(2) Hopmans, J.W. and J.H. Dane. 
1984. Numerical Solution of
the One-dimensional 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:51-58.
(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, 821-886.
(5) Mualem, Y. 1973. Modified Approach to Capillary 
Hystere-
sis Based on a Similarity Hypothesis. Water 
Resour. Res.
9(5):1324-1331.
(6) . 1974. A Conceptual Model 
of Hysteresis. Water
Resour. Res. 10(3):514-520.
(7) . 1977. Extension of 
the Similarity Hypothesis
Used for Modeling the Soil Water Characteristics. 
Water
Resour. Res. 13(4) :773-780.
(8) . 1984. A Modified Dependent-domain 
Theory of
Hysteresis. Soil Sci. 137(5):283-291.
(9) and G. Dagan. 
1975. A Dependent Domain 
Model
of Capillary Hysteresis. Water Resour. Res.
11(3):452-460.
33
(10) Poulovassilis, A. 1962. Hysteresis 
of Pore Water, an
Application of the Concept of Independent 
Domains. Soil
Sci. 93:405-412.
(11) and E.E. 
Childs. 1971. The Hystere-
sis of Pore Water: The Nonindependence of 
Domain. Soil
Sci. 112:301-312.
(12) Talsma, T. 1970. Hysteresis in Two 
Sands and the Inde-
pendent Domain Model. Water Resour. Res. 6(3):964-970.
(13) Topp, G.C. 1971a. Soil Water 
Hysteresis in Silt Loam
and Clay Loam Soils. Water Resour. Res. 7(4):914-920.
(14) . 1971b. Soil-water Hysteresis: 
the Domain
Model Theory Extended to Pore Interaction 
Conditions.
Soil Sci. Soc. Am. Proc. 35:219-225.
(15) and E.E. Miller. 
1966. Hysteretic Moisture
Characteristics and Hydraulic Conductivities 
for Glass-
bead Media. Soil Sci. Soc. Am. Proc. 30:156-162.
(16) Van Genuchten, M.Th. 1978. Calculating the Unsaturated
Conductivity with a New Closed-form 
Analytical Model.
Water Resour. Program. Dept. of Civil 
Eng., Princeton
Univ., 78-WR-08.
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(I-2)]*[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
A-l). 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 A-i. 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 (IPA-1)
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 (IPA-1)
N=N+1
tR(N) =U( I-1)
GO TO 100
GO TO 100
IF(N.EQ.1) A(N+1)=FTH(IPARtJR(N))
IF(N.GT.1) A(N+1)=TH(I-1)
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( I-i)
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=S1-0 .08
M11l. - (1./Ni)
TE = (1.0/(1.0+(A1*UU)**N1))**M1
FTH=R1+(Si-Ri ) *E-o .08
C FTH=(S1-R1)*TE
C Q= (1.0+(A1*UU)**N1)**(-M1-1.0)
C CC= (S1-R1)*M1*Q*N1*(A1**N1)*UU**(N1-1)
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=Si-R1
THS=Si-0.08
M11l. - (1./Ni)
TE = (1.0/(1.0+(A1*tUt)**N1))**M1
FTH=R1+ (Si-Ri) *TE -0.08
C FTH=(S1-R1)*TE
C Q =(i.0+(Ai*UU)**Ni)**(-Mi-i.0)
C CC= (S1-Ri)*Mi*Q*N1*(Ai**N1)*IJU**(N1-1)
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*TT-22 .875297*TT**2+31 .2176*TT**3
TtJO.3i-0.08
TT=TT+0 .08
P1=0.881787+4. 46257*TT-45 .8645*TT**2+71 .920154*TT**3
TT=TT-0 .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 =(TU-FTH(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 =(TU-FTH(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. B-i. Flow chart of predictor-corrector 
model.
43
Appendix Table B-I. 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 domain-dependent 
factor
TEMP provides initial temperature 
distribution
UIN provides initial pressure 
head distribution
47
Appendix Table B-2. 
Required Input Data
1. Input data file
Column Format 
Variable 
Description
1-5 15 
NZ number of space 
steps
11-20 FlO.4 
ZBOT profile 
depth (cm)
21-30 F1O.4 
UTOP top boundary 
condition,
pressure head (cm) 
or
flux 
(cm 
h
1
)
31-40 FI0.4 
UBOT bottom 
boundary condition 
(cm)
41-50 F10.4 
DT(l) initial time 
step (s)
51-60 F10.2 
TEND simulation 
time (s)
61-70 FlO.4 
EPS error criterion 
mass
balance
1-5 F5.1 ALP 
see Appendix Table 
1
6-10 15 
NO number of 
times that output
must be printed (TEND
included)
11-15 15 
IPA see Appendix 
Table 1
1-80 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 B-3. 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 PREDICTOR-CORRECTOR 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=NZ-1
N2Z=NZ-2
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 T-riME 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(I-1)), 
5*CV(I)+V(I-1)), 5*(F(I)+F(I-1
VO(I)= CON*((HO(I)-HO(I-i-))/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(I-1),F(I1l),
1 I-1))/(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(I-1) - 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(I-1)
C(I) =1.0/y
D(I) =(D(I) -D(I-1))/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(- a-1% 1-rI % T -YITT r I T 7
55
D(NZ)=-H0(N1Z)+(2.0-A.(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(1-1) + (2.0-A(I))*H0(I)- HO(I+1) + 2.*B(I)*H1(I-1)
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(I-1)
C(I) =1.0/y
D(I) =(D(I) -D(I-1))/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 SET-UP 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(I-1) ,V(I-1) ,F(I-1),
1 I-i) )/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(I-1) - 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(I-1)
C(I) =1.0/y
D(I) =(D(I) -D(I-1))/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(I-1) ,V(I-1) ,F(I-1),
1 I-1))/(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) *D2-D1-D3-2*CC(1)
D(NZ)=-H0(N1Z)+(2.0-A(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(I--1) + (2.0-A(I))*HQ(I) - H0(I+1) + 2.*B(I)*H1(I-1)
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(I-1)
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(DELM-DELF))/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 n-a999)
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( IPA-1)
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.(U-HY2(IJ,2)).LT.--.O1) GO TO 100
IF(IPA.EQ.1.AND.(U-HY2(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 (IPA-1)
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(FTH-XXX)
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=. 287-SCALF
C WRITE(6,2) N1,Mi
2 FORMAT(2E12.5)
FTE=FTE-SCALF
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+(Si-Ri) *TE-SCALF
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 .283-SCALF
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=(TU-FTE(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=(TU-FTE( 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 AT-2 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 3-0 .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 .06228-24. 54188*T1+168. 526*T1**2-380 .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 .06228-24. 54188*TTT+168 .526*TTT**2-380 .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) ** (-Mi-i .0)
FCC=(Si-R ) *Mi*Q*Ni* (Ai**N ) *UU**(Ni-i)
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.428E-02*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*E-0.000537*E**2+2.2719E-06*E**3
DSIG= -.1328 - 0.001074*E + 6.8157E-06*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.2719E-06*E**3
DSIG= -.1328 - 0.001074*E + 6.8157E-06*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(I-1)=(H2(I)-H2(I-1))/DZ
C01=-FK(H2(I),V(I),F(I),T)
C02--FK(H2(I-1),V(I-1),F(I-1),I-1)
CON=(CO1+C02 )/2.
C CON = -FK( .5*(H2(I)+H2(I-1)), .5*(V(I)+V(I-~1) ),.5*(F(I)+iF(I-~1)))
V2(I) = CON * ((H2(I) - H2(I-1))/DZ) +CON
CO(I-1) =-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*(TIM-7200))/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*(TIM-4000))/86400)
GO TO 400
300 UT =(.../3600)*(TIM-XX)/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.272E-05
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