 

 

 

 

(mu-tutu!

I..‘.

muumunn‘u

"-

m
m
m

 

 

these.

TA NWERSITY LIBRARIE

lllllllllilill liillllllill mu

31293 010510133

ill

 

 

 

 

 

 

 

 

 

 

 

 

 

This is to certify that the
dissertation entitled

A STUDY OF WATER MOVEMENT
IN A)! UNSATURATED SOIL
UNDER THE VELOCITY PERMEAMETER

presented by

Natalie J. Carroll

has been accepted towards fulﬁllment '
of the requirements for

Ph .D. degreein AE
Department of

Agricultural Engineering

tog/Mm

I / rMajor prifessor

 

Date August 3, 1993

MS U is an Afﬁrmatiw Action/Equal Opportunity Institution 0-12771

 

 

LIBRARY
Michigan State
University

 

 

 

PLACE ll RETURN BOX to romovo this ohookout from your record.
TO AVOID FINES rotum on or Moro duo duo.

  

DATE DUE DATE DUE DATE DUE

i
_MAL2_4 ZIOD I
{,1 A

:

| IE—T

| | |

-| IE]
LJE
Jl I l

MSU Io An Affirmative Action/Equal Opportunity Institution
W

   

 

 

 

 

     
 

 

 

 

 

 

 

 

  
  

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

”3-9.1

 

A STUDY OF WATER MOVEMENT
IN AN UNSATURATED son.
UNDER THE VELOCITY PERMEAMETER
By

Natalie J. Carroll

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

in Agricultural Engineering
Department of Agricultural Engineering

1993

ABSTRACT

A STUDY OF WATER MOVEMENT
IN AN UNSATURATED SOIL
UNDER THE VELOCITY PERMEAMETER

BY

Natalie J. Carroll

Careful management of world water resources has become of primary
importance to agricultural, civil and biosystems engineers. The demand for
high quality water is constantly increasing while the depletion of the world
fresh water supply continues at an alarming rate. Agricultural practices are
the major source of this depletion. Agricultural practices also cause
considerable degradation of the water supply in some areas. Surface water
pollution, ground water and sediment pollution and non-point sources of
pollution come from agriculture as well as other origins.

There is much room for improvement of agricultural water use.
Excessive irrigation is no longer acceptable and adding new irrigated fields
must be done with great care. Sustainable agriculture must be practiced in
all agricultural endevors. Careful water use, however, necessitates a good
understanding of soil and water interactions. Agricultural engineers are
particularly interested in the flow of water through soil for both drainage and
irrigation purposes as well as site selection for agriculture waste-water
holding facilities. To predict ﬂuid flow rates in porous media an accurate

measurement of the hydraulic conductivity, k, is needed. The velocity

permeameter, VP, can be used to measure the in-situ hydraulic conductivity
of soil water at the tillage depth. The VP is a portable instrument and has
been found to be very useful for quick and accurate site evaluations.

The objective of this research was to develop a computer program to
model water movement into the soil under the VP. The model will show the
change in soil water potential with time under the VP as well as the shape
and extent of both the wetted and saturated fronts. The model will

accomodate different soil parameters and various equipment conﬁgerations.

Copyright by
NATALIE J O CARROLL
1993

Dedication

this work is dedicated to

Dale,

Adrienne, Eric and John

 

ACKNOWLEDGMENTS

I would like to thank Dr. G. Merva and Dr. L. Segerlind, my major
professors, for their help, patience, insights, and guidance. I appreciate the
advice and direction of the other guidance committee members: Dr. V.
Bralts; Dr. R. Kunze; and Dr. R. Wallace. Dr. J. Steffe, Dr. A. Srivastava and
Dr. R. Ofoli, although not on my guidance committee, aided and enhanced my
professional development and I am grateful. I thank Dr. D. Welch for his
help and motivation and Dr. M. Esmay who was an early positive inﬂuence in
my engineering career.

I would like to thank Ms. C. Kisch for her encouragement and support
throughout my graduate studies.

I appreciate the ﬁnancial support of the Agricultural Engineering
department and of the Kellogg Corporation.

TABLE OF CONTENTS

Bag:

List of Tables .............................................. xi
List of Figures .............................................. xiii
1. INTRODUCTION AND BACKGROUND ..................... 1
1.1 Global Aspects and Objectives ............................ 1
1.2 Hydraulic Conductivity Measurements ................... 3
1.3 Laboratory Measurements of Hydraulic Conductivity ....... 4
1.3.1 Constant-Head Permeameter ........................ 4
1.3.2 Falling-Head Permeameter .......................... 5
1.3.3 Consolidation Test ................................. 5
1.4 Field Measurements of Hydraulic Conductivity ............ 10
1.4.1 Piezometers ....................................... 10
1.4.2 Pumping Tests .................................... 10
1.4.3 Cone Pentrometer .................................. 11
1.4.4 Auger Holes ...................................... 12
1.4.5 Inﬁltrometers ..................................... 13
1.4.6 Permeameters ..................................... 14
1.4.6.1 Guelph Permeameter ............................ 14
1.4.6.2 Velocity Permeameter ............................ 15

1.5 Modeling Saturated Flow ................................ 18

viii

Base

1.6 Modeling Unsaturated Flow ............................. 18
1.7 Parameter Estimation ................................. 19
1.8 The Velocity Permeameter ............................. 25
1.9 The Objectives ........................................ 22
2. ONE-DIMENSIONAL MODEL ............................. 26
2.1 Governing Equations ................................... 26
2.2 Element Matrices ...................................... 28
2.1.1 Hydraulic Conductivity ............ - .................. 29
2.1.2 Speciﬁc Storage, A ................................... 30
2.3 The Model ............................................ 33
2.3.1 The Physical Model .................................. 33
2.3.2 The Computer Model ................................ 35
2.4 The Iterative Process ................................... 37
2.5 One-Dimensional Analysis ............................. 38
2.6 Summary ............................................. 41
3. AXISYMMETRIC MODEL ................................. 42
3.1 Governing Equations ................................... 42
3.2 Element Matrices ...................................... 44
3.3 The Model ............................................ 45
3.3.1 The Physical Model .................................. 45
3.3.2 The Computer Model ................................ 49
3.4 Boundary Condition Determination ....................... 50
3.5 Apparent Hydraulic Conductivity ......................... 53

3.6 Model with Triangular Elements .......................... 53

Base

3.7 Element Matrices, Triangular Elements .................... 56
3.8 Axisymmetric Analysis, Rectangular Elements ......... ' ..... 58
3.9 Axixymmetric Analysis, Triangular Elements ............... 62
4. RESULTS AND CONCLUSIONS ........................... 66
4.1 Versatility of the Axisymmetric Model ..................... 66
4.1.1 Core Saturation and Total Flow ........................ 70
4.1.1.1 Initial Pressure Head ............................. 70
4.1.1.2 Hydraulic Conductivity Curve ...................... 71
4.1.1.3 Water Content Curve ............................. 72
4.1.1.4 Heterogenity .................................... 73
4.1.1.5 Core Diameter ................................... 73
4.1.1.6 Head Tube Diameter .............................. 74
4.1.2 Apparent Hydraulic Conductivity ...................... 74

4.2 Conclusions ........................................... 77
5. ADDITIONAL REMARKS ................................ 78
5.1 Modelling Technique ................................... 78
5.2 Determination of k and A. ................................ 78

5.3 Water Content and Isoheads ............................. 79

Appendices

Page
A. Axisymmetric Computer Program (VP) ....................... 82
B. Soil Water Pressure Head ................................. 99
C. One-Dimensional VP Computer Program ..................... 114
D. Hydraulic Conductivity and Water Content Curves ............ 131
E. Procedures for Triangular Element Analysis .................. 133
F. Pressure Head for Triangular Element ModeL ................. 136

List of References ........................................... 142

 

LIST OF TABLES

Table Bags
1. Input parameter conditions ............................... 67
2. Core saturation and total ﬂow ............................. 69
3. Apparent hydraulic conductivity ........................... 76
Appendices
B1. Head at each node - control parameters ..................... 100
B2. Head at each node - IC1 parameters ........................ 101
B3. Head at each node - IC2 parameters ........................ 102
B4. Head at each node - 103 parameters ........................ 103
BS. Head at each node - K1 parameters ......................... 104
B6. Head at each node - K2 parameters ......................... 105
B7. Head at each node - WCl parameters ....................... 106
B8. Head at each node - WC1 parameters ....................... 107
B9. Head at each node - Kr1 parameters ........................ 108
B10. Head at each node - Kr2 parameters ....................... 109
B11. Head at each node - CD1 parameters ...................... 110
B12. Head at each node - CD2 parameters ...................... 111
B13. Head at each node - HTDl parameters ..................... 112
B14. Head at each node - HTD2 parameters ..................... 113

Table Bags
D1. WP - volumetric WC values ............................... 131
D2. WP - k values .......................................... 132
F1. Head at 15 min., triangular model .......................... 138
F2. Head at 15 min., triangular model, more nodes ............... 138
F3. Head at 30 min., triangular model .......................... 140

F4. Head at 30 min., triangular model, more nodes ............... 140

LIST OF FIGURES

Eism Base
1. Skematic of the constant-head permeameter .................. 7
2. Skematic of the falling-head permeameter .................... 8
3. Skematic of a consolidometer ............................... 9
4. Skematic of the Guelph permeameter ........................ 16
5. Skematic of the velocity permeameter ........................ 17
6. Hydraulic conductivity, k, over time ......................... 23
7. Hydraulic conductivity, k, over time ......................... 24
8. Soil water pressure head versus hydraulic conductivity .......... 31
9. Soil water pressure head versus percent water content .......... 32
10. One-dimensional ﬁnite element model ....................... 34
11. Results, One—dimensional analysis .......................... 40
12. Cross-section of the axisymmetric element .................... 43
13. Axisymmetric ﬁnite element model ......................... 47
14. Plan view of the axisymmetric model. ....................... 48
15. Schematic to calculate element ﬂow ......................... 52
16. Axisymmetric ﬁnite element model, triangular elements ........ 55
17. Triangular element coordinates ............................ 57
18. Soil water pressure head distribution, axisymmetric model ..... 60
19. Results, axisymmetric analysis ............................ 61
20. Soil water pressure head distribution, triangular model ........ 64

xiv

Eisuxe Bass
21. Wetting of nodes 2-5, triangular elements .................... 65
Appenshses
B1. Control isoheads ........................................ 100
BZ. ICl isoheads ........................................... 101
B3. ICZ isoheads ........................................... 102
B4. 1C3 isoheads ........................................... 103
B5. K1 isoheads ............................................ 104
B6. K2 isoheads ............................................ 105
B7 . WCl isoheads .......................................... 106
B8. WC2 isoheads .......................................... 107
B9. Kr1 isoheads ........................................... 108
B10. Kr2 isoheads .......................................... 109
B11. CD1 isoheads .......................................... 110
B12. CD2 isoheads .......................................... 111
B13. HTDl isoheads ........................................ 112
B14. HTD2 isoheads ........................................ 113
F1. Axisymmetric model with triangular elements ................ 137
F2. Isoheads at 15 minutes, triangular elements ................. 139

F3. Isoheads at 30 minutes, triangular elements ................. 141

1. INTRODUCTION

1.1 Global Aspects and Objectives

Careful management of world water resources has become of primary
importance. Assuring an adequate supply of high quality water for all people
in years to come is a concern to engineers, politicians, heads of state and
concerned citizens. By the year 2000 many countries will experience
excessive scarcity of water due to the increasing demand on water resources
for agriculture, industry and domestic use (UN Environment Program, 1992).
The demand for water varies markedly from one country to another and
depends both on population and on the prevailing level and pattern of
socioeconomic development. Conspicuous disparity exists between developed
and developing countries. The average per capita domestic use of water in
the United States is more than 7 0 times that in Ghana (UN Environment
ngram, 1992). Furthermore, world wide water use is increasing rapidly. In
1950 1,360 km3 water was used, in 1990 that ﬁgure had risen to 4,130 km3 and
by the year 2000 it is expected that 5,190 km3 will be used each year (UN
Environment Program, 1992).

Agricultural practices are the major source of depletion of the water
supply. Averaged globally, 69% of water withdrawn is used for agricultural
purposes, 23% for industry and 8% for domestic purposes (UN Environment
Program, 1992). Agricultural practices also cause considerable degradation

of the water supply in some countries. Surface water pollution, groundwater

2

and sediment pollution and non-point sources of pollution come ﬁom
agriculture as well as other origins. Discharge from untreated or
inadequately treated wastewater into rivers, lakes and reservoirs is a
problem. Eutrophication of rivers and lakes caused mainly by the runoff of
fertilizers from agricultiural lands is signiﬁcant in some areas. Acidiﬁcation
of lakes is common in North America and in some European countries that
are heavy fertilizer users. Excessive irrigation has led to water logging and
salinization thereby accelerating land degradation. There are fears that the
rapid expansion of agriculture in desert areas may lead to over-exploitation of
groundwater for irrigation and irreparable depletion of this resource. For
example, the level of the Aral Sea is retreating because excessive irrigation
withdrawals have been reducing inﬂow ﬁ‘om the catchment area. The Aral
sealevel has dropped by 3 m since 1960 and, if the trend continues, will drop
another 9-13 m by the year 2000. With the reduced inﬂow from irrigation
returns the salinity of the Aral Sea has already increased threefold to 1
g/liter and by the year 2000 this is expected to rise to 3.5 g/liter (UN
Environment Program, 1992).

There is much room for improvement of agricultural water use.
Excessive irrigation is no longer acceptable and adding new irrigated ﬁelds
must be done with great care. Sustainable agriculture must be practiced in
all agricultural endevors. Careful water use, however, necessitates a good
understanding of soil and water interactions. Agricultural engineers are
particularly interested in the ﬂow of water through soil for both drainage and
irrigation purposes as well as site selection for agriculture waste-water
holding facilities. To predict ﬂuid ﬂow rates in porous media an accurate
measurement of the hydraulic conductivity, k, is needed. Civil engineering

applications are generally concerned with aquifer transmissivity, T, values (k

3

may be determined from T) and relatively large tracts of land. Irrigation and
drainage applications are on a smaller scale and practiced ordinarily on
unsaturated soils near the soil surface so many traditional methods for
measuring k are not useful to an agricultural engineering. The k value for an
aquifer is very different ﬁ'om the k value for the unsaturated soil in the
tillage depth. The velocity permeameter (VP, developed by Merva, 1979) can
be used to measure the in-situ hydraulic conductivity of soil water at the
tillage depth. The VP is a portable instrument and has been found to be very

useful for quick and accurate site evaluations.

The objective of this research was to develop a computer program to
model water movement into the soil under the VP. The model will show the
change in soil water potential with time under the VP as well as the shape
and extent of both the wetted and saturated fronts. The model will

accomodate different soil parameters and various equipment conﬁgerations.

1.2 Hydraulic Conductivity Measurements

Hydraulic conductivity, k, is an important parameter which is used in
determining soil water ﬂow rates for the design of reservoirs, ﬂood and
erosion control structures, channel improvements, drainage systems,
irrigation scheduling, drainage design, runoff rates and volumes,
groundwater recharge, emulating leaching and other agricultural or
hydrological processes. The measurement of k is fundamental to the design
of irrigation systems which comprise the largest single group of water users
in the United States, expending over 40 percent of the total annual water
usage (Skaggs, Monks and Huggins, 1972). Laboratory and ﬁeld
measurement of soil hydraulic properties are time consuming, often costly
and subject to large error. Also, ﬁeld soils exhibit large spatial variabilities
in their hydraulic properties, particularly unsaturated hydraulic conductivity

4

values (Nielsen et al., 1973 and Stockton and Warrick, 1971). In fact the
natural soil variations may be larger than differences between methods of
measurement (Taherian, Hummel & Rebuck, 1976; Bouwer and Jackson,
1974). Bosch (1991) analysed the errors associated with point observations of
matric potential in soil proﬁles and suggested an autocorrelation function of
the matric potential data and of the hydraulic parameters of the soil proﬁle to
be used in conjunction with an error function to better design ﬁeld

experiments.

A large number of ﬁeld measurements are required to determine k due
to the many variables involved (J abro, 1992). But a knowledge of the spatial
variability in the hydraulic conductivity is essential for understanding and
modeling of water and chemical movement (Rogers et al., 1991). Some
authors prefer to determine soil hydraulic properties from easily measuable
soil properties such as particle size distribution, bulk density, effective
porosity and carbon content and then to estimate thehydraulic conductivity

from these measurements (J abro, 1992).

1.3 Laboratory Measurements of Hydraulic Cconductivity

There are three principal types of laboratory measurements to
determine the saturated hydraulic conductivity. The constant-head
permeameter and the falling-head permeameter are similar with the
difference being that the water level is allowed to drop in the second
apperatus. The hydraulic conductivity of the soil may also be determined
from consolidation tests (Freeze and Cherry, 1979).

1.3.1 Constant-Head Permeameter

The constant-head permeameter (Figure 1) consists of a soil sample of

length L and cross-sectional area A between two porous plates. A tube with a

5

reservoir (and overﬂow) maintains the water supply level at a height H. The
outﬂow, Q, is measured from the apperatus and the hydraulic conductivity, k,

is determined by the equation

_.Q_L

_AH (1.1)

This method has been shown to work best for soils with k > 0.01 cm/min
(Klute, 1965).

1.3.2 Falling-Head Permeameter

The falling-head permeameter (Figure 2) also consists of a soil sample of
length L and cross-sectional area A placed between two porous plates. The
initial height of water in the supply tube is Ho and the value H 1 is the water

aL Ho
k-Atl'{H,) (1.2)

This method has been seen to work better for samples with k < 0.01

height after time t.

cm/min (Klute, 1965). Replacing the soil air with carbon dioxide and covering
the soil surface with sand (a practice that is sometimes employed to alleviate
slaking and dispersion) is not recommended (McIntre et al., 1979). The
authors in this paper contend that the k found is not the saturated k, as is
generally assumed, particularly for low stability soils. They note that the
value found, however, is still a useful measurement in assessing the

suitability of surface soils for irrigation.
1.3.3 Consolidation Test

The consolidation (Figure 3) test measures soil compressibility with a
consolidometer. A soil sample with cross-section A is placed in a loading cell

and a load L is applied which creates a stress, 0', where o = LIA. The

6

compressibility, a, is determined from the slope of the void ratio, 9, versus
effective stress 0,. The rate of consolidation is dependent on the
compressibility and the hydraulic conductivity, k, and the coefﬁcient of
consolidation, c,. Lambe (1951) describes the determination of c, and k using
the decline in sample thickness for each loading increment. The

consolidation test is seldom used by agricultural engineers.

 

 

4- Continuous
supply

 
   

flow

Volume v Cross-sectional
in time t: area A

Osvn

 

Figure 1. Skematic of the constant head permeameter

(Freeze and Cherry, 1979; after Todd, 1959)

 

 

 

  
  

4-—- Head toils irom
Ho to H, in iimct

r‘ Cross -sectionoi
area a

 

 

 

 

 

 

Cross -sectionol
area A

 

Figure 2. Skematic of the falling head permeameter
(Freeze and Cherry, 1979; after Todd, 1959)

 

 

 

L

  
        
  

 

  

Ring
Sonune of I
cross-secnonoi b
area A g

 

  
  
 

o o’a'o'h‘c 4

Bose/ix . 1 .

 

 

 

Figure 3. Skematic of a consolidometer

(Freeze and Cherry, 197 9; after Lambe, 1951)

 

10

1.4 Field Measurements of Hydraulic Conductivity

The hydraulic conductivity is very sensitive to any disturbance of the
soil (Bouma, 1981). Field methods have the advantage over laboratory
methods in that the soil is minimally disturbed. The scientist need not be
concerned that moisture content or other soil parameters might be altered
during transport to the laboratory. Disadvantages of ﬁeld testing methods
include contention with natural weather conditions and the need to transport

all equipment to remote sites.
1.4.1 Piezometers

Piezometers are widely used to determine soil hydraulic properties and
there are many variations on the main theme. Point piezometers are open
only over a short interval at their base while a screened (slotted) piezometer
is open over the entire thickness of a conﬁned aquifer. Measurements are
made by a quick introduction (slug test) or removal (bail test) of water (or soil
volume) and observation of recovery time. Tests are dependent on high
quality intake and the piezometer tubes are subject to corrosian and clogging
problems which may give highly inaccurate k values. Backwashing to clean
the piezometer tube, however, may also give highly inaccurate values.
Piezometers give in-situ values that are averaged over a relatively small

volume of aquifer.
1.4.2 PumpingT'ests

Pumping tests give in-situ values that are averaged over a large aquifer
volume. These are labor intensive tests since the engineer must drill a test
well and one, or more, observation piezometers. Next, a short-term pumping
test, at a constant rate, and the application of predictive formulas (graphical

time versus drawdown data; generally the Theis or Jacob methods) are used

11

to determine soil hydraulic properties. Pumping tests give values that are
averaged over large volumes. They are time consuming, costly and it is often
diﬁicult to evaluate aquifer geometry which is needed for determining what
curves to use in the Theis or Jacob analysis. Furthermore, proper pump
testing involves much art in both set-up and analysis (Freeze and Cherry,

1979).

1.4.3 Cone Penetrometer

The cone penetrometer with pore pressure measurements was ﬁrst
introduced in 1975 and is generally regarded in the geotechnical community
as one of the most efﬁcient tools for stratigraphic logging of soft soils
(Robertson et al., 1986). A coupled system was developed with a porous probe
ground-water sampler to perform soil logging, evaluate k and collect
ground-water samples in an effort to determine the extent and preferential
ﬂow pathway(s) of a soluble hydrocarbon plume in a Texas aquifer (Chiang,
Loos & Klopp, 1991). When steady penetration is stopped the excess pore
pressure decay with time may be used to calculate the coefﬁcient of
consolidation, “(gl in cohesive soils. Then the coefﬁcient of consolidation
may be used to calculate the hydraulic conductivity, (k in cm/sec), based on
the equation

k = c,*mv*yw (1.3)

where m, is the coefﬁcient of volume compressibility, in cm 2/kg and 7,, is the
speciﬁc weight of water in kg/cm’. The method does not work for granular
soils because of the rapid dissipation of the excess pore pressures. An
empirical correlation between k and soil relative density was developed for

granular soils by Chiang, Loos & Klopp (1991).

12

The cone penetrometer was designed for in-situ sampling of liquids and
gases of an entire aquifer. It’s use is reserved to large companies (or possibly
governmental agencies) due to the expense and the need for a testing vehicle

able to hydraulically advanced a coring device via a 25-ton reaction.

Pumping tests and the cone penetrometer are beyond the means, or
needs of the agricultural community. Aquifer analysis is not generally
pertinent to this group’s concerns. The hydraulic properties of interest are
those near the soil surface where tillage, drainage and irrigation take place.
Waste holding tanks, however, may need analysis at some depth below the
tillage layer. This discussion will focus on in-situ analysis methods that

apply near the soil surface and are economically feasible.
1.4.4 Auger Holes

The auger hole method of measuring hydraulic conductivity was ﬁrst
proposed by Diserens and later by Hooghoudt (vanBavel & Kirkham, 1949).
The method consists of an auger hole which is made in soil extending below
the water table and subsequently emptied. The rate of reﬁll in the hole is
dependent on soil permeability, hole dimensions and the height of water in
the hole (van Bavel & Kirkham, 1949). The hydraulic conductivity values are
averaged over a large volume (both directionally and between horizons),
which may, or may not be desirable, depending on the application.
Advantages include the use of the naturally occuring ﬂuid (soil water) not an
unknown or introduced ﬂuid, such as is the case in laboratory experiments.
Furthermore experiments are not unduly time-consuming as reported by
vanBavel & Kirkham, 1949. The auger hole method has the disadvantage
that a water table is needed which is preferably not too low, this often may

occur just in the spring when water tables are high.

13

A two auger hole method was proposed by Childs (Luthin, 1966). The
method utilizes pumping ﬁ'om one auger hole, at a constant rate, to another
auger hole until the elevations in both holes stabilize. This method provided
consistent results under ﬁeld conditions (Taherian, Hummel & Rebuck,
1976). The authors also tested different screens (aluminum, plastic and
steel) to prevent puddling of sloughing of the sides of the auger hole. No
logging of the screens was noticed, although the ﬂow rate through the soil
decreased slightly. The k values for the two hole method averaged 1.2 m/d as
compared to an average of 0.13 m/d for a single auger hole (Taherian,
Hummel & Rebuck, 197 6). These authors also compared the k values from
other methods of k determinations. Results for a single auger hole were 0.14
m/d, the piezometer gave 0.24 m/d and a permeability tank (1.2 cubic soil
block that was saturated and Q measured from under a constant horizontal
ﬂow) gave 0.16 m/d. The primary drawback of the two auger hole method is
the long time period, about 5 hours, to attain steady state conditions.

Rogers and Carter (1987) showed that large variations in k values may
be introduced by non-systematic use of the auger hole method in layered

soils.
1.4.5 Infiltrometers

Hydraulic conductivity may also be measured using an inﬁltrometer
which applies water at the soil surface. The inﬁltration equation was given
by Philip (1957)

I=Stm+At (1.4)
where I is the cumulative inﬁltration, S is the sorptivity, t is time and A is a

constant related to the saturated hydraulic conductivity. Both drip

inﬁltrometers as well as inﬁltrometers utilizing ponded rings have been used

14

to determine soil water ﬂow.

1.4.6 Permeameters

Recently two permeameters, the Guelph permeameter and the velocity
permeameter have been developed that offer further improvements over
previous methods of in—situ hydraulic conductivity measurements. Both
methods are simple to use and easily portable, and both produce results in a
relatively short time (15-20 minutes for the velocity permeameter and 60-90

minutes for the Guelph permeameter, Kanwar, et al., 1989).

1.4.6.1 Guelph Permeameter

The Guelph permeameter determines in-situ measurements in the
unsaturated zone of ﬁeld-saturated hydraulic conductivity, sorptivity and the
hydraulic conductivity-pressure head relationship (Reynolds and Elrick,
1985). These authors report measurements may be made fairly quickly (ten
minutes - two hours) depending on soil texture, initial soil wetness, water
depth in the equipment and cross-sectional area of the hole. The Guelph
permeameter, GP, was designed to operate in uncased wells which may be
dug with a soil auger or probe to a maximum depth of 8 meters. The GP
appears to average vertical and horizontal anisotropy in k which is
particularly useful in the design and monitoring of drainage and waste
disposal systems that rely on three-dimensional saturated-unsaturated ﬂow
(Reynolds and Elrick, 1985). The Guelph permeameter has the advantages of
speed, portability, low water use, equipment which is easily operated by one
person, the capacity to use any wetting liquid for inﬁltration (potentially

contaminates and lechate rates may be measured) and it may be used on

15

areas that do not have a high tolerance for disturbance (i.e. between crop
rows, on sod farms, golf courses, lawns, caps and liners of landﬁlls,

reservoirs, canals, etc.) (Reynolds and Elrick, 1985).

A disadvantage noted by the authors is that digging implements tend to
smear the walls of the well in moist-to-wet porous media containing a
signiﬁcant amount of clay and may give lower ﬂow values. A small spiked
Wheel was used to break up the smear layer and negate the effect. Another
disadvantage is that the analysis assumes an exponential relationship
between saturated k and unsaturated k, an assumption which will have
varying degrees of validity for different porous media. J abro (1992) stated
that the GP may be used on homogeneous horizons, structurally stable soild
and sandy or course loamy textured soils but cannot be used on ﬁne,
textured, organic, wet and heterogeneous soils and that the GP will give
negative and positive values and produces substantial variability within a
given soil. The method gives essentially a "point" measurement and

therefore usually requires replication.

1.4.6.2 Velocity Permeameter

The velocity permeameter, VP, was developed by Merva (1987) and is
similar in apperance to the Guelph permeameter with two major differences.
First a coring device is gently pushed into the soil, rather than using an
auger hole. Secondly the VP utilizes a falling-head method rather than a
constant-head. The rate of fall of the water column is monitored to determine
the time required to fall through distances 5h . From this information k is
calculated using Darcy’s equation. Results agree with conventional soil
coring methods. The device is transportable, easy to use, requires little water
to operate and provides horizontal as well as vertical values of k with equal

ease (Merva, 1987; Kanwar, Ahmed and Marley, 1989).

16

 

 

 

 

 

 

 

 

 

Air-inlet
tube
oufggervoir
Measuring scale
~-- inner reservoir
H
Tripod PM

Hell —-dl

 

 

 

 

 

 

 

~L» Perneuneter tip

 

 

 

 

Figure 4. Skematic of the Guelph permeameter (from Kanwar et al., 1989)

17

 

 

 

 

 

 

 

 

 

 

 

Figure 5. Skematic of the velocity permeameter (from Rose, 1988)

18

1.5 Modeling Saturated Flow

Early work in the area of saturated groundwater ﬂow was analytical,
with models that were homogeneous and had simple boundaries. Only very
simple problems with idealized conditions may be solved analytically.
Numerical models are necessary for more realistic models. The ﬁnite
difference method and more recently the ﬁnite element method are in general
use today. The ﬁnite element method allows more ﬂexibility than the ﬁnite
difference method and was the method chosen for this research. Kinzelbach
(1986) examined saturated ﬂow at length with both the ﬁnite difference
method (explicit and implicit methods) and ﬁnite element method. He also
discussed groundwater management, regional pollutant transport models
and parameter estimation. Zienkiewicz, Mayer and Cheung (1966) used a
ﬁnite element model to look at anisotropic seepage. Current groundwater
text books discuss saturated ﬂow in detail (for example: Freeze and Cherry,

197 9).

1.6 Modeling Unsaturated Flow

Research in unsaturated ﬂow is generally considered to date from 1931
when L.A. Richards proposed the concepts of water ﬂow in unsaturated soil.
Models are generally based on a two-term solution of Richard’s equation
(1931). The one-dimensional form is

’9 6 (0(9)???) "—629 (1.5)

a: 6x
where 9 is the soil water content, t is time, x is the vertical distance (positive
upwards), D is the diffusivity and k is the hydraulic conductivity. In 1960
Gardner summarized research in soil water relations and proposed hydraulic

conductivity curves based on soil suction. His approach was used in this

19

research.

Current research in modeling unsaturated ﬂow takes many forms.
Kunze and Nielsen (1982) used a ﬁnite difference solution to model vertical
inﬁltration with time. A later paper (Kunze and Nielsen, 1983) used the
model to show wetting proﬁles in the soil. J abro and Fritton (1990) used a
ﬁnite element software package (TWODEPEP, 1983) to solve Richard’s
equation in cylindrical coordinates and model water ﬂow ﬁom a percolation
test hole in homogeneous layered soil. Pressure head distribution and rate of
water ﬂow were used to compare data and the ﬂow rate was calculated with

Darcy’s equation
6(1)
q =-k(1iJ)A'5'r' (1-6)

where q is the ﬂow rate, A is the cross-sectional area and r is the radial
distance. In the unsaturated model k is a function of the matric head, it.

Neuman and Witherspoon (1971) evaluated nonsteady ﬂow with ﬁ‘ee
surfaces by redeﬁning the ﬁnite element mesh. Desai (1976) was able to
avoid some of the problems created by a variable mesh by modeling ﬂow with
an invarient mesh using the residual flow procedure. Neuman (197 3) used
the ﬁnite element method to solve the equations of transient seepage in
saturated-unsaturated porous media.

Phillip (1988) gave an overview of analytic and quasianalytic approaches
to unsaturated ﬂow. Raats (1988) commented on the same topic with further

examples.

1.7 Parameter Estimation

Water ﬂow models generally use Richard’s equation (1.5) to deﬁne the
change in moisture content with time at different locations and use

relationships such as Darcy’s law (1.6) to deﬁne ﬂow velocity. Inherent in

20

modeling these relationships is the need to quantify the hydraulic
conductivity, k. Gardner (1958) proposed an exponential equation

km.) we“) (17)
where a and alpha are constants and 1p is the suction head. Gardner also
gave the empirical power equation

a

k=
(w"+b)

where a, b and n are constants for steady-state, one-dimensional ﬂow. w" is 0

 

(1.8)

at saturation so

k =- (1.9)

Gardner (1960) gives n as between 2 and 4 for ﬁne textured soils, discusses
equation (1.8) in more detail and gives graphs of the hydraulic conductivity
as a function of water potential.

Many authors have used the exponential equation (1.3) or the power
equation (1.7) to quantify k. Raats (1983) described the theoretical
background for both forms. Unlu et.al. (1990) noted that experimental
results show the exponential relationship holds only over a limited range of
water potentials and these authors suggest using statistically estimated
values for k. Toledo, et.al. (1990) used theories of fractal geometry and
thin-ﬁlm physics to provide a basis for the power law and to deﬁne the
exponents. The relationship was written as

lease-“‘3‘“ (1.10)

where D is the Hausdorff dimension (fractal dimension, varying from 0 to 3)
of the surface between the pore space and the grains and m is the exponent in
the relation of disjoining pressure 11 and ﬁlm thickness h mark"). These
authors found values of m < 1 and 2.1 < D < 2.7 for length scales of 5 pm - 20

21

pm and D = 2 for smooth pore walls based on data for compacted sand. The
authors note that these values are empirical, that accurate measurements of
1p and k at low moisture contents is challenging and reliable data is rare.
Neuman (197 3) made the assumption that the hydraulic conductivity
may be expressed as a product of a symmetric positive-deﬁnite tensor, Km,
and a scalar function of the degree of saturation, K,
k “KmKr (1-11)

where K”, represents the conductivity at saturation and K, is based on the
volumetric water content.

J abro and Fritton (1990) used the exponential equation (1.7) in their
work. These researchers had good results using the geometric mean of k.
They noted that the wetted front continued to advance in a predominantly
horizontal direction and suggested this was due to a less permeable layer
below two more permeable layers in their model. They also recorded a thin
(approximately 10 mm) saturated zone near the water source.

Kunze and Nielsen (1983) found that arithmetic and geometric mean
values for these parameters gave unreliable results but they had success
using an integrated mean value.

The models of Burdine (1953) and Maulem (1976) which describe the
unsaturated k values are commonly used and are discussed by Schuh and
Cline (1990). These are often called the microscopic models as they are based
on microscopic pore-radius distribution. These models were not considered
for this research as they require soil parameters which are not generally

available in ﬁeld situations.
Bosch reports taht most of the methods for measuring K(<l>) using

discrete measurements are inadequate for estimating effective K(<I>) for

22

heterogeneous proﬁles unless the measurements are used to interpret the
mean behavior of the system. An analytical expression for predicting the
error which can be expected when using point observations of the matric
potential to determine the mean matric potential in a heterogeneous soil

proﬁle was derived (Bosch, 1991).

1.8 The Velocity Permeameter

The velocity permeameter, VP, consists of a sharpened coring device
which is pushed (driven) into the soil and connected to a small diameter head
tube. The head tube is ﬁlled with distilled water which percolates into the
soil core. The small diameter of the head tube acts to magnify the entry
velocity of the water into the soil encased in the core which allows more
accurate measurements. The rate of fall of the water column is monitored to
determine the time required to fall through distances 6h. A Hewlett Packard
"C" series calculator equipped with a timing module is used to accurately
time the process from which k is calculated. The head tube is reﬁlled after

the water level drops approximately 1000 mm.

Figures 6 and 7 show the plot of hydraulic conductivity over time
measured with the VP from Merva (1987 ) and Rose (1988). The hydraulic
conductivity value is initially relatively high but decreases with time and

appears to approach a constant value.

23

 

 

R (mm/hr)

 

 

time (min)

 

Figure 6. Hydraulic conductivity, k, over time (from Merva, 1987)

 

24

 

 

250_
A
E 800..
\ t
E
v 150 __ ,
l

100 _ , ‘ ‘

‘ L
50
l l l l l l

 

time (min)

 

 

Figure 7. Hydraulic conductivity, k, over time (from Rose, 1988)

 

25

1.9 The Objectives

The goal of of this research was to improve the understanding of ground
water ﬂow into the soil under the velocity permeameter, VP, by developing a
computer program to model the phenomena. Speciﬁcally, the model will
show the change in soil water potential with time under the VP as well as the
shape and extent of both the wetted and saturated fronts. The model will

accomodate different soil parameters and various equipment conﬁgerations.

2. A ONE-DIMENSIONAL MODEL

The ﬂow through the soil sample enclosed by the sample cup of the
velocity permeameter was initially simulated using a one-dimensional model.
The problem was basically one of introducing a column of water to a dry soil

and monitoring the movement of water through the soil.

2.1 Governing Equations

The three dimensional governing equation for ﬂow in an unsaturated
soil is given by the Richards equation (Richards, 1931; Freeze & Cherry,
1979)

6 611’

a 6‘11 6 6‘11 6‘1!
a—x[K(\p)a—x + $[KO‘I’05 4' 52—[K(‘y)(32_+1)] - 16—! (2'1)

where ‘1' is the soil water pressure head (pressure head), K is the hydraulic
conductivity and x, y and z are the coordinate directions, 2 vertical and
positive upwards. The variable t is time while A = tie/6‘11 is the speciﬁc
storage and where 6 is the volumetric water content. The one-dimensional

equation for ﬂow vertically into the soil is given by
a aw aw
_ _ .. — .2
az[K(W)(az+1)] "a: (2 )
for the unsaturated zone. Flow is directed downwards in the negative 2
direction. The saturated equation is

624!

K
622

- o (2.3)

 

26

27

The general solution to the equation for unsaturated ﬂow (2.2) in ﬁnite

element notation is a system of ﬁrst order differential equations
[C ] {‘1’} + [S ] {‘1’} - {F} - {0} (24)

where [C] is the capacitance matrix and {‘11} contains the nodal values, with
{‘1'}T - {‘15, ‘11,, ...,‘I’,}. The matrix [S] is called the stiffness matrix (from
structural analysis) and is usually denoted [K], with the element stiffness
matrix denoted with the lower case [k]. The variable name [S] was used to
differentiate the stiffness matrix from the hydraulic conductivity, universally

denoted by k in soil physics literature.
The global force vector is {F} and

aw, am, am,

{‘1’}= a: at 7:— (2.5)

is the time derivative of the nodal values.

There are no point source, no point sinks nor any derivative boundary
conditions in the velocity permeameter problem. The force matrix, {F},
appeared in the models because the known values of {‘11} on the upper
boundary produce values that are held in {F}.

There is one boundary condition in the one-dimensional model - the
height of water in the head tube. It was considered a ﬁxed value during each
time step in the ﬁnite element solution. The drop in head was calculated
aﬁzer each iteration with the new value used as the ﬁxed value during the
next iteration.

The ﬁnite element solution to the saturated equation, (2.3), is

[Sim-{F} (2.6)
since the time derivative vector {(1)} does not change with time in the

saturated solution.

28

The unsaturated solution looks quite different from the saturated
solution because of the addition of the capacitance matrix, [C], which
contains the time dependent parameters. Equation 2.4 gives the unsaturated
terms and is solved by separating the unknown values from the known
values and using the central difference solution technique (Segerlind, 1984)
to yield

([C]+%[SJ) {W}. - ([6] $551) {W}. +92-‘-({F}. + {F}.) (2.7)

where {‘11,} contains the known nodal values at time t, {\P,} contains the

unknown nodal values at time t+1. This equation is generally written as
[A]{‘1’}b = [P] {‘1’}. (23)
where [A] is decomposed into upper triangular form using Gaussian

elimination and the system of equations is solved using backward

substitution.

2.2 Element Matrices

The element stiffness matrix

k.- 1 - 1
= — .9
[s] 1. [_1 1] (2 )
was used where k,- is the hydraulic conductivity at node 1. Node i is the node
closest to z = 0 for each element and L is the element length. The lumped

form of the element capacitance matrix was used. This matrix was

AIL 1 O]

[c] - 2 0 1 (2.10)

where A, is calculated at node i, as described above. A, is the time dependent

multiplier (often called speciﬁc storage) and is calculated by determining the

29

slope of the ‘P — 9 curve at a known soil water pressure head.

The element force matrix

m=-Q;A-{1 1 1 1}’ (2.11)

2.2.1 Hydraulic Conductivity

The hydraulic conductivity at each node was calculated from the known
soil water pressure head, ‘1’, to form the stiffness matrix. Various methods
were tested to determine the value of the hydraulic conductivity.
Exponential equations which are often used to deﬁne the ‘l’ - k relationship
were found to be inadequate for the range of soil water pressure head
modeled. Extensive efforts were made to utilize an equation which would
model the work of Gardner for soils ranging in moisture content from
saturation to a water pressure head of -100 meters with k varying ﬁve orders
of magnitude. The equations were only able to model a portion of the k - ‘1'
curve so curve ﬁtting was not a viable solution due to the large variability of
the hydraulic conductivity (ﬁve orders of magnitude). Breaking the curve
into a series of equations reduced the order of magnitude problem but did not
give satisfactory results. The most consistent results were found by
assuming k had the same relation to water pressure head as is shown in
Gardner’s curve (Figure 8), using the curve for the Pachappa sandy loam soil.
For this approach, nineteen data points from Gardner’s graph were entered
into a data array, and a full logarithmic interpolation was used to calculate
the hydraulic conductivity at each node ﬁ'om these data. Since the saturated
k value from the graph (2* 10'3 mm/sec) was signiﬁcantly lower than values
currently being used in ﬁeld work in Michigan, a multiplication factor was
used (2.8 to give a saturated value of 0.56 mm/sec or about 20 mm/hr). The
assumption was made that the shape of the graph would not change.

30

2.2.2 Specific Storage, A

The speciﬁc storage, A, the time dependent multiplier is the ratio of the

change in water content to the change in presure head (A = ae/aW). The value
of A. at a node was calculated from known data by determining the slope of
the ‘1' - 9 (where 9 is the volumetric water content) curve between the known
node and the next data point. The graph of the ﬁeld data, similar to that
shown in Figure 9, was approximated using 14 data points and semi-log
interpolation. The actual data used was from the US. Department of
Agricultural Soil Conservation Service (1986) for a Ziegenfuss soil for a depth
of 0 - 230 millemeters. The values may be found in Appendix D.

31

 

   

-3 ‘ Indm3loon
10 _ (V

A -4 A
U 10 L
m
m
\ -s
E 10 _. Pachappa sondyloom
v tiny
x

'6
>1 10 _
f
2

-7
’J 10 I
3
e
S -a
u 10 _
Q
g -9
k 10 _.
c
>
I -n

10

 

 

Tl: l-z l l T l
-1xK)-1xH)-Ol ~10 -100 -1000

Son water pressure head (cm)

 

 

 

Figure 8. Soil water pressure head versus hydraulic conductivity

32

 

 

 

5L
4
P 31L
C
m
*c’
0 2L
u
t
9 1L
§
l l I 1
00101 1.0 100 100
SOil water pressure head (tn), ZlgenFuss soil

 

 

 

Figure 9. Soil water pressure head versus percent water content

33

2.3 The Model
2.3.1 The Physical Model

The one-dimensional ﬁnite element model is shown in Figure 10. The
nodes are numbered downward (the negative 2 direction). The element
numbers are in parenthesis. The ground surface node, node 1, is set at the
ﬁxed value of the height of water in the head tube (1.8 meters). The other
nodes are numbered incrementally in the -z direction. The elements are
numbered in the same manner. An element length of 50 millimeters was
chosen so that the time step would not be too small. The maximum time
step, At, cannot be exceeded to avoid oscillations with time dependent
problems is given by Segerlind (1984)

A A
At < m (2°12)

where A is the area of the element, 9 = 0.5 (as required for the central
difference method of solution) and k and it were calculated using data shown
in Figures 8 and 9. Analysis of At showed that the saturated values for k and
A gave the maximum At that could be used to avoid oscillations. The ﬁrst

node was saturated and drove the analysis.

34

 

 

 

node
o___ 10W
(1)
éZ—element
-50 __ 2’
A (a)
E
5 —100__ 3.
(3)
3 4
+3 —150__ J/
/
C
o -300__ 71
O
U (7)
-350__ 8 .
(8)
—400__._ 9 .
(9)
-450__ 10 I
\/

 

 

 

 

Figure 10. One dimensional finite element model

35

2.3.2 The Computer Model

The computer model consisted of three major portions: the saturated
analysis; the unsaturated analysis and the calculations to determine the drop
in head. The input values were the coordinates, the initial pressure head and
the element conﬁguration. The time dependent loop was begun after all data
had been read from a ﬁle. The time step and number of iterations were

speciﬁed with the input data.

The program initially performed some one-time calculations, such as
elemental volume and initial moisture at each node, prior to beginning the
time-step loop. The ﬁrst step in the iterative process was to check for entries
in the saturated node array, which consisted of the nodes initially
unsaturated but which became saturated as the wetting front moved past
them. Pressure head values dictated when the soil was saturated. The
saturated nodes were calculated ﬁrst because, although the pressure head
changed as the water column in the head tube dropped, they are considered
ﬁxed nodes during the unsaturated analysis.

The second major part of the computer model was the time-dependent,
unsaturated analysis which calculated the pressure head for the unsaturated
nodes using equations (2.7) and (2.8). Lastly, volumetric water ﬂow
calculations were made to determine the drop in the water surface in the
head tube. This value was subtracted from the ﬁxed nodes to give the new
boundary value for the next time step.

The global matrices were recalculated during each time step in the
unsaturated analysis because as the pressure head changes the parameters k
and 1., found in the global matrices, change. The global stiffness matrix, [S],

and capacitance matrix, [C], were created from the element matrices by the

36

direct stiffness method as described by Segerlind (1984). Note that in the
saturated solution only the stiffness matrix is needed since i. (as/as!) is zero

for saturated nodes (66 = 0).
The parameters k and i. vary with the water content and could not be

deﬁned for an element since the soil water pressure head was known only at
the nodes. Using relationships for ‘IJ-k and ‘P-G, k and A could be calculated
at the nodes. Attempts were made to deﬁne a representative water content
value for each element but the extreme variation in pressure head, ‘1’,
between adiacent nodes in the area of the wetted front caused complications
in estimating k and l. for an element. For example, initially the ﬁrst element
has a saturated soil water pressure head of 1.8 meters at the upper node and
an unsaturated pressure head of -20 meters (the initial conditions used) at
the lower node. It was impossible to accurately deﬁne a value of k or A for an
element with such a range in values (k is on the order of 4'12 m/sec for a
pressure head of -20 m. and about 1.1‘6 m/sec for the saturated node. Using
an average ‘1' value as representative for the element was unsuccessful since
it slowed the analysis to the point that the head tube was not reﬁlling within
ﬁve minutes - a criteria based on ﬁeld experiments. Noting that although the
disparity will occur at the wetted front as long as it is moving through the
soil, the interaction is probably brief and not critical to the overall problem.
To avoid the problems caused by averaging, the value of k used for each
element was calculated based on the soil water pressure head at the upper
node (nearest the ground surface) for the element. This method gave results
that were physically realistic, but may be a cause of the lack of complete

agreement between calculated and assumed k values.

37

2.4 The Iterative Process

The ﬁrst step in the computer model for each iteration was to calculate
‘1' at the saturated nodes. The saturated analysis was run on all the
saturated nodes, the ﬁxed node (node 1) and the ﬁrst unsaturated node (in
the z-direction). The hydraulic conductivity was calculated for each node (k =
k“, for the saturated nodes). The soil water pressure head, ‘11, was obtained
for each of the saturated nodes by determining the height above the
saturated front (where ‘1' was considered to be zero) and assuming a linear
relationship between the saturated front and the ﬁxed node above it.
Calculations proceed until the new pressure heads have been determined for

all the saturated nodes, except the ﬁxed node.
Next, the parameters k and A. were calculated and the global arrays [S],

[C] and {F} written. The global matrices were then modiﬁed for the ﬁxed
node (node 1). The central difference method was used to write [A] and [P]
which were then decomposed into upper triangular form. The pressure head

at each node was found with backwards substitution.

The subsequent step in the model determined the soil water pressure
head values at the unsaturated nodes. The water content at each node was
found using ﬁeld data (Figure 9) and semi-log interpolation. The change in
water content was determined by subtracting the water content at the
previous time step from the current water content at each node. This value
was multiplied by the element volume to give the elemental change in water
content. The summation of the elemental volumes yielded the total ﬂow for
the time step. The total ﬂow for the iteration divided by the area of the head
tube gave the drop in head for the iteration. The drop in head was subtracted
from the head at the ﬁxed node to give the height of the water column for the

38

next time step.

2.5 One-Dimensional Analysis

The one dimensional analysis was run for simulation time of ﬁve
minutes with a time step of 0.25 seconds (1200 iterations). Figure 11 shows
results from the one—dimensional VP model: the drop in head of node 1 and
the wetting of nodes 2 - 4.

The soil water pressure head (head, shown by node 1) was initially set to
1.8 meters in the sample cup. Node 1 remains saturated with a drop in head
and increases when the head tube is reﬁlled. After each iteration the drop in
head is calculated and the ﬁxed nodes are reduced by that value. When the
head falls below 1.0 meter, node 1 is reset to 1.8 meters to model the reﬁlling
of the head tube. The one-dimensional model reﬁlled 4 times in the ﬁrst 0.1
minute, again at 0.88, 1.7, 2.2, 2.4, 2.7, 2.8, 2.9 and 4.2 minutes. The reﬁlls
cause the spikes in the graph of node 1.

Node 2 is the ﬁrst node to become saturated as indicated by the curve
crossing the time axis (about 2 1/2 minutes) with the other nodes becoming
increasingly wetter (approaching zero). At 5 minutes the height of water in
the head tube was 1.12 meters. The soil water pressure head at node 2 is
0.56 m, at node 3, -0.78 m, at node 4, -7 .55 m and -19.99 m at node 5. All
other nodes are still at the initial condition of -20 meters. These values may
be estimated from Figure 11 where it is seen that nodes one and two are
saturated (above 0 m pressure head), node 3 is about -1 m and node 4 is about
-8 m. The data given here was taken from the program output. The slopes of
node 1, between reﬁlls, are of interest. The slopes appear shallow, and
similar, for time less than two minutes and time greater than three minutes.

Between two and three minutes the model reﬁlled often and the slopes are

39

consequently much steeper. This occured just before node 2 became
saturated and transpires because as node 2 becomes wetter the hydraulic
conductivity increases rapidly an causes greater amounts of water to be

drawn into the model.

4O

 

 

 

 

 

b-

wator potential (m)
:‘s

 

m r I I 7 7
0 1 2 3 4 5

 

 

Figure 11. Results, One-dimensional analysis

 

41

2.6 Summary

The one-dimensional model was written to study the parameter
determination and interaction as well as the general solution method. It was
particularly useful in determining the method to be used in calculating the
parameters k and 1.. Once this was accomplished an axisymmetric model was
studied in greater detail. The one-dimensional computer program was
written in Turbo Pascal for a personal computer (PC) and is given in
Appendix C. Comprehensive results and conclusions are given for the

axisymmetric model.

3. AXISYIVIMETRIC MODEL

The model of the velocity permeameter (VP) in two dimensions is
analogous to the one-dimensional model. The three-dimensional model was
written using axisymmetric elements. Axisymmetric elements are created by
rotating a two dimensional element around an axis of symmetry, in this case
the z-axis (Figure 12). Rectangular elements were used with the nodes
labeled i, j, k and m beginning in the lower left-hand corner and proceeding

counter-clockwise.

Using axisymmetric elements reduces the ﬂexibility of the regional
geometry, as the soil inhomogeneity and anisotropy can not vary in the y
direction relative to r or 2 although it may vary from element to element. The
solution of the axisymmetric problem is much less complex than is a full

three-dimensional solution.

3.1 Governing Equations

Writing the axisymmetric problem in cylindrical coordinates and noting

that the problem is independent of a rotation about the z-axis gives

Kmaw 6 BW 6‘11
r 6_r + a—Z[K(\P)(3?+l)] - x; (3.1)

 

a 6‘?

— K —

ar[ E 3 6r

where r is the radial distance from the z-axis and the other variables are as

deﬁned in Chapter 2. The saturated axisymmetric equation is

6%! K,a\ll 62w
Krzr—z + 7-67 + Kz-éZ—z '3 0 (3-2)

42

43

 

 

/\Z OXlS

 

 

 

 

 

 

/ p” \

 

/\

 

 

Figure 12. Cross-section of the axisymmetric element

 

44

3.2 Element Matrices

The element matrices for the axisymmetric model are determined from

the volume integrals as given in Segerlind (1984)

[s] - [ViBﬂDHBJdV (3.3)
[c]- f MNi’iNJdA (3.4)

and
01- [WNW (3.5)

The solution method for the axisymmetric model was the same as discussed
for the one-dimensional model (Chapter 2). The element stiffness matrix in
cylindrical coordinates is written in two parts, the radial term, r, and the

vertical term, 2
[S] = [5,] + [5.] (3.3)

The matrix in the radial direction, [3,] was

 

 

 

i 2 - 2 — 1 1'
237k,a — 2 2 1 - 1
[3'] " 6b — 1 1 2 - 2 (3‘4)
1 — 1 - 2 2‘
and the matrix, [9,] in the lateral direction was
’ (;+R,-) 7 -7 '(FTRiy
7 r +R- - 7+R- -7
[s.] - g r .) ‘ ,> .. (3.5)
30 -r -(7+Rj) (7+Rj) r
-(7+R,.) -7 7 (7+R,)

 

 

The parameters R,- and R]. are the shortest and longest distances, respectively,
to the element sides parallel to the z-axis. The variable a is half the element

width, b is half the element length and 7 = (R,- +R,-)/2.

45

The lumped formulation of capacitance matrix must be used when the
central difference method of solution is employed (Segerlind, 1984). The

lumped capacitance matrix for the axisymmetric element is

'2(F+R,.) 0 0 0 1
“M 0 2(F+RJ-) O 0
[c] ' 6 0 0 2(‘F+R,.) 0 (3'6)
0 0 0 2(7+R,.)

 

 

The parameters k and 1. were calculated as described for the one—dimensional

element (Chapter 2) using data from the graphs shown in Figures 8 and 9.

The equation for determining At for the rectangular element was given

by

AA

< —-——4k(1 _ 9) (3.7)

At

where A is the cross-sectional area of the element and 0 = 0.5 for the central

difference method.

3.3 The Model
3.3.1 The Physical Model

A cross-section of the axisymmetric model is shown in Figure 13. There
are 80 elements, each 20 mm x 20 mm, and 102 nodes. The core has two
elements in the x-direction and three in the -z direction (a total of six
elements in the core). Nodes 1, 12 and 23 are ﬁxed at 1.8 meters, the water
surface level in the tube, as the analysis begins. The core circumference is
modeled by a double set of nodes at x = 40 millimeters, with no element
between them. Figure 14 shows a plan view of the x-y plane for this model.
Since the z axis was a reﬂective boundary each element forms a torus around

the z axis, except the ﬁrst column (elements one to ten) which were

46

cylindrical (or torus with inside radius of zero).

47

 

 

2 0X55

/

60

80

100

180

140

160

 

(mm)

-20

 

-4O

\1
\

 

-60

\\
\

 

-80

 

-100

 

-120

 

~140

 

-160

 

-180

 

‘200

 

 

 

 

 

 

 

 

 

 

7‘ OXiS

(mm)

 

Figure 13. Axisymmetric finite element model

 

49

3.3.2 The Computer Model

Data was input from a ﬁle and initial calculations of water content and
element volume were determined, and then the time dependent solution was
begun. The ﬁrst portion of the time dependent solution was the saturated
analysis. The node above and below a saturated node in the z-direction were
ﬁxed and the arithmetic mean of the soil water pressure head of the two ﬁxed
nodes was assumed to be the soil water pressure head for the saturated node.
This procedure was performed on each saturated node except the ﬁxed nodes:
1, 12 and 23.

The next step in the model was the unsaturated solution. All saturated
nodes were considered ﬁxed during the unsaturated portion. The global
matrices were redeﬁned at each time step using the WP at the nodes to
calculate k and A. The [A] and [P] matrices were calculated using the central
difference method and then [A] was decomposed into upper triangular form.
Backwards substitution was used to determine the pressure head at each

node with both solution methods.
The backward difference method (0 = 1) was used to ﬁnd [A] and [P]

when the initial soil water pressure head was less than -20 meters, and the
central difference method was used in all other cases. The numerical solution
using the central difference method caused negative numbers on the [P]
matrix diagonal when the initial soil water pressure head was below -20

meters which caused an oscillating solution.

Appendix A gives the computer program written to run the VP analysis
with rectangular elements. The program was originally written with Turbo

Pascal for a personal computer (PC) but was altered to run on a Unix to

speed up processing.

50

3.4 Boundary Condition Determination

The nodes in the VP core at the ground surface were initially set to 1.8
meters pressure head. These were ﬁxed during the ﬁnite element analysis.
The other boundary conditions were insulated boundaries which did not
require any input data (g- - 0%:- - 0 across the axes, at the model extents and
at r = 40 mm, the radial edge of the cylinder). The drop in pressure head had
to be calculated after each iteration. This was accomplished by determining
the total ﬂow in the model for the iteration and dividing by the pressure head
tube cross sectional area (over which the flow had to occur) to give the drop in
pressure head, dh, in the pressure head tube. The drop was subtracted to
give the new boundary condition and this value remained ﬁxed during the

next ﬁnite element iteration.

The ﬁnite element analysis calculates the pressure head at each node of
an element. The corresponding water content was determined based on the
pressure head from the USDA Soil Conservation Service curves (as shown in
Figure 9, with actual data in Appendix D, and the program, Appendix A) for a
given soil. A representative volume for each node was calculated by
determining the torus area (since this was an axisymmetric problem) with
the inner radius equal to the node coordinate less 1/2 the element width and
the outer radius equal to the node coordinate plus 1/2 the element width.
The volume was determined by multiplying by the element depth. The
dashed lines in Figure 15 indicate the volume calculated to be represented by
each node shown. The solid lines indicate the elements. The nodes on the
boundaries were multiplied by 1/2 the element depth since there were only

two elements affected. The interior corner nodes affected an area of m2 since

51

they were on the axis of symmetry, and the exterior corner nodes had the
inner radius as previously described and an outer radius equal to the model

extent.

52

 

 

node 1

 

1
1

node 18

 

 

 

 

 

 

node 23

 

 

 

 

 

 

 

 

 

 

 

 

 

1
1
l
1
' H
1
l
|

 

 

 

 

 

 

 

J
l
elemerft 1 element
node 8 l
I node 13 node 24
I
1 Y — 1 _
element 2
i elemetc 12
node3 l l
node 14 node 85

 

Figure 15. Schematic to calculate element flow

 

53

3.5 Apparent Hydraulic Conductivity

The apparent hydraulic conductivity, K was calculated as described by

app)
Merva (1979) for each set of input parameters. The volume ﬂow was
determined by multiplying the nodal moisture content by the representative
nodal volume, as described in the previous section. The velocity in the core

was then calculated using

 

Q
vel =- 3.8
nr2*dt ( )

where Q was the volume ﬂow for the time step, r was the core radius and dt
was the time step. The height of water above node 1 in the core, (pressure
head), and velocity g) values were summed during each run in order to
perform a linear regression analysis to determine the slope of the line, 2%. Soil

water pressure head was the independent variable and velocity the

dependent variable. Multiplication of the slope, %, by the core length, 8, gave
Kw.
dv
Kapp = is (3.9)

The programming of the regression analysis is found in Appendix A under
the axisymmetric computer program procedure "FLOW".

3.6 Model with Triangular Elements

The axisymmetric model can utilize triangular elements as well as the
rectangular elements described previously. Triangular elements are useful
when studing irregular model conﬁgurations or when it is desirable to change
element size. In order to study the saturated front movement as it left the
VP core a model was initially written which used 2 mm x 3 mm elements in a

90 mm x 90 mm grid. This gave a total of 1350 elements and 1432 nodes.

54

Elements with an area of 6 mm’ necessitated a time step of 0.001 seconds
giving 300,000 iterations for a ﬁve minute run. A run of 90,000 iterations (90
seconds) with this model on the MSU Case Center computer never was
completed. An application was made, and accepted, to run on the N CSA
Cray-2. Although the Cray was much faster the time needed to run a model
with 1432 nodes (requiring 1432 simultaneous solutions each time step) was

not available for this research.

The area at the bottom of the cylinder whas the area of interest and
required a ﬁne grid, however, at the model extents the grid could be quite
course. Triangular elements may be used to change the size of elements in
ﬁnite element analysis and were tested to determine their feasibility. The
model with triangular elements is seen in Figure 16. Small elements, 2 mm x
4 mm, are concentrated around the bottom and outside of the VP. Triangular
elements were used to increase element size to 20 mm x 20 mm as was used

in the control data ﬁle. This model has 1 10 elements and 101 nodes.

55

 

80 40 60 80 100 120 (mm)

P axis

-20

-4o

-1

 

2 0X55

 

 

 

Figure 16. Axisymmetric finite element model, triangular elements

56

3.7 Element Matrices, Triangular Elements

The triangular element matrices are similar to the rectangular element
matrices except they are 3 x 3’8 (rather than 4 x 4’3) since there are only

three nodes per element. The stiffness matrix was given by

[S] = [8,] +152] (3-10)
The matrix in the radial direction, [9,] was (Segerlind, 1984)
b? 12,1), b,b,

1

12,17, b? bjb, (3.11)

1

12,1), bjb, bf

21t7k,
4A

 

[5,] =

and the matrix, [3,] in the vertical direction was (Segerlind, 1984)

2

_ c,- c,cj cick

21trkz 2 2

4A cicj cj cjc,‘ (3.1 )
2

cc, cjck c,

l

 

[5.] =

where 7, k, and k2 are as previously deﬁned for the rectangular element, A is
the triangle element area and the b and c variables depend on the distance to

each node

bi=zj'zk ci=Rk'Rj
bj=Zk-Zi cjzRi-Rk (3.13)
bk=zi'zj ck=Rj'R'i

where Z,, Z,, Zk, m, m, and R, are the radial and vertical distances to each node

(i,j,k) as shown in Figure 17.

57

 

 

 

 

 

 

 

 

 

\/

 

Figure 17. Triangular element coordinates

 

58

The lumped form of the capacitance matrix was given as (Segerlind,

personal correspondence)

(37+Ri) 0 0
[c1 - 3%3 0 (37.12,) 0 (3.14)
0 0 (374-12,)

The equation for determining At for the triangular element was given by

2M

A‘ < 9k(1-0)

(3.15)

where A is the element area, 0 = 0.5 and k and A are for the saturated values.

The computer subroutines which calculated the triangular element
matrices may be found in Appendix E. These subroutines were added to the
basic VP computer program (Appendix A). A check was made to determine if
the element was triangular or rectangular, the proper subroutine used and
the element contributions were put in the global matrices [K] and [C]. The

solution then proceeded as discussed for the rectangular elements.
3.8 Axisymmetric Analysis, Rectangular Elements

The model using axisymmetric rectangular elements proved appropriate
for observing water ﬂow through unsaturated soil. The data ﬁle for this
model had an element size of 20 millimeters x 20 millimeters. A time step of
0.1 seconds was used (based on equation 3.7 and using saturated values for
soil parameters). The extent of the wetting proﬁles may be seen Figure 18.
The water potential at each node was initially set at -20 meters (except for
the ﬁxed nodes) and this water potential is still seen at a radial distance of
140 millimeters (100 millimeters outside the core) and a vertical depth of 180

millimeters. Water potentials ranging from -1 m to -10 m are also shown.

59

Gravitational flow predominates as seen by the depth of the isoheads below
the core. The isohead of ~10 m radiates about 50 mm from the core and is 80

mm below the bottom of the core.

The wetted front moved vertically downward (one-dimensionally) until it
reached the bottom of the core (node 4 at a depth of 60 mm). Flow became
three dimensional as the wetted front moved out of the bottom of the core and

the saturated front moved more slowly than with one dimensional ﬂow.

Figure 19 shows results from the axisymmetric VP model: the drop in
head of node 1 and the wetting of nodes 2 - 5. The three dimensional graph is
similar to the graph of the one dimensional analysis. Node 2 becomes
saturated at about 1 minute, node 3 at 3 minutes and node 4, at the center of
the model and the bottom of the core, at about 12.5 minutes. This was
similar to the result seen with the one-dimensional model when it is recalled
that the one-dimensional element length was 50 mm and the axisymmetric
element depth was 20 mm. Comparing node 2 from the 1D model (2 = -50
mm), which saturated at about 2.5 minutes, with node 3 from the
axisymmetric model (2 = -60 mm), which saturated at about 3 minutes, shows
similar results with the two models. One signiﬁcant difference, however, is
that the axisymmetric model appeared to be reﬁlled with some regularity, as
shown by the spikes in the curve of node 1, whereas in the one-dimensional

model ﬂow into the soil slowed after 3 minutes.

60

 

80 40 60 80 100 120 140 160 (an)

1 l l l l l \
A /

\

~80 n

 

 

 

 

 

 

 

 

Figure 18. Soil water pressure head distribution, axisymmetric

model

 

water pressure (m)

61

 

 

 

time (min)

Figure 19. Results, axisymmetric analysis,

62

3.9 Axisymmetric Model - Triangular Elements

The axisymmetric model with triangular elements was run for 15
minutes with a time step of 0.01 second. The time-step was one-tenth the
time-step used for the rectangular model because of the small triangular
elements that were used. Figure 20 shows the equipotential curves for this
conﬁguration. The model is smaller than the axisymmetric model utilizing
only rectangular elements to keep the number of nodes and elements to a
workable number. Comparing the equipotential curves for the model with
triangular elements with the model using rectangular elements (Figure 18)
shows that the equipotential lines are similar for the two types of elements.
The models will be referred to as triangular and rectangular to differentiate
between the cross-sections of each. The elements in both models are
axisymmetric.

Figure 21 shows results from the axisymmetric triangular VP model: the
drop in head of node 1 and the wetting of nodes 2 - 5. The wetting proﬁles
are similar to those graphed from the 3D results, except that node 3 crosses
the x axis earlier in the triangular element model. The difference may be due
to the affect of the concentration of triangular elements beginning at the

depth of node 3 and the effect they have on the numerical calculations.

Notable differences in the model using axisymmetric rectangular
elements and the model using triangular elements were observed. One
dissimilarity was the frequency of reﬁlls. The model with axisymmetric
elements reﬁlled 11 times in 15 minutes while the model with triangular
elements reﬁlled only 6 times (see Figures 19 and 21). Another variation in
the models was the time needed to saturate the VP cylinder. The cylinder in
the axisymmetric model became saturated after 10 minutes but the cylinder

with the triangular elements did not saturate by 30 minutes run time. After

63

a run of 10 minutes the sample cup of the axisymmetric model was saturated
and the total volume inﬂow of water was 3.08 E-4 m3. The model with
trianular elements had a smaller total volume inﬂow of 2.06 E-4 m3 after 15
minutes. Since the input parameters were the same for both models these
differences are suspected to be due to the element size difference in the two
models or the difference in the total number of elements (80 elements versus
110 elements) or a combination of both. To use the model with triangular
elements for comparison to the axisymmetric model it would be necessary to
increase the rate of inﬂow. The following section (4.4) discusses changes in
the input parameters which affect the rate of inﬂow.

Appendix F gives the soil water pressure head at each node at 15 and 30
minute run time and isohead plot at 15 minutes for the model using
triangular elements. These plots show the shape and extent of the wetted
front. It should be noted in this data that some nodes appear to become dryer
during the analysis. This effect is not physically possible and appears to be
an anomaly due to the methodology in calculating the hydraulic conductivity
for the element matrices.

The model with triangular (axisymmetric) elements was developed to
study the movement of the saturated front outside the core. Since the core
did not saturate in an acceptable run time this model was abandoned.
Knowledge gained in further study of the axisymmetric element model would
be useful in determining adjustments to the triangular element model that
might make it suitable for further study. The model may be useful for
studying the saturated front movement as it is, but the rate that the core

saturates would need to be increased.

64

 

40 60 80 100 120 (an)
1 \
7

r- QXIS

 

2 0X15

 

 

 

 

Figure 20. Soil water pressure head distribution, triangular model

65

 

 

 

 

 

 

 

 

pressure head (cm)

a
O
_—L

D
i
(I

I

 

 

 

 

time (min)

Figure 21. Wetting of nodes 2-6, triangular elements

4. RESULTS AND CONCLUSIONS

The axisymmetric VP model was useful in monitoring water ﬂow
through soils with different parameters and for studying different head
tube/core size conﬁgurations. The 1D model was essential in determining the
calculation method to be used for the parameters k and A. The subroutines
which allowed analysis of triangular elements did not prove as useful as was

anticipated.

4.1 Versatility of the Axisymmetric Model

The versatility of the axisymmetric model of the velocity permeameter
was examined by altering six parameters. The parameters that were varied
were: the initial soil water pressure head; the ‘P—k curve used; the ‘P - 9 curve
used; the value of r which describes the anisotropy between the horizontal k
value and the vertical k value (the k value in the y direction was assumed
constant for the axisymmetric element); the core diameter and the head tube
diameter. Heterogenity (the spatial variation of k) was not investigated but
could be studied by assigning different soil parameters on an elemental basis.
Table 1 gives the parameter conditions for each analysis. The input ﬁle name
is given in column 1 of Table 1 with each run’s input parameters given in the
corresponding row. The ﬁrst set of parameters listed was designated the
’control’ data set. The values chosen to act as the control data set were
somewhat arbitrarily chosen, particularly the ‘P—k and ‘11-9 curves. The
control parameters were necessary, however, in order to have a point of
comparison while varying the six parameters. Throughout the remainder of
this work these values will be referred to as the control values or control

parameters. Each set of parameters make up a hypothetical soil type.

66

67

Table 1. Input parameter conditions

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

initial k 6 core head tube time step
ﬁle head :- curve curve mm W r (see)
(um) (mm) (mm)
control ~20 1.0 1 1 80 7.1 0.1
[01 ~10 1.0 l 1 80 7.1 0.1
IC2 ~50 1.0 1 1 80 7.1 0.1
103 ~100 1.0 1 1 80 7 .1 0.1
K1 ~20 1.0 2 1 80 7 .1 0.01
K2 ~20 1.0 3 1 80 7 .1 0.6
WC1 ~20 1.0 1 2 80 7 .1 0.1
WC2 ~20 1.0 1 3 80 7.1 0.06
Krl ~20 10.0 1 1 80 7.1 0.1
Kr2 ~20 0.1 1 1 80 7.1 0.1
001 ~20 1.0 1 1 40 7.1 0.1
CD2 ~20 1.0 1 1 120 7.1 0.1
HTDl ~20 1.0 l 1 80 12.7 0.1
HTD2 ~20 1.0 1 1 80 3.2 0.05
where:
head - soil water pressure head
r deﬁnes the relationship:
It. - "’9
1: curve, deﬁnes the W-k relationship used:
1 ~ Pachappa sandy loam (k_ - 20mm/hr)
2 ~ India loam (It... - Zahara/Iv)
3 ~ Chino clay (I:- - Zonal/hr)
6 curve, deﬁnes the we relationship used:
1 ~ Zigenfuss soil (6 - 46 % at saturation)
2 ~ Capac soil (8 a 36.8 % at saturation)
3 ~ Lenawee soil (9 s 43.4 % at saturation)
core, head tube ~ measurements in rmlh' 'meters

 

 

 

The k curve was obtained from data by Gardner (1958), as discussed in
Chapter 2. The 6 curve information was obtained from SCS Soil Water
Retention Curves from the Midwest National Technical Center (1988).
Analysis may have been facilitated by using the same soil types for the k

curves that were used for the W-G curves had this information been available.

The time step was 0.1 second in most cases but had to be decreased for

ﬁles K1, WC2 and HTD2. In the ﬁrst two cases the smaller time step was

68

necessitated because the parameters in the equation governing the choice of
the time step (eq. 3.7) had changed. This was also the reason that a larger
time step could be used for ﬁle K2. The time step for ﬁle HTD2 had to be
smaller due to the physical restrictions. The head tube was so small,
compared to the soil core, that the drop in head exceeded allowable

boundaries with a time step of 0.1 seconds.

The control value of initial pressure head was ~20 m and the soil was
considered isotropic. The ‘P-k curve approximated that of a Pachappa sandy
loam with a saturated value of 20 nun/hr. The ‘ILG curve approximated that
of a Zigenfuss soil. A sample cup (also referred to as the core) diameter of 80
mm with a head tube diameter of 7.1 mm was used (Table 1). Figures 18 and
19 (Chapter 3) give results from the VP analysis using the control

parameters.

Table 2 shows the time for the center of the core (node 4, at 60 m,
depth) to saturate and the total ﬂow volume after 5 and 10 minutes
simulation time. The center node of the core became saturated before any
other nodes at the same depth. The author’s hypothesis is that this node
became saturated ﬁrst since it is on a reﬂective boundary and ﬂow was only
in the downward (vertical) direction. All ﬂow was initially downward, in a
one-dimensional manner, at the beginning of the analysis since all nodes at
the ground surface (2 = 0) in the core had a water potential of 2 meters with
all other nodes in the model at an initial value of ~20 meters. This resulted in
faster wetting of nodes towards the interior of the core, particularly the nodes
along the z axis. Once the wetted front reached the soil core bottom ﬂow
became three-dimensional as water began to move in a radial direction out of

the core.

69

Table 2. Core saturation and total ﬂow

 

volume inﬂow

volume inﬂow

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

core volume inﬂow
file saturated at saturation at 5 minutes at 10 minutes
(min) (* 10’ mm’) (* 10’ mm’) (“ 10’ m3)
control 10.0 3.08 1.79 2.80
101 6.0 2.41 2.09 3.78
IC2 24.0 3.82 1.09 1.93
IC3 42.9 4.10 0.97 1.28
K1 1.4 4.12 13.32 30.07
K2 80.0 2.39 0.38 0.59
WCl 3.4 2.53 2.73 4.94
W02 3.6 2.19 2.61 4.69
Krl "‘ * 2.22 3.81
Kr2 3.3 1.39 1.56 2.34
CD1 * * 0.51 0.83
CD2 5.8 3.96 3.62 6.00
HTDl 9.7 2.68 1.71 2.68
HTD2 11.6 3.17 1.80 2.81
"' indicates that the center node of the core did not saturate
after 60 minutes running time.

 

 

 

These results are, generally, what was expected and indicate that the
model is applicable for studying water movement through unsaturated soil,
different soil parameters and different VP equipment conﬁgurations.
Explanations of these results may not be immediately apparent, however.
Each alteration in the initial conditions (from control values) will be

individually discussed in the following sections.

The total ﬂow at the time of the core becomming saturated was within
30% of the value for the control parameters except in the case of ﬁles Krl and
Kr2. File Krl never reached saturation and Kr2 became saturated at a much

lower ﬂow rate due to the gravitational ﬂow predominance. It would appear

70

that a ﬂow of about 3.1 E~4 m3 is required for the core to saturate. This value
is higher for soils with drier initial pressure heads (I02 and IC3) and for soils
with higher saturated k values (K1). The slope of the W-G curve (from which
i. is determined) also affects the total ﬂow at saturation. Curves that are
more linear (W Cl and WC2) seem to lead to lower total ﬂow values when the

core saturates.

4.1.1 Core Saturation and Total Flow

4.1.1.1 Initial Pressure Head

The time needed to saturate the core (node 4) was longer for the
conﬁguration ﬁles with drier initial pressure heads. The hypothetical soil
with the highest initial pressure head, ICl, had node 4 saturated after 6
minutes, the control parameters, 10 minutes, with the drier initial pressure
head taking 24 minutes (ﬁle IC2) and 43 minutes (ﬁle IC3). The same result
is reﬂected in the volume inﬂow results. With an initial pressure head of
~100 m (ﬁle: IC3) the total ﬂow after ﬁve minutes was just over 1/2 of the ﬂow
from the control values (initial pressure head of ~20 m). At 10 minutes the
volume inﬂow for the hypothetical soil IC3 is less than 1/2 the volume inﬂow
for the control soil. The volume inﬂow at core saturation, however, was
greater for the hypothetical soil IC3 than for the control soil. The wetter
initial pressure head (ﬁle: 101) had a 16% increase in total ﬂow after 5

minutes and a 35% increase after 10 minutes over the control parameters.

The initial pressure head of the soil is inversely related to the time to
saturation of the soil core. The time for saturation increases as the initial

pressure head decreases. The volume inﬂow at a given simulation time is

71

directly related to the initial pressure head of the soil. The volume inﬂow
increases as the initial pressure head increases. This is a consequence of the

increased number of paths for ﬂow in the wetter soil (less air pockets, etc.).

4.1.1.2 Hydraulic Conductivity Curve

The ‘1' - k curve relationship chosen to determine k affects the analysis

profoundly. Three curves were used. The control k curve approximated that
of a Pachappa sandy loam adjusted to have a saturated conductivity of 20
mm/hr. The curve for ﬁle K1 corresponded to that of a Indio loam, adjusted
to use a saturated k of 200 mm/hr. The curve for ﬁle K2 corresponded to that
of a Chino clay, adjusted to use a saturated k of 2 mm/hr. Appendix D gives
the actual data points used for determination of the hydraulic conductivity

from a known water potential.

The total ﬂow for K1, the ﬁle with the highest saturated k, was 1.332 E6
mm3 water after 5 minutes. This is an unreliable result, however, since the
wetted front had reached the model limits. Recalling that the model limits
are either reﬂective or impermeable boundaries it becomes clear that when
any signiﬁcant change in pressure head at the model limits occurs the run
should be considered ﬁnished. A larger model would be needed to obtain
reliable data for this ﬁle at 5 or 10 minutes run time. The results for ﬁle K1
show that much more water ﬂows through soils with higher saturated k
values and, consequently, the core saturates much faster. The control
analysis had a total ﬂows of 1.79 E5 mm3 at 5 minutes and 2.8 E5 mm3 at 10
minutes. The ﬂow for ﬁle K2 was 3.78 E4 mm3 at 5 minutes and 5.88 E4 mm3
at 10 minutes.

The time to saturate node 4 is similarly affected. The interior node of

the core was saturated after just 1.5 minutes with ﬁle K1. The core was

72

saturated after 10 minutes for the control parameters and K2 took 80
minutes to saturate the interior node of the core. Clearly the hydraulic
conductivity curve chosen, particularly the saturated k value, affects this

model remarkably.
4.1.1.3 Water Content Curve

The data used for the determination of A. and water content were taken

from the USDA Soil Conservation Service (Midwest Technical Center) for 3
soils (top 230 mm). The control ‘1’ - 9 relationship approximated that of a
Zigenfuss soil, W01 used a Capac soil and W02 used a Lenawee soil.
Appendix D gives the data used for the 3 soils.

The different curves appear to affect both the total ﬂow, as well as time
to saturate the core. Both curves chosen increased the total ﬂow and
consequently decreased the time to saturate the core. The control values had
resultant total ﬂows of 1.8 E5 and 2.8 E5 mm3 (5 and 10 minutes), whereas
W01 had total ﬂows of 2.7 E5 and 4.9 E5 mm", and W02 had total flows of 2.6
E5 and 4.7 E5 mm’. The core saturated at 10 minutes for the analysis with
the control parameters, W01 saturated at 3.4 minutes and WC2 saturated at

3.6 minutes.
The results for W01 and W02 are similar because the slope of the ‘1’ - 6

curve is similar for these two soils. These three soils were chosen because of
the difference in their saturated volumetric water contents (the Zigenfuss soil
(control) had a saturated volumetric water content of 45.3%, the Capac soil
(WC 1) had a saturated volumetric water content of 37.2% and the Lenawee
soil (W02) had a saturated water of 43.5%). The ﬁnal water does not effect
the ﬁnite element analysis, but A, the slope of the ‘1’ - 6 curve that is used in

ﬁnite element analysis. The slope of the curves for W01 and W02 were

73

similar and more linear than that for the control parameters. This lead to
larger total ﬂow values and faster core saturation times. The data used to

calculate A and WC are found in Appendix D.

4.1.1.4 Heterogenity

The study of heterogenity, where k, is numerically related to k, by the

multiplication factor, r (Table 4.1), showed the expected results. When radial
ﬂow dominates (k, 10 times k,, ﬁle Krl) the core never saturated and the total
ﬂow was greater than that for the control parameters (2.2 E~4 ":3, 23.9%
increase at 5 minutes simulation time and 3.8 E~4 m3, 36.1% increase at 10
minutes simulation time). The ﬂow was greater due to the increased k value,
over that for the control analysis, but the core did not saturate due to the
strong radial inﬂuence which pulled water out of the core. Once the wetting
front escaped the core and ﬂow became axisymmetric the larger radial k
values resulted in considerable radial ﬂows, with little vertical ﬂow. The
wetted front moved out of the soil core so quickly that the soil never became

saturated in the cylinder.

Conversely when k, = 0.1 * kz (ﬁle Kr2) ﬂow is preeominately vertical and

the volume inﬂow was decreased by 13% at 5 minutes and 16.5% at 10
minutes simulation time. The core became saturated in 3.3 minutes, much
faster than the 10 minutes for the control parameters, because of the

predominately vertical ﬂow.
4.1.1.5 Core Diameter

Changes in the core diameter produced total ﬂows roughly associated to
the cross sectional area of the core since this was the source of water for the

model. The core diameter of ﬁle CD1 was 25% the area of the control core.

74

At 5 minutes the ﬂow for CD1 was 511 mm3: or 28.4% of the ﬂow for the
control analysis. At 10 minutes the ﬂow for CD1 was 29.6% of the ﬂow for
the control analysis. Similarly the area of CD2 was 225% the area of the
control core diameter and the total ﬂow was 201.7% at 5 minutes and 214.5%

at 10 minutes.

The core never saturated for ﬁle CD1, the small diameter core. This ﬁle
was run for 150 minutes and still node 4 did not saturate. This size core did
not seem adequate to put enough water in the soil to saturate the core. This
may not be a problem, however, as a steady state condition appeared to occur

and will be discussed in Section 4.4.2 (apparent hydraulic conductivity).

4.1.1.6 Head Tube Diameter

The head tube diameter had no effect on the total ﬂow. This was
expected. After 5 minutes the total ﬂow for HTD1, control and HTD2 ﬁles
was 171 cm’, 179.4 cm3 and 179.8 cm’, respectively. After 10 minutes the
total flow was 269 cm’, 280 cm3 and 281 cm’, respectively. The time to
saturate the core was 9.7 minutes for HTD1, 10.0 minutes for the control
parameters and 11.6 minutes for HTD2 (the smallest head tube). The
differences in the time to saturate the core may be signiﬁcant but is probably

due to the ability of the head tube to supply ample water to the model.

4.1.2 Apparent Hydraulic Conductivity

Calculating the apparent hydraulic conductivity, kw, allows scrutiny of

the model and comparison with Merva’s (1979) ﬁeld analysis of the VP. The
apparent hydraulic conductivity, kw, was determined from a linear
regression analysis of head and velocity to give dv/dh which was multiplied by
the core length to give k,” (discussed also in Chapter 3). The time step was,
generally, 0.1 seconds but values of velocity as a function of head to be used

75

in the linear regression were taken every 6 seconds. This was done to more
closely duplicate ﬁeld readings and to smooth some of the irregularities that
the elemental model presented. The head and velocity data had to be taken
more often in a few cases (ﬁle K1 for example) because the core reﬁlled so
quickly.

The apparent hydraulic conductivity values were calculated before each
reﬁll. The results of this analysis on the ﬁles listed in Table 1 are shown in
Table 3 for two to three reﬁlls after the core saturated. The coefﬁcient of
linear regression is also given. The regression coefﬁcient was best just after
the core saturated. Much lower regression coefficients occured while nodes
became wetted. The numerical inﬂuence of unsaturated parameters,
particularly k, may have resulted in a nonlinear dv/dh relationship. After the
core saturated the model appeared to reach a quasi-steady state condition,
when the k,” values shown in Table 3 were recorded. The dv/dh relationship

became more nonlinear as nodes below the core became signiﬁcantly wet.

The saturated hydraulic conductivity used to determine k for all ﬁles,
except K1 and K2, was 20 mm/hr. File K1 had a saturated k of 200 mm/hr
and K2 had a saturated k of 2 mm/hr. The apparent hydraulic conductivity,
kw, values are generally smaller than the saturated value because although
the sample cup becomes saturated the soil under the cup does not, which
slows the movement of water through the model. The apparent hydraulic
conductivity values are generally within 50% of the saturated value excluding
ﬁles 102 and 103 (initial pressure head was less than ~20 meters) and ﬁle
Kr2 (radial ﬂow one-tenth of the vertical ﬂow). The soils that were initially
drier than the control value took signiﬁcantly longer to saturate (24 minutes
and 43 minutes as opposed to 10 minutes) and a steady state condition may

have become established with a k,” value that was lower than the saturated

Table 3. Apparent hydraulic conductivity

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

file time to saturate core Kappansnt (1’)
(minutes) (mm/hr)
11 (83.8%)
control 10.025 27 (87.9%)
11 (‘)
10 (82.7%)
IC1 5.9769 20 (90.8%)
8 (‘)
6.1 (94.5%)
IC2 23.9616 9.6 (89.3%)
16 (’)
5.2 (94.9%)
IC3 42.9043 6.4 (')
13 (‘)
16 (96.2%)
W01 3.4034 19 (95.4%)
16 (96.8%)
WC2 3.5851 20 (81.9%)
58 (‘)
K1 1.4492 69 (81.6%)
140 (92.3%)
1.3 (79.5%)
K2 79.95 2.5 (‘)
1.0 (‘)
last 3 the :
Krl core not sat at 60 min. 20 (91.5%)
27 (99.2%)
28 (99.1%)
6.1 (.)
Kr2 3.3451 8 (‘)
17 (94.3%)
HTDi 9.6833 1.5 (94.8%)
13 (93.1%)
HTD2 11.550 10 (97.7%)
14 (97.5%)
CD1 150 13 (98.7%)
(core not saturated) 13 (98.4%)
5.8 (‘)
CD2 5.760 9.8 (‘)
12.0 (‘)
* denotes a correlation coefﬁcient < 80%

 

value due to the inﬂuence of the unsaturated k value. File Kr2 also had

lower k,” value than was expected. The predominantly vertical ﬂow caused

the sample cup to be saturated quickly (3.3 minutes).

77

During the analysis nodes just below the saturated front seemed to
reach a critical pressure head (around -5 m) when they would cause a
signiﬁcant jump in the total ﬂow for the iteration which, in turn, affected the
velocity calculation and the drop in head in the head tube. This effect was
particularly noticable in the initial stages of analysis (at the beginning of the
run) which is why the km, values are given only after the core saturated. The
saturated hydraulic conductivity ~ apparent hydraulic conductivity
relationship merits further study.

4.2 Conclusions

The axisymmetric model proved useful in comparison of the effects of
varying soil parameters on the ﬂow of water through a soil. The VP core
diameter and head tube diameter choice also affect ﬂow rates and the model

may be used to estimate proper VP conﬁguration.

The triangular element model may be useful in monitoring the saturated
front movement. Another k curve (perhaps the India loam curve) would
speed up the analysis and make this model more practical. Further research
would be necessary to determine the usefullness of the triangular element in

the VP model.

Anisotropic ﬂow in a soil is complex. The method chosen to determine
the elemental hydraulic conductivity is critical to the success of the analysis
using the VP ﬁnite element model. Many other factors aﬁ‘ect the extent and
proﬁles of both the saturated and wetted fronts. The VP model is a useful
tool for studying water movement through the soil under the velocity
permeameter. The model increases our understanding of how water moves

through the soil under the VP, in manner and extent.

5. ADDITIONAL REMARKS

Some observations were made during the course of this research which,

while not germane to the objectives of this research, merit consideration.
5.1 Modelling Technique

A different model conﬁguration might help in the study of the movement
of water outside the core. The model utilizing triangular axisymmetric

elements might be useful, as has been discussed, in Chapter 3.

Another change which might be considered is using a double node at the
radial edge (outside) of the bottom of the core. Early models in this research
used the double node but the ﬁnal model was changed so as not to have
discontinuity at the radial edge of the core. However, after the change in the
model was made it was observed that the core had saturated faster with a
single node at the radial edge, and the node at the radial edge had a pressure
head closer to other nodes at the bottom of the soil core. The approach may
give smoother isoheads and more realistic (compared to ﬁeld data) results,
although the justiﬁcation of using a single node at the radial edge of the soil

core may be unclear.

5.2 Determination of k and A

The time-dependent soil parameters k and it had to be calculated

between each time step so the element matrices could be determined. The
elemental value had to be determined from the known water potentials at the
nodes. All four element nodes (axisymmetric rectangular element) might

78

79

have different water potentials so ﬁnding a representative k or i. for the
element presented great difficulties. The extreme variability of k for the
water potentials involved in this research made proper selection of the k and
1. curves imperative. The work of Gardner (1960) was chosen as a model for
the ‘I’ - k curves as he presented these curves over the desired range of ‘I’.

Lambda, 1., was calculated from actual data.

Both exponential functions and power functions (such as those used by
J abro & Fritton (1990) and Unlu (1990)) were used in the model. Coefﬁcients
were used that approximated Gardner’s curve for a Pachappa sandy loam as
closely as possible. Results from both methods were disappointing due to
ﬂuctuations in the water potential at the nodes. Nodes would begin wetting
and then get drier, a physically unlikely condition. With a smaller range of

water potentials the functions appeared useful.
The best results were achieved by using a 19 point ‘I’ — k curve and full

logarithmic interpolation between points. This method was also used for
calculating the 9 for a given ‘I’ and i. (slope of the ‘I’ - 9 curve), with the data
from the USDA Soil Conservation Service and semi-logarithmic interpolation.

Data for both calculations (k and A) are given in Appendix D.
5.3 Water Content and Isoheads

The total volume of the soil in the model was 1.6085 E7 mm’. The initial

pressure head of the control model was 20.9 % giving a volume of water of
3.362 E6 mm’. The volume of soil in the core was 3.02 E5 mm3 and at
saturation the pressure head is 45.3% (control parameters) with 1.37 E5 mm’
water in the core. The total inﬂow, when the core was saturated, was 3.08 E5
mm3 water indicating that more than half the inﬂow had escaped the core
(1.71 E5 mm3 water). In fact, the wetted front had moved to a depth of 140

80

mm (80 mm below the core) and radially 100 mm from the z axis (40 mm
outside the core). The isoheads (Figure 18 and Appendix B) show this more
clearly. Appendix B gives run outputs (for all ﬁles) of the water potential at
each node after the core saturates as well as a plot of some isoheads. These
plots show the shape of the wetted front for each set of parameters and are
useful for visual comparison of the extent of the wetted fronts. It should be
noted that that some nodes appear to become dryer during the analysis. This
effect is not physically possible and appears to be an anomaly due to the
methodology in calculating the hydraulic conductivity for the element

matrices.

Flow was one-dimensional while the wetted front was inside the core.
Once the wetted front reached the bottom of the core ﬂow became
axisymmetric. The node at the center of the model (node 4) is on a reﬂective
boundary (2 axis) so ﬂow must be directed downwards (gravitational) at that
point and isoheads must be perpendicular at that boundary throughout the
analysis. Nodes that are not on the reﬂective boundary, however, experience
ﬂow that is both radial as well as gravitational. Therefore as the wetted
front moves out of the core ﬂow has both gravitational and matric
components. The result is that the node at the center of the core became

saturated before nodes away from the center in all ﬁles.

The isoheads shown in Appendix B for the ﬁles are quite similar to the
control parameter isoheads. A notable exception is when the soil is
anisotropic (k varies with direction of measurement). Figures B9 and B10
show considerable variation for these parameters when compared to each
other or the control parameter isoheads. For a radial k that is ten times the

value of the vertical k the wetted front extends farther radially than the

81

isoheads for the control values. For a radial k that is one-tenth the value of
the vertical k the wetted front is severely limited in the radial direction, and

similar to the control isohead in vertical depth.

APPENDICES

APPENDIX A

Am'symmetric Computer Program (VP)

The com uter program was written in Pascal and run on a Unix system.
Procedure [INDATA read the input data from a ﬁle, procedure FEM ran the
ﬁnite element analysis, FLOW_CALC calculated the ﬂow velocity and drop in
head between iterations and RESULTS printed water potential results.

AAA LAAAA;;AAAAA-A AAAAAAAAAAAAA ALAAA

rogram vpt(input,output);
for the mainframe only - it won’t compile with turbo pascal}
type dlsiz = arra [1..200] of real;
smsiz = array 1..5] of real;
d38iz = arrayll..60] of inte er;
nlsiz = arrayll..200,1..4] o integer;
nmsiz = arrayIl..200] of integer;
nssiz = arrayIl..4] of inte er:
dZsiz = arrayIl..200,1..60 of real;
(note ~ if you change the array size be sure to change
the initialization in fnt_elemi
var data__hp: text;
data_name, save_name: acked array[1..15] of char;
out_name, gph_name, he _name: acked array[l..15] of char;
Wp_name: packed array[ 1..15] 0 char;

outﬁle: text,
graph: text;
ead: text;
wp: text;
ques: char;
title: packed array[1..40] of char;
phi. rf.
wcp, volel,
xc, yc: dlsiz;
dt , xav: smsiz;
nmtl, _n _array: nmsiz;
sat_nd_array,nel_sat: nmsiz;
nel: nlsiz;
fm: dlsiz;
c,k: dZsiz;
iptl, iteration,
iwt, i_main, i_siz,
j_main, kw,

k_choice, wc_choice,
num__x_core, num_y_core,
num_x_nodes, num_y_nodes,
, nbw, ncoef,
nd , ne, nsteps,
np, writ_mult,
num_fx_nds, num_sat_nds,
prnt.a,prnt_b,prnt_c,prnt_d,
prnt_e,prnt_f,prnt_g,prnt_h,
pmt_i, run_term: integer;
aa, ahtube, bb,
core_radius, core_length,
delta, dhtube,
gelta_q_total,

0".

length, lamda,
PL.

r, r1, r2, reﬁll,
sat_dist,

sum_x, sum_x2,
sumJ, sum_y?.,
auany.

sum_n,
run_time,

82

83

 

thta, timeavg,
total, unﬁw,
{dﬁnitiondthein parameters

~ a Wefmtment of the problem

:lfmcfeqiutionshlsonumberofnodes)
ne-numherofelements

nood- numberofsetsofequationcoeﬁcients
maxmumofﬁve

iptl-integercontrellingtheoutput
(4- dobuswtion)

kw-flag whichallowsprintmrtd'dataatvarims
pointsinthep

program
(kw- 0 cmitwrite,kw-1 write)
ahtube-areacftheheadtube
dhtube- diameterdtheheadtube
the number ofthe following sets ofparameters must

dtﬁl~ mam roperty that multiplies the time
derivativega 1-4)

q[i1~ constantcodﬁcientinthediﬁ'.eq.nodal
coordinate '

coordinates must be m numerical sequence relative

to node numbers

element data
nmt.ll~ t n m the coeﬁi

~integer equation cient set

nelIn,1]~ numerical value of node 1
nelln,2]~ numerical value of nodej
Martian ofthe variables read by Jintval
invl- ~integer controlling the input of the initial values

1 ~ input a node at a time

in
Mumbrswhichdonotchangewith
ibla]tirr:e(terrninatewithaserio)
theta~ ~thethetavalueusedinthesinglestepmethod
0 ~euler’sforwarddiﬁ'erencemethod
1/2-centraldifferencemethod

delta-thetime
run_time~ totaltimethemodelhasnm
~integercontrellingthetypeofanahrtical
solution (1 for the velocity permeameter problem
~numberoftimestsps
iwt-integercontrellingthemtputofthecalculated
valuestoGPHHdat valuesarelﬁ‘mtedeveryiwtsteps
writ_nrult~ contmlsthewritetoO datandtothescreen
ndhc- number of element sides with a derivative boundary
condition
r~ theratiocyddxaiuemonductivityinthey
direch'on is r’conductivity m the x direction
nelIn,1]~ numericalvalueofnodei
nelln,~4] ~numericalvalueofnodei
nel[n,3] ~numsrical value of node
nelln,4]~ numericalvalueofnodem
nelln,4]issstequaltoseroforthe
triangularel

numx_core~ numberd’telanentsinthecore
mthexdirection

mam_y_ core-numberdelementsinthecore
mtheydirection

mmwxelun numberofelemntsinthexdirection

num:y_ elem- mimberofelemntsintheydirection

r” m.”02?0:10(3;’i¢?0130n Om ‘ 3 3 gégzggzi C C C C

thiamhsmitineeitherreadstheinitialvahies

ddnitiond‘thevariablesreadbymval
~thevaluecltheinitialsoilmoistureisassignedtoall
} nodeathentheﬁxednodesaregiventheirvalue

VII

 

84

r: . integer; .
bep'n ’
iﬂiptl >- 4) then writeln(outﬁle,’ entering intval,
’np - ’=.np 6);
madln(dd.a _hp,head);
fu- r :- 1 to up do
M] =- heed:
readln(data ,
until (i I1); -hp'i’phi[in

-AAAAAAAAA AAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAA-AAAAA‘AAAAAA AAAAAAAA-

‘vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv vv vvvvvvvvvvvvvvvvvvvvvv’

eeeeeeeeegeeeeieeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
thisnbroutine andintvalreadsintheinputdata
inputdthetitlecard andcontrolparameters

 

coast idnn - 50;
var
na: nssiz;

i1; ii. : hm

basin
readln(data_ ,np,ne,r;coef, ,kw) ,
nts and the nodal coordinates}

readln(data:hp,core radius,core length);
ahtube :3 pr ‘(dhtuFe‘dhtube)/4-;
writeln(outﬁle);
writeln(outﬁle,’ r: k x - ’,r:0:3,’ ‘ h}; element size: 2 x 2 cm2’);
for 1.: 1 to ncoef do
readln(data hp, dtti].q[i]);

intval( ,iptl,p)
writelnIiPmrtﬁle outﬁlefthe initial soil moisture ll .phil2]:0:1,’ an');
readln(data hp,thta,delta);
writeln(outﬁle, ,thta:6: 2,’ delta: ’,’delta:0:2, secondd);
(input the number of time steps, nsteps,
and the write control, iwt (for nsults»
readln(data hpmateps ,nstepsj wt);
readlnfdataﬁk_ choice,wc _choice);m
readlnda ,num_x_core,num um_ x _nodeqnum nodes' ,
writalp(outﬁle,’ diameter head tu J
’cm. andthediameterdthe cors- ’,nuinx_:,’core‘40 cm');
readlnﬁa?_ _hrz,writ_ rmrlt,t h, ppztntbgprnt __b,prnt c,prnt_ d,prnt_ e,
.13"! P”
if (hw> 3') than-L

writeln(outﬁle,’ nodal coordinated);
writeln(outﬁle,’ node x y');

and:
bi. - Itonpdo
madam _:hPJcﬁD

fa'i - 1tonpdo
read(data_h [i]),
(outputoftﬁequatroncod‘icients)
(cutputolthenodal ' tes
rﬂhw>3)then

fai. .
Maladies? 4,ac[i]: 15: Gycﬁ]: 16: 6);
echo print (1' the element nodal data}
writeln(outfile,’ element data');
writeln(mrtfila.’nel nmtl node mrmbers’);

and:
nid:- 0;
forhk :-1tonedo

begin
"£61381.“ hpmnmﬂlkklmdlmllmomlmollnﬁlmllnﬂl);
if (kw > 3) then
writeln(outﬁle,n. '4,nmtl[th:6, nelln,1]: 7,nel[n,2]:6,
ndln.3]£.nelln.4]£);
iﬂ(n-1) <> nid) then
muldanfdg’elemt ’,n:4,’ not in sequence');
:- n;

(inputthemrmbersofthenodeswithvaluesthatremain
constantdm-ingtheﬁniteelemanalyais}

85

i :I 0;
repeat

Marika. hpﬁ-nd_ arrayﬁl);
gnfloﬁ_ nd _Wfi] <- 0);
repeat

ii:- sii+1;
until((fx__ndarray +1]<- 0)or(iimod6= 0));

untiufx nd" _arF-yl”' +1 <- 0),
um_?x_ nds > orthen

um_ I: _nds>idnn)then

write(«rtﬁle,’ number d boundary conditions enceeds’);
writeln(«rtﬁle,’ the allowed numberof’ ,idnn:6);

and:
and; .
{analysis of the node numbers
initialization of a check vector}
(creation and initialization d’ the a
nb vector calculation of the bandwidth}
w :2 0;

forhkm 1tonedo

for-i. . 1to4do
.nelIkkjl'
for-i. . 1to3do

':_-i+1;
j :niito4do

-:abs(na[i]-nsﬁl)

iﬂnb- 0) then wrrteln(«rtﬁle, element’,k.k:3,
’hastwonodeswiththesamenodenumber’);

iﬂnb>nbw)then
nbw:-rib;
and:
and;
and:

nbw:;-nbw+1

if(nbw>60)then

writ?n(autﬁle,’ w!nbws' ,nbw:0,’ (max-=60”;

calculation emacereqinred'}

if(kchoice-1)then
writan(«rtﬁle,-’hcurve-1 Pachappasandyloamﬂhat: 2cm/hr)');
if(h_chcics=2)then
writsln(outﬁle,’hcurve 2- lndioloamaisat- 20cm/hr)');
if(k_ choice-3)th«r
writeln(«rtﬁle,’ hams-3o Chinoclay(Ksat- 0.2cm/hr)’);
if(wcchoics-1)then
writeln(«rtﬁls,'wp/wc«irve-1-Zigsnﬁrsssoil');
if(wcchoics-2)then

wri,’teIn(«rtﬁle wp/wccurve-2- Capac soil');
if(wcchoics- 3)then

wri,’telh(«rtﬁle wp/wccurve-3- -Lenawessoil');

rocedureindata

{seeeeeeeeeeeeeeeeeeeeeeeeeeee ”0......”O...”O”OOOO”O”O”O0.0000}

fumeueeeeeeeeeeeeeoeeeeeeeeeeseeeeeeeeeeseeee

thisprecedurepu'intethewaterpotentialsatthenodestoﬁleouttxt

 

vet 1|.me

eounter,it: integer;
basin

ﬁg),

e«mta- :- num_y {used fordo }
' JJcIa]: 8: 1.xc[b] :8 1,8ch :1,xc[d]:8:l,xc[e]:8:1,
xdﬂﬁdxkl 8: mm: 8: 1.xcﬁl: 8: 1);

86

writdn(«rtﬁls);
for it :s 1 to cormter do

5,5211%“ 0) then
if (yclal 21.10.“ 0) then

mxﬂmda]:2z ',1,’ ’ hila]. ':8 1, hilb] :.8 1, hilc]:
11111181831). hi1:1:8:1pi:m::8 1.p' 8 1.115111 :8:

ama+1;
b:—b+1;
c:-c+1;
d:-d+ 1;
e:-a+1;
f:-f+ 1;

g‘; 111‘
i and.i+l;

writeln(«rtﬁle);

en ts
{000;00.‘0000000OOOOOOOOIOOIOOOOOOOOIOQttttttttﬁttototttttttttt0......)

8:1,
1

eat elem;
.0...OOOOOO‘;OOOOOOOOOOOOOOOOOO.

procedure directs the saturated ﬁnite element analysis

5.5. hi.
vary“: .

 

forhk :- 1 to num_sat_nds do

I0;
for nltonumﬁ ndsdo
ifga‘t ndﬁarrayﬁkrsfx_ nd _arrayIkiDthen

if (ﬂag- 0) then
bep'n
minfﬁfmcmliﬂmn]
«id; .
and:
and; snoosdure sat elem}
0.. .O...‘.....‘.....‘O...... 0“.O‘.....".‘.0......‘0‘00...O.‘0‘.0..}

procedure we calc(hk:in r;varmc:real);
‘C..‘..‘0.“‘...OOOO0.0QOOOOOO...
thissubrmrtinecalculatesthewatercontent,wc
atanodefor}aknownwaterpotential,wp

i

wc_type.-array[1.14;]ofreal

wc[14]. - 0.112:
and:

] I00
$13340”;

wp[6] :- 403.3;

87

:.m'ﬁ'l 3.23::
wp[8] E: 61928;,

$11] 2. $166.0:

”[12] :- 40330.0;

wp[13] :. 45495.0;
{wpl14m1031320.0; (0 23 »
CID“ BIO - cm
if(wc_choics:p2)thenn

wc[1] :2 0.372;
wc[2] :- 0.370;
wc[3] :- 0.366;
wc[4] :- 0.362;
“[5] :I 0.353;
wc[6] :- 0.338;

wc[13] :- 0.137;
wc[14] :- 0.093;

and.
if (wc_choice I 3) then
(lenawee soil, ap horizon (0—23 cm)}

well] :11 0.436;
wc[2] :- 0.433;
“[3] :I 0.429;
wc[4] :- 0.426;
wc[6] :- 0.416;
we[6] :2 0.4;
wc[7] :- 0.384;
wc[8] :11 0.368;
we[9] :- 0.333;
wc[10] :2 0.296;
wc[11] :8 0.247;
wc[12] :- 0.216;
wc[13] :- 0.199;
wc[14] :- 0.144;
and:

. the ordinate values - water content - vohrmstric}
wp - the abcissa values - water potential- in cm}
: 1033 cur/bar was used}
node :- phiIhh];
“(9048 < 'PlMD then
begin

writeln(«rtﬁle,’ value (',nodezl2z2,

') less than least absima value: node - ’,hk:0);
writsln(outfrle,’ value (’,node:12:2,

') le- than least absissa value: node - ’,hk:0);
writeln(«1tﬁle,’ the soil can not be this dryl');

lateto theta- ate tet
if(node<-10mn get w rcon n}

for i_wc :- 2 to 13 do
if (node <s i_wc]) and (node > wpfi wc+1]) then
me :- (( _wel-wcli wc+1])‘ln(node7wP[i_wc+1])/ln(wp[i_wc]
lwpﬁ_wc+ 11))+ wc-[i_wc+1];

else
me:- well];
' rocedurewc calc}
{OOOCOOOOOOOOOOOOOOCOOOOO... 0..O....000‘3000OOOOOOOOOOOOOOOOO0.0.0.0}

' r;var lam l:real;var ﬂagzinteger);

OOOOOOOOOOOOOOOO‘OOOOOO 00.00.... -

this «rbroutine calculates the al of the
head-water confent relationship,

we type - arrayll..14] of real;
var i'lam: integer;
, me: real;
we: we_type;

88

yr: we_type;
' ' ho ' (0.23 )
‘W 95.2”“ m”

wen: :- 0.463;

 

wc[14]. 'I 0.112;

'9“) 0. 0‘

"£21:- 40. 33;
wp[4] :- .51 .66;
wp[6] :I 403. 3;
wp[6]. :- .206 .'6,
”{7}; :s .340. 89;

8 a -619 .;8
wp[9] :s .1033. 0'
wp[10] :- -2066. 0;
wp[ll] ;. -5165. o;

12] :I 40330. 0;

wp[13] :- 46496.0;

[14] :I 403300. 0;
lapse soil, ap horizon (0-22 cm)}

if (we_ l”ﬁchoicsm 2)then

begin
we[1] :- 0.372;
wc[2] :- 0.370;
we[3] :- 0.366;
we[4]. 'I 0. 362;

wc[9] :I 0.269;

“[1013- 0 .229;
we[11]: I 0.181;
we[12] :I 0.161;
wc[13] :- 0.137;
“[14] :I 0.093;

{Lenawee soil, ap horizon (023 cm)}
1’ (we_ choice- 3)tha1

well] :I 0.435;
"[2] :I 0.433;
we[3] :I 0.429;
we[4] :11 0.426;
wc[6] :- 0.416;
wc[6] :- 0.4;
we[7] :- 0.384;
wc[8] :I 0.358;
wc[9] :11 0.333;
we[10]. 'I 0.296;
wc[11] 'I 0247;
we[12]: I 0.216;
wc[13] :11 0.199;
we[14] :- 0.144;

- the ordinate values- water content- volumetric}
(note-ﬂ” abci-avaluee owaterpotsntial- incrn}
:1033 cm/bar was used}

23:“: "° mm,

if(nodns4p < wp[14l) then

write(«1tﬁle,’ value lessthan least abeissa value. ');
writsln(«1tﬁle,’ nodd- ’,.node4 ',’0 wp[14] I ’,wp[14]: 0' 2);

89

write(«1tﬁle,’ the soil can not be this dryl');
writeln(«1tﬁle,’errer found m inter-p’);
1f (nede4 < 40.33) then
for 1_ lam :- 2 to 13 do
if «1 wp[i_ laml) and (node4 > wp[i_ lam+1]) then
begin

me :- ((wc[i lam]-wc[i lam+1])‘ln(node4/wp[i lam+1])l
ln(wp[i__ T [i_1am+1])) + wc[i_ lam+1;
lam _.l '11 ((mc- wc[1_lam+1])/(node4- wp[i_ lam+1]))

and:
and:
end
else
lam_l :- 0.0;

_ .. -1---A-AAA.A-A--A...---- -..-.-A........AA-A-.....1--.---.Ai--i.-
‘ -vvvvvv- vvvvvvvvvvaVVVVVJ v‘vvvvvvvvvv .. -vvvvvv- Viv Vtvvvv '1}
an 1 Vt?” Om“.
’ ‘ O

f”.~.”‘” ”0”.”OOOO”O”OO

this procedure «15:11.15. the hydraulic conductivity

type I arrayﬂ..19] of real;
teger

3:14} :I 420',
”[1513 .500;
wp[16] :- .7oo'
‘1)[17]. 'I .900;
”33}; ' 2330
(fo‘i'Ii’ ppaeandy loam: }

if(hP PdroieeI1)then

h_curve[1] :I 6. 66e-4
k _curve[2]. 'I 6. 66e-4'
k _.«nve{3] 'I 6. 66e-4
k _eurve{4] :- 6. 06e-4'
k _«1rve[6] :I 4. 04e-4'
k _airve[6] :I 3. 64e-4'
k _curve[7]. 'I 2. 62e-4'
h _curveIB]. I 3. 64e-6
h _.arrveIQ] 'I 1. 01e-6;
k_cur've[10 :I 6.3e-7;

h curve[11: :I 3. 63e-7;
h_«1rve[12: :- 2.3e-7;

h_«uve[13: :- 1.01e-7;
k curveluj :- 6.3e-8;

h curve[16j :- 4.04e-8;
h_curve[16j :I 2.02e-8;
k curve[17]:I 1. 62e-8;
h _curve[18‘. :I 2. 02e-9;

 

(k choice- 2)tha1

90

k _WJ '8 6. 280-3;
k _cm-ve[6.] I 6e-3;
k _.cur'vel61 I 4.17e-3;
k _:cur-vel71 I 2. 780-3;
k _:cur'veIB] I 1.39e-3;
Ecru-v49]: I 2. 92e-6;
h_curve[10] :I 1. 39e-6;
k _curveﬂl] :- 6. 66e-6;
k _cur've(12] :- 2.92e-6;
k_curve[131:I 1.11e-6;
k _curve[14] :I 8. 33e-7;
k _cur've[16] :- 2.9297;
h_cmve[161:- 1 .94e-7;
k _curve[17] :I 1. 39e-7;
h_cunn[18] :I 8. 33e-9;
k cr1rve[191:I 6. 66e-10;

cley-eeturetedvelueof02un/hr}
if(k choice-3)then

k_.cm've[1] I 6. 66e—6;
k _:curvel21 I 6 .M
k _curveﬂ'l]: I 2. 78e—6;
h _:curveI4] I 2. 6e-6;
k _cur've[6]: I 1 6.7e-6;
h _.ameIG] I 1 .39e-6;
k _:curve[7] I 8. 33e-6;
k _:curveIB] I 2. 78e-6;
k_:curve[91 I 8. 33e-7;
k_curve[101:I 6. 66e-7;
k:curve(11] :I 2.92e-7;
k_cur've[121:I 2. 22e-7;
k _curveIl3] :- 8. 330-8;
k_cur-ve[141. I 6. 94e-8;
k _curve[161 :I 2.92e-8;
h_curve[161. I 2. 6e-8;
k _curve[17] :- 1. 39e-8;
k _curve[18] :- 2. 77e-9;
k cur-0419]:- 8. 33e-10;

{h_a1rve-the ordinete veluee - weter content volumetric
tehn from richer'de but multiplied by e oonetent
to meet reehﬁc eetureted conductivity veluee)
{wp -the ebcieee veluee - weter potentiel- 1n cm}

”halal,“ I 0) then
heed =' philnoltkl. 3]]

nN=-Phi[nellk1.4]];
312°“) -8)thenkl :Ik_curve[1]
for-i Phym 1t018do

if(heed<-_kxy1)end(heed> [i_hythhen
Ielp(( curveﬁ l/k_ curve _hxy+1])
“'95 l/wrili +11)
_kx,y+1])) + ln(k_ curveﬁ _hxy+1]));
Ikl ‘ 1';
I0) thar

ifwr'it:llr(oul:ﬁle,’ error-l, kl- 0, node ’,kl :0);
end; (edeorptionﬂowiegovernedbythe upper nodee}

{-O”M:O”O”O”O”O”O” ”COOOOOO” ”C”.”O”OOOO”O”OOOOO‘O”O)

mtx(ele1n:integer;kl,klx,lem_ lzr-eel);

{m “0...“. WO”O”O”O”O”

thieprﬁueaeeteethe‘bbdmetrioeqcendk

.m-fnud .41 .4;]olreel
uis-erreyl dreel;

ver ecrrr, eern,et,
ee: eeie;
«any: nix:

”a":

91

3.!

*3
.1,: g
3.

5

$151111}? '
105-23>?th
a '- a " v a
.w .H g.»

£13.99.
0 O 0'. I O
O
3"

I

khhnemh
.0 0'. O. 0.. O. C. I.
. 5: r

 

gfﬂiiléééiii

telnCentering glb_ mix, let e1em.');
element without the
civetive condition)
{retrievel ofnodel coordinetee end node number-e}
{0:29; 1 to 4 do

.i =- nelldunil:

:5] - '1;

1'6] :- 1'1

‘51. I 1.00;

a :- y[41- [11;
:' 8121.111].
er :- ee‘bb
=- uni + rum/2. o;
if (ebe(er) <I 0. 00001)then
begin

writeln(ouﬂ'1le,’the eree of element ’ ,elem:,4

’ie leee then 0. 00001');
en'iteln(ar.1tﬁle,’the node number-e ere in the wrong order’);
writeln(outfile,’ or the node form e etreightline

end;
{'mitielizetion: et[ ] rraetr'rx n31“. thembl }
- time epe art p em
m }

et[1,1] :I the:- + 3(1);
0611.21 :- rber;
“193] :' "hm
“[194] :3 “[191];
“291] :. ”'3
et[2,2] :- the:- + x[2];
$12,313 412.21;
$12.4 1:- Ir;
#3, 1] :I «but;
#13313 $12 31;
$13.31:- $12.21;
04.413 rber:
et[4,11:- et[1,4];
«14.21:I 4'5
“[‘131 :3 l'bll';
et[4 ,41 :- et[1,1];

M metrix, time depende t at}
c - n p
1ft! c-1umped formuletion}

for 1 :I .1 to 4 do

or-kdr-
if ((1ptl >I 4) end (elem I 1)) then

92

heg'n
forj:I1to4do

bean
:I(1‘rber'klx‘ ‘ 6‘1)
eemﬁj] (2"p1 1_"kl‘b‘et[1,j]/?6:'rmy( )1
061116.51" 8 (pi_ ‘er’dte‘lem _l’lmpﬁgm;

$1 :I pi_ ‘qe‘ee‘bb‘dIil/4;

whgoneluetoredmcolumn1

bytheminordiegoneleﬁonbw)

incolumns2- nbw. emeereusedtoﬁll

thecolumntonp,themetrixsizeienpxnbw}

{the ﬁerce veftorzfm, hes fixed nodel veluee}
or1 I to

basin

’5'6'1'3ufm561l‘dn
+ 1;

3333.11.44.

:InellelemJ];
J1: Ij6+1~ 6;
£61>0)en (ji<Inbw)then

(1:53:11: 8 8:155:51 + “[1,;1;

> n
gainfmtﬁle,’ ’j6 > up in glb_ mtx! ﬁx it!’);
and:
' Wye b
eeeeeeeeemeeeeeeeeeeee ‘0‘.”ﬁ£“m}.‘.”.ﬁt.ﬁtﬁO..OOOOOOO..0.0
{ }

{.e.-222222.112.

thismbroutinemodifreetheglobel cepecitence end
Mnemmetriceewhenthereerenodelveluesthet
remeinwnshntwithﬁmgusingdeletiond’mwsendcolumne.

M}
ver i_mod,j_mod,jm_mod,
m:nnd,n _mod. integer;

‘iﬂiptl >- 4) then

wr-iteln(outﬁle,’ executing modify);
writeln(outﬁle,’ up I ’:,np 6, nbw I ’,nbwz6)

and:
modify dkmetrr <1de endcolumns
‘ for-if; :Iltonmeet_nde rows }

begin
mod .Ieet__nderrey[i mod];
jr'r'rrrod': °JI modjl;
{thieeememrtt'herewfortheﬁxednode t
theﬁntcdummwhichisectuelbrthedi'emJ
fol-jm_:Imod 2tonbwdo

m mod :1. j _mod+jm_mod-1;
iffm_ mod<I np)then
mod hi mod ,
{161111 1:13:13- ﬁnlxrmmodl-kl mod-1m. 1‘? [i_ ]

{thiseeroeouttheeohmnfortheﬁxednodm but

Remember-the ' elsmeheup thecolumnseo

celeuletionsmust doneeteboute46degreeengle}
if(n_mod>0)then

21::Lmod2nI fm[n_::r(r)r:;d1-k[n_mod,im_mod1‘phi[j_ mod];
nmod :_Inmod-1

and;

and:

:1

93

if matab;

”“0.“0”0”.“.~.

this ubroutine gem-ates the a and p matrices
fa- the n'ngle step methods using the equations

.(,).Ic[,1+(d.1umhwkt. 1
1x. ) Ic[, Hdelta)’(l thetb‘kl. l

 

 

var i_ma mat: integer,

aFmat, _mat,tac,th: real;

if (iptl > 4) then wﬁtelnhutfile,’ entering matab');
a_:tamatItb delta;

b:mat :-(1thta)‘delta;

fori_mat:I1tonpdo

begin
forj mat:I1tonbwdo

:I r t t
a. ‘°3.iir$.£’,_ $3,
{a matrix}
oIi matj_mat] :I tc+a_mat‘tk;

{p matrix}- _
H1_mst,)_mat] :- to-b_mat‘tk
end;

AA-A‘.AA‘A‘

{magnum

thispmoeduredeoompesthekandcmatrioesintouppertriangular
formusing}gauuianeliminatio n

var1dnn, Jdmc,

m. 1:271:11]; n1. n91: integer.

1f (iptl 1ptl>I 4) then writeln(outfile,’ entering danpbd');

np1:Inp-;1
fori dine:- 1tonp1do

begin
.I i_dmc+nbw-1;
$1» i> up) than mi:= up;
mk:I l"nbw;1
if(( i _dmc+1) < nbw) thenmk :I np-i_dn1c+1;

{c.1312 :- 1;tomkdo
:- nd+k2
if(o[i dmc,1]I 0) then

'ﬁWG-ltﬁb,’ cl’.i_dmc: 0111] ' .cﬁ dmc, 13])
oﬁ W] 3' cﬁ.dm€.k21-cﬁ_dmc.nll‘cti_dmc.'_nk1/cﬁ dmc 1];

A-AAA-A‘ -AnAAAA-AAAAAAAAAAAAh‘AAAAAAAA AAAAAAAAAAAAAAAA‘AA‘AAAAAAAA-

f”.~.“.”.”.”..

thispmoedurepuformsthematﬁxmltiplication
var i k{mu-id, E mbd,

 

sum:m
iﬂiptll >- 4) then writeln(outﬁle,’ entering multbd');

94

 

if (h2: > 0)then
:21“; :. a11;1+k[k2,i_ mtbd]“phi[k2];
and;
sad;
.;fd’ILmtbd] I a1m+k[i_mtbd,1]‘ph1[1_mtbd],

“-4: AAAAAAAAAAAAAAAAAAAAAAAA (2 m 921?st AAAAAAAAAAAAAAAAAAAAAAAA
‘l """""""""""""""""""""""""""""""""" I
919995131! 91114; .......

 

‘vvvvvvvvvvvvvvvvvvvvvv

this plocedune decompes the global force vector, f, and solves
the system in! equations using backwast substitution

var flag,
i _slv, j_ slv, k2,
1:91". mi.
no. npl 3. 11': integer.
91111: real;

hunP-l;
{decompositionofthecolumn ecto rﬂ)}
for1__:slv-1tonp1do V 1'

I1 slv+nbw-
mi>npithenmi=- up;
:i1;-_llv+

v:I1;

jslv :Injtomido

:I l_ siv+1
iii-1v] :- IflLllVI-(CIL Ilv,l_ slvl‘rﬂi_ alvl/di _slV.11);

and;
{backwaxdaibstitution for determination of phi[ ]}

fhﬂglg =-1 mildnpdls

333' 5;.

'I

91%

i_slv:Inp-k2;

3(1_ slnv+nhw-1) > np) then mi:= np-i_ slv+1;
am: I 0.0;
for J slv. I 2 to mi do

n_ slv: I 1 alw-
aim :- afmwei' _slvj :slvl’philn_ slv];

_slv] :I (1ﬁi_slv]-a1m)/c[i_siv,1];

{End Finite Element Analysis - new heads (phi) found)
{the code following adds a saturated
nodetotheﬁxednodeset}
fa- slv :Iltonpdo
hi[i_ slv]> -10. 33) then

:IO
:8 m_ndl’
is“ 11:2?“ i:slv do
if(;at_n arrayﬂ_slv]Ij_ slv)thenflag:=1;
)then

95

 

 

aid;
whom-1M} AAAAAAAAAAAAAAAAAAAAAAAAA
................................................................ I
fmcedme fnt elem;
”.“O”.~“‘O”.”O”OO
this pmcedure directs the ﬁnite element analysis
vari fnt.i_fnt,elem_fnt,ﬂag: integer;
lam_l: real;
fori fnt :I 1 to 200 do
fm[i__ fnt]. I 0. 0;
{oci-fat :- 1 to 60 do
fnt]. -I 0. 0;
fut"31ml: - 0.0
and;
forelem_ﬁ1t:I 1 tonedo
hy(elem fnLkLklx)
lambda(elem _amfnt,l Lang)
glh_ mtx(elem_ fnt ,klﬂmlam_ l);
(calculates the element matricies
calculation of the element capacitance (c),
stiﬁ’nessa), matrices and the element
. fence vector (0}
matah;
nm time. I um_ time + delta;
multhd;
for 1 fnt: I toting
fnt]:I1rf[i t]+delta’fm[i _fnt].
alvbd;
99d: ..................... {(194131}! AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
1" ' ' ' " ' ' """"""""""""""""""""""""""" I

mcedureﬂcwcalc

eeeeeeeeeeeeaeeee'e
thisprocedmedeterminesthevolumedwaterintmduced inthelaet
iteration andperforms the reg1'eeeion anab'sis to determineK.

 

caldi ﬂow. [i_ n l);

wcpl- i_flow ﬂ'cg‘volzlYL flow] + voll;
writelﬂhead,’ initial water volume I’,v011:12:3);
tcd'law :I 0.0;
{initial set-up}

96

delta_wc:Iv012-voll,
Voila-V012,

:Itctflow+delta we,
hlz- '[1]:
dh:I ta_wc/ahtube;
fcriﬂowmltonumfx fxndsdo

nIdIarrqIL ﬂown. I philfx_ nd __an'ay[i flow]]- dh;

“(aptilﬁitﬂlof warning dh 0, change
m - < mch’

delta wc04’at,runt1me02, seconds’);

vel:- delta _wc/(pi_ ‘core _radius’core _radius‘delta);
if (iteration mod ith 0)then
writeln(he:l ”,run_ time/60: 0:4} ’,totﬂow:0:6,’ ’,hl:0:2,
v .

item-.1301. mod iwt - 0) then
begin

a1m_x :Ia1m_x+h1;
a1m_y:-a1m:y+vel;
a1m_x2-sumx2+(h1‘h1)
aim_y2 :Isum:y2+(vel’vel);
aim_xy :_Isumxy+(h1"vel);
writaliif amh}.iIlli]:8.1’ ’phi[12]'8:1’ ’ hi[3]::81’ ’phi[4]::81,
9p 9 D 9 0p 9 9
31:55]: ':8 1.’ ’.phi[6]=8=
if(r1_in time> 290) and (run time< 310) then
wri (graph,’ at ,run_ time/60::03,’ minutesthetotal ﬂowwas’ ,
totflow:0:,’1 cm3');
if(run t1me> 590) and (run _time<610) then
m teln(graph,’ at’ ,run__ time/60: 0: 3,’ minutes the total ﬂow was’,
totﬂow:0:,’1 cm3’);

1f(1teration mod writ_mult I 0) then

begin
writeln(autﬁle);
writeln(outfile, ’time: ’_,run time/60: 0: 2,’ minutes (itr. ’,
iteration: 0,’)’);
(autﬁle,’ head at each node (in centimeters):');
results;
writeln(outﬁle );
wnteln,’(mitﬁle the total inflow 1s ’,:totﬂow 10: 3,’ cm3');

$111“) > .10) and (run__ term - 0) then
bean

runterm:_Irunta-m+1;
w:iteln(outfil’e, node4iseaturated; at’ ,run_ t1me/60::,’04 min');
wr1 h,node4issaturated;at’_,mntimel60::,’04 min');

awn < 100.0) then

{ﬂyi'lmm‘ _hnndsdo

writeln-(ouﬁlen mvﬂlill?refillzz30,’ ’cm)at’ ,run_ time:5: 1,
'seconds (’,run _time/60: 7:2,’ min)’);

if(phi[4]>-10)thenrun_ term :Irun _term;+1

Wand line regress:
if(a11{n_n>2)then m ar on}

an :I sum_x2-((a1m_x‘a1m_x)lsum_n);
syy :- sum _y2 ((sum _y‘sum _y)/sum n);
my. I um_-1y ((sum x‘a1m_y)/sum_ n);
mat _:r I sly/aqua“ syy); {correlation coeﬁ'icient}
out '32 :a stat _m t_r;
ﬂat) :- sly/sax; b
stat a :I (sum _y/a1m_n) - (stat_ ‘(sum_x/sum_n));_
writeln '

(graph).
writeln(grap:)’ the head tube rd'illed at ’,run_ time/60: 0:2,’ minutes’);
write1n(graph,’ the linear coeﬁcients cf mansion ');
writeln(a-aph,’ for this drop arez') ;
writelnw-ph);
writeln(graph,’ a I ’,stat_a:0:6);
'ﬁtolnb'aph);
writeln(graph,’ b I ’ ,etat _b:0:10, ’ (slope, cm/ch);
writeln(graph),;
writeln h, r2I ,stat 1204’ or’,

stat_100::01,'%');

97

(grap ’ thenKa ; ’ .0. ’cm/hr‘);
wﬁtehm'aph’ m.ult by ntcouffﬁngfh)? ');
whmh):

aim_x:;I0
aim:y:;I0
m_:-x2 0;
aim:y2:I I0;
aim_xy:I0;
aim_n:I0;
and;

and:

rocedure ﬂow calc}
{M”.“.”m.“.“.”.0...3...”‘”..OC“....”..O..OOO..OOC‘0‘...)

{iuebeeeeegmheeeeeeeee

inpu}tsectionofthepma'am

 

 

data name: I 'vp.da
reseﬂdsta_ ,data name);
nadln(data_t1tl-),
out_ name :- ’out.txt’
3:31am name: I’mhtxt’;

:- 'hed.ut’;
wp_ name :I ’.wp txt’;
1ewrite(outf1le,out name);
writeln(mitﬁle,’Vefocity Permeameter Analysis ');
writeln(outf11e );
w1iteln(outﬁle,’ data ﬁle name: ’ ,title);
rewri ph,gph_ name);

1mmM saddled 3111):»)- .

rewn wp,wp name ;

writeln h,'Reaem1o'

writeln ap ’time totaljilnflaw head vel');
writeln(head, ’(min) (m3) (cm) (cm/ Deer);

mm . mm mm hil3 him 51');
( magiﬁim £1qu 1’ pm

{callproceduretoneaddatgcm}
“MFG;
iteration :80;
num_ sat nds :_Inumfx nds;
fori_ma1n :I1tonum E ndsdo
_nd_' i _mainl': '_I-h nd _array_[i main];
fori mainleto4do

_main] :I would _minn;
.yg;__main]:ch[nel[1,i_mainD;

total:- 0. 0;
u . (3(21- nun/2;
bb=-(y[ﬂ-y[1])/2:

fcri main. 'I 1tonpdo

:Ixcli _mainl-aa;
12:Ixc[i _main]+ +aa;
ifucﬁmsin] ulnp-D:thenr2=xc{i_ main];

98

volelﬁ_ main] :- pi_ ‘2‘hb‘((r2‘r2)-(1'1‘r1));

if (xcli min]: 0)then
velel[i_ main] := )2‘bb‘pi_ ‘12‘12;
if (ycli main]: 0) or (:5 _main]- ycInpl) then
$5653") .. W121i“) yin bow; minl) _length) th
I an a < core on
ifvolollignam]fumirtﬂel[i_ma1:1n]/2

ﬂow, calc;
writeln(«ﬁle );
writeln(«1tﬁle,’ (initial vah1ee)’);
writeln(«ﬁle);

run _ta-m :I 0;

um_ 8. :-g;

«nn I

«m32::I 0;

um_-:72:- 0;

«1m _xy :- 0;

«imn :- 0;
”.“ﬁnittiiiifiiCCCCiltiiithtiittﬁiﬁ‘
eolutien of the time dependent problem}
meat

iteration :I iteration + 1;

if (iteration mod writ _multI 0) then

writelnCita'ation ’,iteration: 0);

if(num_ eat _nde>num _fx _nde)then

eat elem;

ﬂow calc;
until-(iteration I natepe);
{ until (iteration I netepe) or (run WI 4);
} wr1teln(’ 11m completed. ﬁle: ’,title);
aave name. 'I ’aave’;
clue-(data hp,aave_ name);
cloae(h

ve _name);

 

h,aave _name);
cloae( ,aave _name);
cleee(wp,aave_ name);

Appendix B

Soil Water Pressure Head

Appendix B lists the soil water pressure heads just after the core
saturates and shows the isoheads for each set of parameters studied with the
axisymmetric rectangular element VP model.

The pressure head in some cases seems to decrease. The effect seemed
to be due the numerical problems resulting from the hydraulic conductivity
calculations.

99

100

MW

r:k_r= 1.0‘k_z; elementaiaez20mmx20mm
theinitialeoilwaterpmaaum headia -20.0m

theta: 0.6 deltaI0.1eeconda

diameter, headtuhe=7.1mm. andthe diameter, eomISOmm
kcurve-l-Pachappaaandyloammaatz20mm/hr)
‘P-ecurve-l-Zigenfuaaeoil

 

'1‘ahle Bl. fmanum head at each node- control parametere

time at aaturation: 11.0 minutea f l l I
headateachnocleﬁnmetere): _L thetotalvolumeinﬂowin3.07986 m3

         

 

 

 

 

 

 

 

 

 

 

 

 

 

z\r 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160. 0
(mm)
0.0 1.231 1.231 1.231 49.990 49.996 -20.009 4999.6 -2000.3 4999.2
-20.0 0.616 0.616 0.616 49.976 -20.004 -20.003 4999.4 -2000.3 4998.3
I -40.0 0.308 0.308 0.308 49.778 49.897 49.968 4999.3 2000.3 4998.3
-60.0 0.206 -0.194 -0.429 -0.910 4.718 48.676 4999.3 -2000.3 4998.3
-80.0 -0.609 -0.666 -0.748 4.021 -2.103 48.666 4999.3 -2000.3 4998.3
_ 400.0 -0.848 -0.943 4.138 4.696 -4 .268 49.216 4999.4 -2000.3 4998.3
420.0 4.616 4.860 -3.226 -9.937 47. 720 49.870 4999.4 -2000.3 4998.3
440.0 48.200 48.687 49.171 49.743 49.974 49.999 4999.4 -2000.3 4998.3
| 460.0 49.993 -20.016 -20.021 49.980 -20.004 -20.003 4999.4 -2000.4 4998.3
480.0 49.993 -20.019 -20.023 49.977 -20.004 20.003 4999.4 -2000.3 4998.3
-200. 0 49.993 -20. 011 -20.000 49.994 49.994 -20.010 -2000.0 -2000.0 4999.0
0 0 120 140 160 (nm)
' L 1 1 \
‘ /
I f‘ (1)05
-20 ._ I
|
l
-40 J I
l
-60 _.4 :
l
-80 _ l
I
l
‘100 _. I
l
_ I
180 _l |
|
-140 _ l
l
l
-160 _ l
I
—190 _ :
|
-800 _. _____________________ .1
(mm)
\/ z axis

 

Figure Bl. Control ieoheada

WM

r: k_r a 1.0 ‘ k_z; element nine: 20 mmx 20 mm
the initial aoil water pmanum head in 40.0 m

theta: 0.6 delta: 0.1 aecondn
diameter, head tube = 7.1 mm,

W-ecurve-l-Zigenﬁinnnoil

 
 
 

  

10

1

diameter, core a- 80 mm

kcurve-l-Pachappaeandyloammaata 20mm/hr)

 

Table B2. Preanum head at each node - lCl parametem
. r
time at aaturation: 6.0 minuten

 
 

I the total volume inﬂow in 2.412 E6 m’

 

 

 

 

 

 

 

 

 

 

 

 

 

. head at each node (in meters): 1 _

z\r 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0
(mm)

0.0 1.911 1.911 1.911 40.000 40.000 40.000 40.000 40.000 40.000
-20.0 0.967 0.967 0.967 -9.999 -9.998 40.002 -9.997 40.000 -9.996
-40.0 0.478 0.478 0.478 «9.884 -9.942 -9.987 -9.998 -9.990 -9.996
-60.0 0.319 0.127 -0.406 -0.961 .2118 -9.807 -9.998 -9.999 -9.996
-80.0 -0.611 -0.678 -0.763 4.044 ~2.469 -9.780 -9.998 -9.999 -9.996

400.0 -0.899 -0.983 4.146 4.630 4.416 -9.884 -9.998 9.999 -9.996
420.0 4.766 -2.146 -3.469 -7.016 -9.698 -9.983 -9.998 ~9.999 -9.996
440.0 -9.667 ~9.732 -9.854 -9.964 -9.996 40.003 ~9.998 -9.999 -9.996
460.0 -9.999 -9.997 40.001 40.000 -9.999 40.003 -9.998 -9.999 -9.996
480.0 -9.999 -9.997 40.001 -9.999 -9.999 40.003 -9.998 -9.999 -9.996
~200.0 -9.999 -9.996 -9.997 40.001 -9.998 40.001 -9.999 40.003 -9.996
120 140 160 (mm)
0’0 1 1 1 \
' /
I 1" axis
-20 _ I
I
|
-40 _ I
I
—60 e I
I
-80 _ I
I
|
- l 00 _, I
I
~1eo _ I
I
- l 4 0 _ I
I
- l 0 n I
|
-180 a I
I
-aoo 1L ____________________ _J
( )
nn 4/ z axis

 

Figure 32. [Cl iaoheadn

102

W

r:k_r- 1.0‘k_z; elementnixe:20mmx20mm
theinitialaoilwaterpmanum headin-60.0m

theta: 1.0 delta=0.1aecondn

diameter, headtube=7.1mm. diameter, com880mm
kcurve-l-Pachappananhloam(Kaat=20mm/hr)
cm've-l-Zigenfuanaoil

Table B3. Preanum head at each node - IC2 parameters

 

___..— _ _eiv ,,
I

time at saturation: 26.6 minutea I

headateachnodeﬁnmetern): I thetotalvolumeinﬂowia3.824E6 m’

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

z\r 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0
(mun)

0.0 1 .360 1 .360 1.360 -60.000 450.000 -60.000 -60.000 «60.000 -60.000
-20.0 0.676 0.676 0.676 -60.000 49 .999 -60.000 ~60.077 -60.003 49.911
40.0 0.338 0.338 0.338 48.686 49.247 49.767 60.061 -60.003 49.911
-60.0 0.226 -0.176 -0.4 14 -0.896 4 .743 40.386 -60.036 -60.003 49.911
-80.0 -0.498 -0.667 -0.739 -0.999 -2.018 40.821 60.033 -60.003 49.911

400.0 -0.836 -0.926 4 .147 4 .684 -3.293 43.820 -60.046 -60.003 49.911
420.0 4.666 4.803 -2.892 43.100 -36.417 48.219 60.066 -60.003 49.911
440.0 -39.663 41.268 44 .236 47.777 49.744 49.962 -60.076 -60.003 4.9911
460.0 49 .962 49.907 49.907 -60.023 49.976 60.016 -60.077 -60.003 49.911
480.0 -60.003 49.924 49.909 -60.034 49.974 60.013 -60.077 -60.004 49.912
-200.0 -60.001 49.939 49.988 60.07 7 450.000 450.022 -60.066 -60.013 ~60.000
0 0 100 120 140 160 (nm)
' 1 1 \
I 7
I r axis
- 20 I
I
I
- 4 0 I
I
-60 I
l
- 80 I
I
I
-100 I
I
-120 I
I
—140 e I
I
I
- 160 ._ I
I
-180 _ I
l
- 200 _. _____________________ _I
(run)
/ z o xis

 

Figure 33. IC2 1.0m

W

r:k_r= 1.0‘k_z; elementniae:20mmx20mm
the initial nail water pmanum head in 400.0 m

theta: 1.0 delta 8 0.1 was

diameter,headtuhe= 7.1mm,
kcurve-l -Pachappaaandyloam(Kaat= 20 mm/hr)
W-ecm-ve-l-Zigenﬁiaaaoil

ThlB4. Whead node - IC3 parametern
, , r 1
time at naturatian: 44.0 minutes

103

diameter, com s 80 m

I 7. 7 W . , . e.
I

 

 

 

 

 

 

 

 

head at each node (in metal-a): 1 total volume inﬂow in 4.102 E6 nun’
z \ r 0.0 $1.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0
(min)

0.0 1 .642 1 .642 1 .642 400.000 400.000 400.000 400.000 400.000 400.000
-20.0 0.821 0.821 0.821 ~99.780 -99.915 400.079 400.014 -99.994 -99 .999
40.0 0.411 0.411 0.4 1 1 -94 .661 ~97.034 -99.021 -99.963 -99.971 -99.992
-60.0 0.274 -0.147 -0.405 ~0.927 -2.087 -71.457 -99.870 -99.965 -99.985
-80.0 0.497 -0.560 -0 .754 -1 .030 ~2.43 1 -73.799 -99 .866 -99.963 -99.983

400.0 -0.868 -0.%3 -1 . 197 4 .745 -3.942 -82.370 -99.942 -99.961 -99.981
420.0 4.871 -2.180 -3.499 48.974 -58.552 ~94.358 400.014 -99.961 99.980
440.0 -74.028 -78.049 84.800 -93.676 -98.829 -99.901 400.077 -99.963 -99.981
-160.0 -99.821 99.944 -99.610 -99.801 -99.990 400.216 400.046 -99.967 -99.982
480.0 ~99.995 400.050 -99.625 -99.817 400.000 400.197 400.045 -99.972 -99 .985
I -200.0 -99.999 400.006 ~99 .849 -99.943 400.129 400.000 -99.923 -99.992 400.045
— — — - —————— ==l=u— — —- -- -
100 180 140 160 (mm)
0.0 \
4 g A 1
I /
I r o xis
- 20 ._ I
I
I
- 40 _ I
I
— so __ I
I
" 80 _ I
I
I
" l00 ._ I
I
— 120 _ I
I
- 140 .1 I
I
I
- 160 _ I
I
-190 A I
I
- 200 .4 _____________________ .1
(mm)
/ z o xis

 

 

Figure B4. IC3 iaoheadn

 

 

 

104

W

11k}: l.0‘k_z; elementsize:20mmx20mm
thein'nielsoilwsterpmreheedis -20.0m

theta- 0.5 delts=0.01seoonds

diemder, hesdtube=7.1mm, diameter, core=80mm
kcurve-2-Indiolosm(Kssts200mm/hr)
W-ecurve-l-Zigenfusssoil

Table 85. Pressure heed st eech node - K1 peremetere

 

 

 

 

 

 

 

 

 

 

 

I T 7 ,, ‘
tune' st saturation: 1.5 minutes
hesdst eechnodeﬁnmeters): J. thetotslvohmemﬂow' is4.1221'35 m’
8 \r 0.0 20.0 60.0 80.0 100.0 120.0 140.0 160.0
1.599 1.599 -20.000 -20.000 -20.000 -20.000 -20.000 -20.012
0.800 -20.000 49.971 ~20.018 49.994 49.995
0.400 49.865 49.890 49.974 49.988 49.992
0.267 -0.916 4.350 -6.183 49.881 49.992
-0.498 -0.995 4.429 -6.867 49.840 49.992
-0.822 4.285 4.706 40.378 49.893 49.992
4 .229 4 .88. 1 -6.169 44 .216 49.966 49.992
41.713 46.429 49.885 49.992 49.992
49.920 49.967 -20.015 49.995 49.992
-20.005 49.973 -20.019 49.995 49.992
~20.007 49.993 49.997
0 0 120 140 160 (mm)
' 1 1 1 \
I /
I r- axis
-20 _. I
I
|
" 40 _ 1
I
-60 _. :
l
'80 _ I
I
l
- 100 __. I
I
-120 _ :
I
- 140 ._ I
I
l
I
-180 _ :
I
'200 -I_ ____________________ .J
(nn)
4/ z axis

 

Figure 35. K1 isohesds

 

 

 

 

105

W

rzk_r= 1.0‘k_z; elernentsisez20mmx20mm

the initial soil water pressure head is -20.0 m

theta: 0.5 delta=1.0seeonds

diameter, headtuhe=7.1mm. diameter, oore=80mm
kcurve-3-Chinoclay(Ksat-2mmfhr)
W-ecur've-l-Zigenﬁrsseoil

Table 86. Pressure head at each node - K2 parameters
time at saturation: 79.9 minutes I

 

 

7<

 

 

 

 

 

 

 

 

 

 

head at each node (in meters): I the total volume inﬂow is 2.388 E6 m’
I z\r 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0
0.0 1.984 1 .984 1 .984 49.996 49.998 -20.005 -20.001 49.998 -20.001
~20.0 0.993 0.993 0.993 49.924 49.967 49.995 49.999 49.9% -20.004
-40.0 0.497 0.497 0.497 48.309 49.161 49.756 49.989 49.9% ~20.005
-60.0 0.0 0.141 -0.391 4.359 -5.316 48.634 49.972 49.9% -20.006
-80.0 ~0.546 0.659 -0 .986 4 .594 -6.052 48.207 49.965 49.996 -20.005
400.0 4.338 4.545 -2.147 4.143 -8.031 48.844 49.972 49.9% -20.005
420.0 -5.523 -6.438 -8.604 42.232 48.050 49.4% 49.987 49.9% -20.005
= 440.0 48.462 48.648 49.069 49.556 49.865 49.984 49.997 49.9% -20.005
460.0 49.970 49.969 49.992 49.994 -20.005 -20.009 ~20.002 49.9% -20.005
480.0 49.999 49.992 -20.010 -20.003 -20.009 -20.011 -20.002 49.9% -20.005
-200.0 49.999 49.994 49.996 -20.003 49.994 -20.011 -20.002 49.995 -20.004
7 7 7 7 777 7 7 7 — 7
0 0 120 140 160 (mm)
' 1 1 1 \
I /
I r o xis
_ 20 l
I
I
- 40 I
I
-60 I
I
- 80 I
I
I
-100 I
I
-120 ._ I
I I
7140 r I
I
I
- 1 6 0 _ I
|
-180 _ I
I
“300 _ _____________________ _I
(nn)

 

\/ z axis

Figure 86. K2isoheads

 

 

 

W

106

nk_r- 1.0 ‘k_z; elementsise:20x20 mm
the initial soil water presarre head is -20.0 m

theta . 0.6 delta = 0.1 seconds
diameter, head tube a 7.1 mm,

diameter, core = 80 mm

heurve-l-Pachappasantbloam(Ksat=20mm/hr)
W-chr've-2-Capacsoil

     

time at saturation: 4.6 minutes

 

 

 

 

 

 

 

head at each node (in meters):
1 \r 0.0 20.0 40.0
(mm)

0.0 1.626 1.626 1.626
-20.0 0.814 0.814 0.814
-40.0 0.407 0.407 0.407
-60.0 0.271 0.271 -0.296
-80.0 -0.286 -0.412 -0.61 1

400.0 -0.692 -0.762 -0.960
420.0 4.231 4.317 4.693
440.0 43.380 46.113 47.668
460.0 49.996 -20.002 -20.006
480.0 49.998 -20.006 -20.007
-200.0 49.998 -20.003 -20.000

 

 

 

100

180 140 160

___——_.__———__————_———————_—1—

-200 4.. ____________________ .J

 

\/ z axis

Figure B7. W01 isoheads

 

 

(rm)

r 0X15

 

 

 

W

r: k_r a 1.0 " i_z; element size: 20 x 20 mm
the initial soil water pressure head is ~20.0 m
theta = 0.6 delta = 0.06 seconds
diameter, head tube a: 7.1 mm,

107

diameter, core a: 80 mm

kcurve-l -Pachappasandyloam(Ksat= 20mm/hr)
W-chrve-3-Lsnawessoil

 
   
 
 

 

Table BS. Pressure head at each node - WC2 parameters
time at saturation: 4.0 minutes

I __

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

head at each node (in meters): I the total volume mﬂo‘ w is 2.193 E6 m’
— - -— 7 77 7 7
0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0
-20.000 .20_000 -20.000 49.991 49.993 49.997
49.989 -20.001 -20.001 49.984 49.985 49.992
49.934 49.971 49.991 49.987 49.987 49.991
-0.918 4.813 49.4% 49.986 49.986 49.991
4.040 -2.264 49.364 49.987 49.986 49.991
4.636 -5.534 49.697 49.987 49.986 49.991
40.432 48.443 49.971 49.986 49.986 49.991
49.881 49.999 49.999 49.986 49.986 49.991
49.991 -20.0(B -20.001 49.985 49.986 49.991
49.990 40.0% ~20.001 49.985 49.987 49.991
49.998 49.999 49.999 49.994 49.990 -20.001
20 60 100 120 140 160 (mm)
0.0 \
1 1
I /
I .
- 20 — I I" 0 XIS
I
|
- 40 _ I
I
-60 _ I
I
-80 _ I
I
I
- 100 ._ I
I
7120 _. I
l
'1 4 0 _ I
l
|
I
-180 _ I
|
‘8 00 ______________________ .1
(mm)
\/ z axis

 

108

W
Note: This file did not saturate so water potentials are shown at 10 minutes. The file was
run for 60 minutes, by which time the radial extents had been reached so that data would be
incorrect.

r: k_r a 10.0 ‘ k_a; elemsm size: 20 x 20 mm

the initial soil water presure head is -20.0 m

theta = 0.6 delta :2 0.1 seconds

diameter, head tube = 7 .1 mm, diameter, core - 80 mm

kcurve- 1 -Pachappa sancb'loam (Ksat = 20 mm/hr)

‘P- 8 curve - 1 - Zigenfuss soil

 

Table 39. Pressure head at each node - Krl parameters

.-..........:..~.. 1W...““*T " “ ' ‘
head at each node (in meters): | the total volume inﬂow is 3.811 E5 m3
» i H i , ‘

 

 

 

 

 

 

 

 

 

 

 

 

z \r 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0
(mm)

0.0 1.382 1.382 1.382 -20.000 -20.000 -20.000 49.996 49.998 -20.000
~20.0 0.692 0.692 0.692 -20.000 49.999 ~20.000 -20.001 49.997 49.991
-40.0 0.346 0.346 0.346 49.829 49.849 49.933 49.9% 49.906 49.980
-60.0 ~0.609 -0.637 -0.613 -0.822 4.026 4.346 -3.429 47.806 49.948
-80.0 4.086 4.106 4 .148 4 .174 4.374 4.804 -6.373 47.860 49.940

400.0 4.628 4.727 6226 6.196 -7 .633 42.096 44.930 49.369 49.970
420.0 49.690 49.709 49.740 49.766 49.889 49.860 49.966 49.964 49.990
440.0 49.998 -20.004 49.991 -20.004 49.993 ~20.002 -20.002 49.996 49.992
460.0 49.997 -20.000 49.989 -20.006 49.993 -20.002 -20.000 49.994 49.990
. 480.0 49.997 -20.000 49.989 -20.006 49.993 -20.002 -20.001 49.994 49.990
-200.0 49.993 49.986 49.998 -20.006 49.989 -20.001 49.996 49.997 -20.001

0,0

 

-20
-40
-60
-80
-100
-120.s
-140 _
-160 ._
-180 _

-200 _

(mm)

 

\/ z axis

 

Figure 89. Krl isoheads

 

109

W

r:k_r- 0.1‘k_z; elanentsise:20:20mm

the initial soil water presure head is -20.0 m

theta: 0.6 delta=0.1 seconds
diamder,headtube-7.1cm. diameter,coie-80mm
hm-l-Pachappasanlhloam(Ksat=20mm/hr)
V-Ocurve-l-Zigenfusssoil

able BIO. Pressure head at each node - Kr2 parameters
time at saturation: 4.0 minutes

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

headateachnodeﬁnmeters): thetotalvohimeinﬂewisl.392 E6 m’
8 \r 0.0 20.0 4 0.0 60.0 80.0 100.0
(mm)

0.0 1.490 1 .490 1 .490 40.000 40.000 40.000 -20.000 -20.000 40.005
-20.0 0.746 0.746 0.746 49.995 -20.007 49.997 -20.001 49.999 -20.002
40.0 0.373 0.373 0.373 49.979 -20.0(B 49.997 -20.002 49.998 -20.002
-60.0 0.249 0.249 ~0.84 1 45.097 40.005 49.997 -20.002 49.999 -20.002
-80.0 —0.303 ~0.427 4 .497 45.147 40.005 49.997 .20_002 49.999 40.002

400.0 -0.978 4.051 4.336 49.734 -20.002 49.997 -20.002 49.999 -20.002
420.0 -8.331 -9.37 3 43.333 49.902 40.006 49.997 ~20.002 49.999 40.002
440.0 49.981 49.987 49.989 49.994 40.007 49.997 -20.002 49.999 40.002
460.0 40.001 -20.004 49.997 49.995 -20.007 49.997 -20.002 49.999 -20.002
480.0 40.001 40.004 49.997 49.995 40.007 49.997 -20.002 49.999 -20.002
-200.0 -20.004 -20.000 49.999 49.998 49.998 40.000 49.998 49.999 -20.005
20 40 60 80 100 130 140 160 (mm)
0, 0 \
1 1 1 1 1 1
I /
I I" 0 XIS
- 20 _ I
I
I
- 40 4 I
I
-60 _ I
I
- 80 _. I
I
I
- 100 ._ I
I
- 120 ._ I
I
I
I
- 160 _. I
I
-180 _ I
I
“200 _u— ____________________ ..I
(nn)
\/ z a xis

 

Figure BIO. Kr2 isoheads

M

r: i_r :- 1.0 ’ h_z; element size: 20 x 20 mm
the initial soil water pressure head is -20.0 m

theta- 0.6 delta a: 0.1 seconds

diameter, head tube a 7.1 cm,

110

diameter, core a 40 mm

kairve-l-Pachappasandyloam(Ksat= 20mm/hr)
W-ecurve-l-Zigenﬁisssoil

  
  

time at saturation: 45.0 minutes
head at each node (in meters):

 

 

77““ 1 , , .
1 thetotalvolume inflowis3.106 E6 m’

 

  
  
  
  

 

 

 

 

 

 

 

 

 

 

 

 

0 0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0
1.169 1.169 -20.000 49.976 49.989 40.017 49.993 40.006 49972
-20. 0 0.680 0.680 40.061 49.902 -20.010 40.006 49.987 40.006 49936
40.0 0.290 0.290 48.953 49.117 49.668 49.832 49.969 -20.010 49936
-60.0 -O.286 -0.467 -0.809 4.106 4.997 46.111 49.948 -20.010 49936
-80.0 -0.664 -0.726 -0.856 4.176 -2.231 44.939 49.936 ~20.010 49936
400.0 -0.876 -0.932 4.081 4.379 -3.387 46.361 49.961 o20.009 49936
420.0 4.241 4.319 4.517 -2.362 -8.093 48.711 49.969 -20.010 49936
440.0 -3.441 4.382 -7.137 41.646 46.943 49.736 49.979 -20.010 49936
460.0 48.400 48.866 49.335 49.636 49.913 49.996 49.979 -20.010 49936
480.0 49.981 40.062 40.078 49.930 40.019 40.012 49.980 -20.010 49936
49.996 49.981 49.983 40.036 40.000 -20.000 49962
20 40 80 100 120 I40 160 (mm)
0,0 \
1 1 1 4 1 1
/ I /
r axis
-20 .4 I
core / I
-40 _/ I
/ I
-60 a I
l
—80 _. I
I
l
-100 _. I
l
-120 a I
I
-140 _ I
l
l
‘160 _. I
I
-180 _ I
I
-200 .41.. ____________________ ..I
(nn)
/ z axis

 

Figure 811. CD1 isoheads

111

W

nk_r= 1.0‘k_2; elementaizez20x20mm
theinitialecilwaterpmeeiueheadie400m

theta: 0.6 delta=0.1eeeonde
diameter,headtube=7.1mm. diameter, cone-120mm
kcurve-l-Pachappaeanxbloam(Keet=-20mm/hr)
W-ecune-l-Zigenfuaeeoil

 

Table 312. Pressure head at each node - CD2 paremetere

3.4-mm... 6.033.... _ * ’ ‘
. head at each node (in meters): the total volume inﬂow is 3.966 E6 m’
- , III-IIIIIII-IIII-I-IIEFi ., ~ w , 77* , .

 

 

 

 

 

 

 

 

 

 

s\r 0.0 20.0 40.0 60.0 f i 8.0 10.0 12.0 14.0 16.0
(M)
0.0 1.449 1.449 1.449 1.449 2000.0 2000.0 2000.0 2000.0 4999.6
20.0 0.725 0.725 0.725 0.725 2000.0 2000.0 2000.0 2000.0 4999.1
.40.0 0.363 0.363 0.363 0.363 4991.8 4996.4 4999.1 20002 4999.1
60.0 0.242 0242 -0.168 -0465 415.2 -986.0 4997.4 20002 4999.1
-80.0 -O.362 0.453 -0599 .0880 427.7 4175.3 4996.7 20002 4999.1
400.0 -0.967 4.029 4.149 4.405 407.7 4562.7 4998.5 2000.2 4999.1
420.0 4.793 -6.878 -9.082 44.917 4876.8 4992.1 4999.5 2000.2 4999.1
440.0 49.906 49.933 49.964 49.975 4999.8 2000.0 4999.7 20002 4999.1
460.0 49.996 20.011 20.013 49.987 2000.2 20002 4999.7 20002 4999.1
480.0 49.996 20.011 20.013 49.987 2000.2 2000.2 4999.7 20002 4999.1
200.0 49.996 20.006 20.000 49.997 4999.7 2000.6 2000.0 2000.0 4999.5
00 100 120 140 160 (nm)
' 1 1 1 1 \
' /
I 9" axis
-20 I
I
I
-40 I
I
-60 :
I
-80 I
I
I
~100 |
I
—120 2 :
I
-140 _ I
I
I
-160 _ l
I
-180 - :
I
~200 2 _____________________ _I
(nn)
/ z axis

 

 

Figure 812. CD2 isoheads

 

 

 

112

W

rzk_r= 1.0‘k_z; elementsisez20x20mm
theinitialsoilwaterpree-ueheadis-mﬂm
theta:- 0.6 delta=0.1eeo¢nds

diameter, head tube a 12.7 mm, diameter, core - 80 mm

kcurve-l-Pachappaeancbloam(Ksat=20mm/hr)
W—chrve-l-Zigenfuseeoil

  

time at saturation: 10.0 minutes
head at each node (in meters):
‘

 

 
 

 

thetotalvohime inflowis2.686 E6 m’

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.0 20.0 40.0 60.0 80.0 — 100.0 120.0 140.0 160.0
1.881 1.881 1.881 49.991 49.996 40.008 49.997 40.003 49.993
0.941 0.941 0.941 49.978 40.004 40.002 49.994 40.003 49.986
0.470 0.470 0.470 49.813 49.920 49.980 49.994 40.003 49.986
0.314 -0.136 -0.418 4.012 4.992 49.633 49.996 40.003 49.986
-0.617 -0.679 -0.790 4.117 4.860 49.616 49.996 40.003 49.986
-0.969 4.068 4.237 4.963 -8.(£4 49.817 49.996 40.003 49.986
4.711 4.667 -7 .280 44.484 48.732 49.967 49.996 40.003 49.986
49.477 49.644 49.838 49.940 49.994 40.000 49.994 40.003 49.986
49.994 40.017 40.022 49.979 40.004 40.003 49.994 40.003 49.986
49.994 40.018 40.021 49.979 40.004 40.003 49.996 40.003 49.986
49.994 40.010 40.000 49.994 49.996 40.010 40.000 40.000 49.991
0 0 100 120 140 160 (mm)
' 1 1 \
I /
I 7" 0X55
-20 _ I
I
I
-40 _ I
I
-60 .. :
I
-80 _.. I
I
I
-100 _ I
|
-120 _. :
l
-140 _. I
I
I
-160 _ I
l
I
-200 _ _____________________ .J
(rm)
\/ z axis

 

Finite 313. H'I‘Dl isoheads

113

W

r:k_r- 1.0‘k_z; elementaisez20x20mm

the initial soil water presume head is 40.0 m

that.- 0.5 ““2005 seconds

diameter, headtube=3.2mm, diameter, core-80mm
kcurve-l-Pachappaeandyloam(Ksat=20mm/hr)
‘P-Scurve-l-Zigenfuseaoil

Table 314. Pressure head at each node - HTD2 parameters
, i T 7, ——
time at saturation: 11.6 minutes ‘

 

 

 

 

 

 

 

 

 

 

head at each node (in meters): the total volume inﬂow is 3.176 E6 m’
I \ r 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0
(mm)

0.0 1.814 1 .8 14 1 .8 14 49 .977 49.988 40.0% 49 .989 40.021 ~19 .991
-20.0 0.909 0.909 0.909 49.986 40.022 40.003 49.972 -20.019 -20.015
40.0 0.455 0.455 0.4 55 49.775 49.891 49.964 49.971 -20.022 -20.015
-80.0 0.0 -0.160 ~0.404 -0.894 4.598 48.327 49.976 40.021 40.015
-80.0 ~0.515 -0.587 -0.744 4.001 4.922 48.328 49.981 40.021 -20.015

400.0 -0.847 -0.933 4 . 1 28 4 .535 -3.621 48.987 49.981 40.021 40.015
420.0 -1 .449 4.710 4.808 -8.981 47.382 49.7% 49.977 20.022 40.015
440.0 47.711 48.175 48.880 49.628 49.983 49.9% 49.972 -20.021 -20.015
460.0 49.993 49.993 49.952 49.987 40.021 40.004 -19.971 40.022 -20.015
480.0 49.996 49.994 49.949 49.987 40.020 40.004 49.970 40.022 40.015
~200.0 49.995 49.997 49.985 49.974 49.999 49.9% 49.978 -20.018 49.990
_
00 100 120 140 160 (mm)
' 1 1 1 1 \
I /
I r- o xis
~20 ._ I
I
I
- 40 _ I
I
-60 _ :
I
‘80 _ I
I
I
I I
~120 _ :
I
- 140 _ I
I
I
-160 _ I
I
—180 _ :
I
-200 _._ ____________________ _I
( )
nn ¢/ z axis

 

 

Figuie 814. HTD2 isoheads

 

 

 

Appendix C

One - Dimensional VP Computer Program.

The one-dimensional VP program is given in this appendix. It was written in Pascal and
run on a Unix system.

***********************************************************************}

PROGRAM VPT(DATA,OUTFILE);
INT 1D

Unsat_lD,
FEM_1D,
Flow_1D,
Printer,

Crt,
TYPE DISiz = ARRAY[1..150] of Real;
SmSiz = ARRAY[1..5] of Real;
D3Siz = ARRAYll..50] of Integer;
NlSiz = ARRAYII..150,1..4] of Integer;
NmSiz = ARRAYll..150] of Integer;
NsSiz- — ARRAYll. .2] of Integer;
DZSiz- — ARRAYIl..150,1..27] of Real;
VAR DATA: TEXT;
OUTFILE: TEXT,
DATSAV: TEXT;
GRAPH. TEXT;
CURVE: TEXT;
SURF. TEXT;
HEAD: TEXT;
Ques: CHAR;
NAME: STRINGIZO];
TITLE: STRINGI4OI;
{DEFINITION OF THE INPUT PARAMETERS
TITLE - A DESCRIPTIVE STATEMENT OF THE PROBLEM BEING SOLVED
NP - NUMBER OF EQUATIONS (ALSO NUMBER OF NODES)
NE - NUMBER OF ELEMENTS
NCOEF - NUMBER OF SETS OF E UATION COEFFICIENTS MAXMUM OF FIVE
IPTL - INTEGER CONTROLLING T OUTPUT (4 - DEBUG OPTION)
KW - FLAG WHICH ALLOWS PRINTOUT OF DATA AT VARIOUS POINTS
IN THE PROGRAM
(KW = O OMIT WRITE, KW = l WRITE)
ahtube - area of the head tube
DHTUBE- DIAMETER OF THE HEAD TUBE
THE NUMBER OF the followm mﬁSETS of parameters MUST EQUAL NCOEF
D’I‘[I]- MATERIAL PROPER THAT MUL'I‘IPLIES THE TIME DERIVATIVE

(1= -4)
Q[I] - VEIOINUEEANTIZOEFFICIENT IN THE DIFF. EQ. NODAL COORDINATE
XClII- X COORDINATES OF THE NODES
YCII] - Y COORDINATES OF THE NODES
THE COORDINATES MUST BE IN NUMERICAL SEQUENCE RELATIVE TO N ODE NUMBERS
ELEMENT DATA
N - ELEMENT NUMBER
NMTL - INTEGER SPECIFYING THE EQUATION COEFFICIENT SET
NEL[N ,1] - NUMERICAL VALUE OF NODE I
NELlN,2] - NUMERICAL VALUE OF NODE J
DEFINITION OF THE VARIABLES READ BY INTVL2
INVL - INTEGER CONTROLLING THE INPUT OF THE INITIAL VALUES
1 - INPUTA NODE AT A TIME
2 - INPUT BY GROUPS
IBII] - NODE NUMBERS WHICH DO NOT CHANGE WITH
TIME - TERMINATE A ZERO
THETA - THE THETA VALUE USED IN THE SINGLE STEP METHOD
0 - EULER’S FORWARD DIFFERENCE METHOD
1/2 - CENTRAL DIFFERENCE METHOD
2/3 - GALERKIN’S METHOD
1 - BACKWARD DIFFERENCE METHOD
DELTA - THE TIME STEP
ITYPE - INTEGER CONTROLLING THE TYPE OF ANALYTICAL SOLUTION (1 FOR
THE VELOCITY
PERMEAMETER PROBLEM
NSTEPS - NUMBER OF TIME STEPS

114

115

IWT . INTEGER CONTROLLING THE OUTPUT OF THE CALCULATED VALUES.
VALUES ARE PRINTED
EVERY IWT
NDBC - NUMBER OF ELEMENT SIDES WITH A
DERIVATIVE BOUNDARY CONDITION
R - THE RATIO DYE/DXE, I.E., CONDUCTIVITY IN THE Y
DIRECTION IS R’CONDUCTIVITY IN THE X DIRECTION
NELIN,3] - NUMERICAL VALUE OF NODE K
NEW.“ - NUMERICAL VALUE OF NODE M
NELINA] IS SET EQUAL 'TO ZERO FOR
THE TRIANGULAR ELEMENT
num_x_ccre - number of elements in the core
in the x direction
num_y_core - number of elements in the core
in the y direction
num_x_elem - number of elemnts in the x direction
num_y elem-numberofelemntsintheydirection

Units were usa to separate procedures into compatible groups and to make reading and checking easier.

The main (controlling) program is VPT. The units are:
MyGlobal - lists the global variables and deﬁnes all variables
(other than those used only in one procedure)
FEM - ﬁnite element method, this unit includes:
Glb mtx, (creates the global stiffness and capacitance matrices)
Modify, (modiﬁes the global matrices for speciﬁed nodal

values)
Matab, (creates the A and P matrices - for Gaumian
l' . .

tron)
(decompose the A and P matricies, Gaussian method)

Mul (more Gaussian method)
Slvbd. (solution by backwards substitution)
Interact - read and write procedures:

Indata, (reads original input data from a ﬁle)

Reallts, (writes to a ﬁle)

Redata. (reads the data written by 'aeate')
Unsat_co - calculates soil/water parameters
Intw, (calculates soil moisture, based on potential)
WP C, (calculates lambda, the change in WC/change in WP)
ny. (calculates K, hydraulic conductivity)

Flow - calculates the ﬂow into the soil

Flowcalc (the only procedure)

JAR Phi, RI,
319?, Qarc: sgéSiz;

: 1:;
min., air”... NmSiz;
sat nd_ . NmSiz;
NEL: NlSiz;
c.1c mss,

NE, NSTEPS, Kapp_num,

NP. num_ﬁx_ndg

mun sat_nds, ne new, nodes: INTEGER;
Abmbe, core_radi,

dyel, Delta,

1
PROCEDURE INTVL2(np,i Lszinte er;varphi:D1siz);
~00............0........ 0......00 .00... 0.00....

S SUBROUTINE EITHER READS THE INITIAL VALUES OR
CALCULATES THE VALUES USING A PROGRAMMED UATION.
THE OPTION IS SPECIFIED BY THE INTEGER INVL WH CH
IS READ BY THE SUBROUTINE.

DEFININITION OF THE VARIABLES READ BY INTVL2
INVL - INTEGER CONTROLLING THE INPUT OF THE
INITIAL VALUES
MP1. - THE SPECIFIED VALUE FOR OPTIONS 3, 4, AND 5
1 - INPUT ONE NODE AT A TIME
INOTE - OPTION CURRENTLY UNAVAILABLE!

 

116

2 - INPUT BY GROUPS
IBEG - THE FIRST NODE IN THE GROUP
IEND - THE LAST NODE IN THE GROUP
VALUE - THE VALUE ASSIGN TO THE NODES
REPEAT UNTIL ALL NODES ARE ENTERED
3 - ZERO AT THE END POINTS AND A SPECIFIED
VALUE AT ALL THE OTHER POINTS
NOTE - OPTION CURRENTLY UNAVAILABLE!)

 

VAR
I IBEGJEND: INTEGER;
vALUE: REAL;
BEGIN
IF(IP'I'L >. 4) Then WRITELN(OUTFILE,’ ENTERING IN’I‘VL2’,’ NP - ’,NP:5);
{INPUT THE INTIAL VALUES IN GROUPS}
WRITELN(OUTI-'ILE);
readln(data,phr[1]);
REPEAT

READLN(data,IBEG, IEND,VALUE);
For I = IBEG To IEND Do

);
{OUTPUT OF THE INTIAL VALUES}

ROCEDURE INTVL2
{00000000000000.00000000000000000000000000000000000000000.00000000000000000}
Procedure INDATANu np,ne,ncoef, numb,i tl,kw:integer;

Var dhtube ,ahtube,core_ radi: reef:J

Var Dt,Q: Smsiz; Var Yc: Dlsiz;

Var Nmtlszsiz; Var Nel: leiz;

Var Ib: NmSiz; Var thta, delta: real;
ar'wvtnst? Pmt: Integer, Var Phi: Dlsiz);
{000.00.00.0000000 .0000?
{THIS SUBROUTINE AND INTVL2 READS IN THE INPUT DATA}
{{INPUT OF THE TITLE CARD AND CONTROL PARAMETERS}
VAR

I IJ,INBW,J ND,
IOLN,NB,NI‘I')E: INTEGER;
BEGIN
REAOLNmATATITLE);
writeln(title);

writeln(OUTFIIEtitl
READLNmATANP, NE, NCOEF ,I;PTL,KW)

(IN’PUT OP EQUATION COEFFICIENTS

and the NODAL COORDINATES}
REAOLNmATADHTUBE,m radi);
Amuse:- Pi'mIITUEE‘DIITUEES/«I;
WRITELN(OUTI='ILE, ' diam» head tube . ’,dhtube:0:4,
’ ooreradius’ ,core_ radi: 0:4);
FOR I . 1 'IO NCOEF DO
REAOLN(OATA.DTIII,

mmhrﬁpﬂkwmhi);

n
{GENERATION OF THE SYSTEM OF ORDINARY DIFFERENTIAL UATIONS
DEFINITION OF VARIABLES USED IN SUBROUTINE NUMO E
THTA - THE THETA VALUE USED IN THE SINGLE STEP METHODS
0 - EULER'S FORWARD DIFFERENCE METHOD
1/2 - CENTRAL DIFFERENCE METHOD
2/3 - GALERKIN’S METHOD
1 - BACKWARD DIFFERENCE METHOD
DELTA - THE TIME STEP
NSTEPS - NUMBER OF ITERATIONS
IWT - INTEGER CONTROLLING THE OUTPUT OF THE CALUCLATED
VALUES. VALUES ARE PRINTED EVERY IWT TIME STEPS.
INPUT OF THTA AND THE TIME STEP)
READLN(DAT THTA,DELTA)
WRITELN(O L,E’ Theta’ ,thta: 8. 2,’ Delta ’,delta:8:1,’ Sec
(’,DELTA/BO ::8 4,’ min)’);
{FORM THE (A) AND (P) MATRICES AND DECOMPOSE (A)
INPUT THE NUMBER OF TIME STEPS, NSTEPS,
AND THE WRTTE CONTROL, IWT}

 

U’I‘ OF THE INITIAL VALUES)

117

READLN(data,NSTEPS,IWT,prnt);
if (kw > 3) then

Writalnbutﬁle,’ NODAL COORDINATES,»
writeln(OUTFILE,’ NODE Y’);

end;
FORlzanONPDO

READ(DATA,YC[I]);
{OUTPUT OF THE EQUATION COEFFICIENTS}
{OUTPUT OF THE NODAL COORDINATES)
IHICW > 3) Then

FORI. °- 1 TO NP DO
WRITELN(OUTFILE, I. 4, ’ ’,YC[I]: 0:2),
NPUT AND ECHO PRINT O)F THE ELEMENT NODAL DATA)
WRITELN(OUTFILE, ’ ELEMENT DATA’);
WRITELN(OUTFILE, 'NEL NMTL NODE NUMBERS');

e;nd
ForKK =:1ToNEDo

BNETMIUO. - 1;

NELIKKJI :8 KK;

NELIKKJ] :s ICK + 1;

IF (KW > 3) THEN
WRITELN(OUTFILE,kkz4,NMTL[KKI:6,NEL[kk.1]:7,NEm2]:6);

End;
{INPUT THE NUMBERS OF THE NODES
WHOSE VALUES DO NOT CHANGE WTTH TIME}

{ANALYSIS OF THE NODE NUMBERS
INITIALIZATION OF A CHECK VECTOR}
{CREATION AND INITIALIZATION OF THE A
Vector CALCULATION OF THE BANDWIDTH}

writeln(outﬁle)
WRITEIJNOII’I‘FILE,’ COMPLETED READ OF ALL INPUT PARAMETERS');

WRITELN(OUTFI LE);
END; AWIBBQQEPHBEINPATN.m--- - - - - “u-

vvvvvvvvvv vvvvvvvvvvvvv vivvvvvvviv‘vvvvvivvvvvvvvvvvvvvvvvvivv'vvvvvvvvvv

IB[1] :- 1;

00000000000000.0000.000.00}

ar i: '
BEGIN
For I :s 1 to np do

dumphiﬁl: 10: 4. ’ ');
if (1 mod 6- 0) then writeln(outﬁle);

{PROCEDURE RESULTS)

Procedure Redata (Var ,ne,nbw,ncoef, numb,numb2,iptl,kw:integer;
var dhtube,r, tube,core radi: real;
Var Dt,Q :Smsis; Var Yc: Dlsiz; Var Nmtl,ib: Nmsiz;
var Nellesiz; Var ib2 :nmsiz;
Var thta, delta ,time:real; Var nsteps, iwt,prnt: integer;
PhLRf,ch:D1siz;var totﬁow, reﬁll: real);

Var
{eeeeeeeeeeeeeeeueeeaspects}

VAR
I,IJ,KK,N,NID: INTEGE -
{THIS SUERO TINE READS IN THE INPUT DATA
FROM A FILE CREATED BY A PREVIOUS RUN}
BEGIN {INPUT 0!" THE TITLE CARD AND CONTROL PARAMETERS}
READLNmATATITLE);
writeh(OUTFII.E,tiﬂ
READLNmATAN'P, NE, NCOEF ,,IPTL,KW NBW);
mm- .1) Then

begin
WOW: NP - ’.NP£.’ NE a ’.NE:6.’ NCOEF = ’,NCOEFz6);
writeln(OUTFILE,’ IPTL I ’, IPTL:6,’ KW - ’,KW:6);

{INPUT OF EQUATION COEFFICIENTS
and the NODAL COORDINATES}
READLN(DATA,DHTUBE,core radi);
READLN(DATA,THTA,DELTA)';
READLN(DATA,TIM E);
READLN(DATA,NSTEPS,IWT,prnt);
Readln(data,totflow,rdill)

e
O

118

WRITELN(0UTFILE,’ EQUATION COEFFICIENTS ');
WMngNgOI'gTFILE,’ diameter head tube ’,dhtube:0:3,’ cm; k: y -- ’,
r: : ,’ x ;
WRITELN(OUTFILE,’ SET UP Q’);
NR} :2 1 T0 NCOEF D0

LN(DATADTIII QIU)
d LNO( UTFILE,’ ’,I: 2, ' ’,DTII]:7:2,Q[I]:7:2)
93 3
AHTUBE :- Pi‘(DHTUBE‘DHTUBE)/4;
if (kw > 2) than
LN(OUTFILE,’ NODAL COORDINATES’);
.:‘ritelMOUTFILE,’ NODE X Y’);
FORI := 1 '10 NP D0
READLN(DATA,YC[I]);

{OUTPUT OF THE EQUATION COEFFICIENTS}
{OUTPUT OF THE NODAL COORDINATES)
IF(I{W>2)ThenFORI: = lTONPDO
WRITELN LN,:,(OUTFILEI4’ ’,zOYCII] '5;)
{INPUT AND ECHO PRINT OF THE ELEMENT NODAL DATA)
if (kw > 2) then

WRITELN(0UTFILE, ’ ELEMENT DATA');
WRITELN(OUTI"ILE, 'NEL NMTL NODE NUMBERS');

and;

smash.)
For Rx -- 1 To NE Do

11" (wgﬁAENNMTUIGQNI-‘ILIN.ll.NEuN.2l.NEHN.3l.NEHN.4D;
>
WRITELN(OUTFILE, N. 4 NMTLIKKI: 6 ,NELIN, 1]. 7 ,NELlN,2]£, NELIN, 3]; e ,,NEL[N 4MB);
Il-‘((N-1) <> NID) Then ’WRITELN(OUTP’ILE, 'EIJEMENTN
Nm NNOT IN SEQUENCE);

Eit'dkw>2)then

WRITEln( OUTFILE,’ COMPLETED READIN OF ELEMENTS AND
THEIR NODE NUMBERS’);
FORI- .= 1 T0 NP D0
READLNmATAJHIIID;
if (kw > 2) then
FORI := 1 T0 NP D0
writeLN(wtﬁle,PHI[I]:10:2);
{INPUT THE NUMBERS OF THE NODES
WHOSE VALUES DO NOT CHANGE WITH TIME)
readln(data,numb,numb2);
writeln(mnﬁlepumbﬂ, .2’ﬁxed nodes, and ’:,numb2 2, ’ saturated nodee’);
Fori. I 1tonumb do
READ(data,IB[I]);
For I I 1 to numb2 do
IJREAD(data,ib2[i]);
write(dv:tﬁle,’ ﬁxed nodes: ');

ICOJ);
WRITE(OUTFILE, IBIIJ]:10
Until ((IBIIJ+1] <- 0) 0R (IJ MOD 6: 0));
writeln(outf' Ile;)
writeoln(outﬁle,’ saturated nodes, including ﬁxed nodee’);
INCQJ);
WRITE(OUTFILE,IB2[IJ]:10);
Until ((IB2[IJ+1] <I 0) 0R (IJ MOD 6 = 0));
writeln(ouﬂ'
IF(NUMB > 0) eThen

NUMB > IDNN) Then

LN(OUTFILE ’ NUMBER OF BOUNDARY CONDITIONS EXCEEDS THE
ALLOWED NUMBER OF ’,idnn:6);
Exit
End;
End;
{ANALYSIS OF THE NODE NUMBERS
INITIALIZE OF A CHECK VECTOR)

119

WRITEhNtQBU'zl‘FILEé Theta ’,thta:8:2,’ DELTA ’,delta/60:8:2,’ min. C,
' : ,’ eec. ;
{FORM THE (A) AND (P) MATRICES AND DECOMPOSE (A)
INPUT THE NUMBER OF TIME STEPS, NSTEPS,
AND THE WRITE CONTROL, iwt,prnt

}

WRITELN(0UTFILE, ’number of itetatione ’,netepe:6,’ print eontrol’jwtprnt: 6);
WRITELN(OUTFILE);
if (kw > 2) then writeIxi(outﬁle,'RF);

FOR 12 1 IO NP D0

READLN(DATA,RF[I]);
if (kw > 2) then
FOR I :2 1 '10 NP D0

neadLN(data,WCP[I]);
if (kw > 2) then writeln(outﬁle,’WCP’);
if (kw > 2) then
FOR I :2 1 T0 NP D0
writeLN(outﬁI [1])
WRITELN(OUTFI TFI”?; COMPLETED READ OF ALL INPUTS ');
WRITELN(0UTFILE);

END; {PROCEDURE REDATA)
0.0“.”0”.”.O..”0000000”000.00....0.0.0..0t0000.0.00”0000000000”0”00
Procedure CREATE(np,ne,nbw, neoef ,iptl, kw, numb, numb2. integendhtube, ahtube,
come radizreal; Dt ,Q: Smaiz ;:Yc D13iz,NmtI,ib: Nmeiz;
NeI: meizdb2 .nmeiz ,thta,delta, time: real ,netepe,
wt,:rnt Intsgerfhi, rf ,wcp: dlsiz;totﬂow,reﬁll: real);
”0”.”0”. .0...“ 0......)
VAR I,J: INTEGER;
{THIS PROCEDURE WRITES THE CURRENT
DATA TO A FILE TO ALLOW RUNS WITHOUT STARTING
FROM TIME 2 0 EACH TIME.)
BEGIN {INPUT OF THE TITLE CARD AND CONTROL PARAMETERS)

.;TITLE)
WRITELN(DATSAV, NP, ’ ’ NE ’ ’ ,NCOEF, ’ ’ ,,’IPTL ’,,KW ’ ’ ,NBIV);
WRITELN(DA'ISAV, ’DH’TIIE E,“ ,eore radi);
WRITELN(DATSAV ,THTA,’ ’ ,DELTAY;
WRITELN(DATSA, ’,TIME);
WRITELN(DATSA ,’ ’,NSTEPS, ’ ’,iwt,prnt);
Wl'ite1n(dateav,’ ’ ,totﬂow,’ ’,reﬁll);
501112 1 '10 NCOEFDO
WRITELN(DATSAV, VIII],’ ’,Q[I],’ ’);
FOR I -- 1 '10 NP DO
WRITELN(DATSAV,’ ',YC[I]);
ForJ :2 1 To NE Do
WRITELN(DATSAV ,J,’ ”,NMTL[J], ',NELIJ,1],’ ’,NELlJ,2],’ ’,NEHJ,3],’
’.NEHJ.4].’ ');

FORI :2 1 T0 NP no

WRITELN(DATSAV,PHI[I],’ ');

writeln(dateavmum ,’ ’,numh2);

For I :2 1 to numb do
WRITELN(DATSAV,IB[I],’ ');

for l. 2 1 to numb2 do
wﬁteln(dateav,ib2[i],’ ');

FOR I :2 1 T0 NPDO
WRITELN(DATSAV,R1'TI].’ ');

FORI :2 1 T0 NP DO
WRITELN(DATSAV ,’WCP[I], ');

- END; - A - U {PROCEDURE CREATE)

‘ rocedure wc ca1c(kk:integer;var me. real);

«eueeeeeeﬂeeeeeee

this subroutine calculates the water content, we
at a node for a WIN,“ water potential, wp

sarrayfl. .;14]afrea1

120

wc[3] .2 0.440;

wc[13]. . 0.146;
wc[14]. - 0.112;

vim]; 2 0. 0;

wp[4] = .51. 66;

2:2; = :32:

WP ;: I340. 89;

wp[8] '2 619.8;

”III”:- -1033. .03

wp[11]. a -6166. o;

wp[12] ;. .10330. 0;

wp[13). -. -15495. 0;
(39114]:1-1032200; (0 23 )}

I) men cm

if (w‘gtcchoice?2)the

well: :2 0.372;
wc[2 :2 0.370;
wc[3: :2 0.366;
wc[4: :2 0.362;
wc[6: :2 0.363;
wc[6? :2 0.338;
wc[7: :2 0.321;
wc[8] :2 0.296;
wc[9] .2 0.;269
”[101:2 0 .229,
wc[11]. '2 0.181;
we[12]: 2 0.161;
wc[13] :2 0.137;
wc[14] :2 0.093;

and;
if(wc _choice23)then
{IenaweesoiL apherizon (0-23 cm)}

wc[1] :2 0 .;436
wc[2] :2 0 .;433
wc[3] :- 0 .;429
wc[4] :2 0.426;

 

wc[13] :- 0.199:
wc[14] :2 0.1“;

wp-the ordinate values- water content- volumetric)
-the abcissa values - water potential- in cm)
1:033 cm/bar was used)

inner-1ep < wp[l4l) then

wnteln(ouiﬁle,’ value (’,node: 12:2,

') 1e- than Ieast absiesa value: node 2 ’,zkk 0);
writeln(outﬁle,’ value (’,node: 12:2,

') Ieuthanleast absiuavalue: node 2 ’,:kk0);
writeln(outﬁle,’ the soil can not be this dryI’);

121

if (node < 40.33) then

for i_wc :2 2 to 13 do
if (node <2 Wp[i_wc]) and (node > Wp[i_wc+1]) then
(interpolate to get th - water content)
Inc :2 ((wcﬁ_wc]~wc[i wc+1])‘ln(node/wp[i_wc+1])/ln(Wp[i_wc]
nd [Wp[i_wc+1]))+ wc-I"1_wc+ 1];
e

else
no :2 wc[1];
um wc calc}

en ; graced
{”0”O”UOOO“O~O”0”O“O OM‘OOtttt$6.66”.“000000000...00.0”...)

pm Wagipfeaenvar lam_l:real;var ﬂagzinteger);

this subroutine calculates the .1 ofthe
head-water content relationship, ambda

 

we type 2 amyIl..14] ofreal;
var “an: integer;
, me: real;
wc= “.3390;
PP: m_type;
beam

(Begum ' hon'zo (0-23 )
ifbewc chug? 1 then n cm)

gm
wc[1] :2 0.463;
wc[2] :2 0.462;
wc[3] :2 0.440;
wc[4] :2 0.420;
wc[6] :2 0.377;
wc[6] :2 0.329;
wc[7] :2 0.296;
wc[8] :2 0.262;
wc[9] :2 0.236;
”[10] :2 0.207;
wc[11] :- 0.176;
we[12] :2 0.166;
wc[13] :2 0.146;
wc[14] :2 0.112;

'2: .
1, :2 0.0;
wp[2j :2 40.33;
wp[3: :2 -30.99;
wp[4: :2 «61.66;
wp[6: :2 403.3;
wp[6: :2 -206.6;
"PW: :2 -340.89;
Wp[8: :2 -619.8;
wp[9: :2 4033.0;
wp[lO] :2 -2066.0;
wp[ll] :2 -6166.0;
wp[12] :2 40330.0;
13] :2 46496.0;
14] :2 4033000;
{ pae soil, ap horizon (0-22 cm)}
if (we_choics 2 2) then

wc[1] :2 0.372;

wc[2] :2 0.370;

wc[3] :2 0.366;

wc[4] :2 0.362;

wc[6] :2 0.363;

wc[6] :2 0.338;

wc[7] :2 0.321;

wc[8] :2 0.296;

wc[9] :2 0.269;
wc[10] :2 0.229;
wc[11] :2 0.181;
we[12] :2 0.161;
wc[13] :2 0.137;
wc[14] :2 0.093;

and:

{Lenawee soil, ap horizon (0-23 cm)}
if (vuc_choice 2 3) than

wc[1] :2 0.436;
wc[2] :2 0.433;

 

122

wc[3] :2 0.429;

wc[4] :2 0.426;

wc[6] :2 0.416;

wc[6] :2 0. 4;

wc[7] :. 0.384;

wc[8] '2 0. 368;

wc[9] '2 0.333;

wc[10] :2 0.295;

wc[ll]: 2 0.247;

we[12]: . 0. 216;

wc[13]. '2 0.199;

wc[14]. '8 0.144;

and:

(wc- the u'dinate values- water content- volumetric)
{wp 2ths abcissa values - water potential- m cm}
(note: 1033 cm/bar was used)

gamma 2 0) then
“1:619! :2 philnelfkk,3]]

nods4 :2 philnell'kkAD;
“(W04 < I'p[14]) then
bean

wn'te(outﬁle,’ value has than least absiesa value: ');
w1'iteln(mrtf11e,'node42 ’,node4:0,’Wp[14] 2 ’,Wp[14]: 0: 2);
write(m1tﬁle,’thes01lcannotbeth1sd1y)
w1-itsln(ouﬂ'1le,’e17or found 1n interp’);

flag :2 1;
92d;
11' (node4 < 40.33) then
foriqlain :2 21:0 13do
if (node4 <2 wp[i_ lam]) and (11on > wp[i_la1n+1]) then

no '2 ((wcﬁ lam}wc[i lam+1])‘ln(node4/wp[i_ lam+1])l
ln(wp [i_ [i _lam+1])) + wc[i _lam+1];
«11:?- l: 2 ((mc- wc[1_la1n+1])/(node4- wp[i_la1n+1]))

and:
end
else
lamb-0;

and; mcedum interp}

AAAA AAAAAAAAAAAAAAAAAAAA-AAAAAAAAA AAAAAa
vvvvvvvvvvvvvvvvvvvvvvvvv‘vvvvvvvvvvv vvv

AAAAAAAAAAAAAAAAA-AAAAAAAAAAA
‘vv vvvvvvvvvv‘vvv'vvvvvvvvvvv

fu.«.».u§«.3§t&§°£" k1,”: real),
this procedure calculates the hydraulic conductivity

Wp type 2 arrayIl..19] of real;
var 13¢."d1’maeger;

WP: "P. ”I”:

k _an've: wp_ type;

WPII] :8 28;
Wp[2]. - .10;
"12(3) :2 -20;
W14]. 2
W6] 2 40.
wvl612-60;
WPI7] '2
w1>[8] 2 400:
wp.[9] '- 460;
“[10] :8 2200;
will] :2 ~250;
wp[12] ;a .300;
111911312400;
wp[14]. - .420;
wp[15] :- .600;
”[16] :2 2700;
wrpll'l] :2 -900;
wp[18] :2 .2000;
«@191. '2 40600;
01' 10h. sandyloam.
if (k choicewl) then

k_curve[1] :- 6.66e-4;

k_curve[2] :.'2666e-4
k_cunre[3] :.26660-4'
k_curve[4] :.2606e-4'
k_cu1've[6] :.'24O4e-4

k _.curvelG] '2 3. 64e-4

k:a11've
k_cu1ve

k_curve

h__cu1ve
k_cu1've
k_curve

k_curve-

k_curve
h_cu1ve
k_cuxve
k_curve
k_curve
h we

 

7] :'2.262e-4
8] .2.'364e-6

:9] :2 1.01e-63
[10] :2 6.30-7;
:11] :2 3.63e-7;
:12] :2 2.3e-7;
:13] :2 1.01e-7;
:14] :2 6.3e-8;
:16] :2 4040-8;
:16] :2 2.02e-8;
:17] :2 1.62e-8;
:18] :2 2.02e-9;
'19]. '2 3. 03e-10;

123

{for indie loam soil (saturated rate 20 cm/hr) :}
if (k_ choioe2 2)th

k_:curve[1] 2 6. 66e-3;
k. _curve[2] :2 5. 560-33
k_:cu1've[3] 2 6. 66e-3;
k _.curve[4] '2 6 .328e23
k_:cu1've[6] 2 6e-3;
k _:curve[6] 2 4.17e-3;
k_cu1've[7]. '2 2. 78e-3;
k _curve[8]: 2 1. 3949-3;
k _.curve[9] '2 2. 9211-63
k _curve[10] :2 1. 39e-6;
k _curve[11] :2 6. 66e-63
k _curve[12] .2 2. 929-63
1: _curv.e[13] '2 1.11e-6;
k __.cm'vell4] '2 8. 33e-7;
k _curve[16]' .2 2.39297
k_cu1've[16] .2 1 94e-.73
k _curve[17] .2 1. 39e-73
k _curve[18] :2 8. 33e-93
k cuwe[19]. 2 6. 66e-10;
.3 3
{chinoclay-saturatedvalued'Q2cm/hr}
if(kgin choice 2 3) then

be]: c11rve{1]:2 6. 6619-6;
k_cu1-ve[2]: 2 6 .3289-6
1: _:curve[3] 2 2. 78e-63
k _curve[4]. '2 2. 6e-6;
k _curve[6]: 2 1 .67e-6;
k _:curvelG] 2 1.339e—6
k _curve[7]: 2 8. 33e-63
k _.curve[8] '2 2. 78e-63
l! __:curve[9] 2 8. 33e-73
h_curve[10] .2 6. 66e-7;
k_curve[11] :2 2.929273
k_.curve[12] '2 23229-7
k _curve[13] :2 8. 33e-83
k__curve[14]. '2 6. 94e-83
k _curve[16]:2 2.92e-83
k _curve[16] :2 2. 6e-8;
h_cu1've[17] .2 1.39e-8;
k _curve[18] :2 2.377e-9
k c11rve[19]. '2 8. 330-10;

{k_curve- the ordinate valuee- water content- volumetric
taken from richards but multiplied by a constant
to Meet mah‘ic eaturated conductivity values}
{wp -the abciua values - water potential- 1n c111}

30101111,“ 2 0) then
e133.“ .2 philnelfkl 3]]

hex):- philnellkl,41];
if(head>-8)thenkl:2 h_cu1've[1]
fori_hy :2 1t018do

basin

1 24
if 2 th
13"“;- 13355?) “n

if(hsad <2 [i _kxyD and (head > wp[i_ lacy+1]) then

kl :2 curveI i ]/k_ curve[1_kxy+1])
lln(wl>ﬁ WP'IkayHD‘
i_kxy+1])) + ln(k_curve[i_1cxy+1]));

end;
klx:2 kl‘r;
if(k12 0)th81
wnte1n(outﬁle,’_en'or1,kl2 0, node’,kl: 0);
ﬂowiegovernedby theup pp:rnodee)

..‘..OOOOOOOOOOOO..‘...OOOQOOOO‘.....’..t...‘ 0.0.0.0...OOOOOOOOOOOO}

Sat elem;
..‘.‘O....O.‘DO0.0000}
VAR kk: integer;
count, epxead,
18" 01.
k1,k2: real;

count: 2 np/l;
spread. '2 philll/count;
For. kk:2 2 to up do

bean
level. '2 philkk-l]- greed;
philkk] 21221;

-AAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA.AAAAAA-AAA A AAAAAA AAAAAAAAAA-

..................................................................... ,
l’ROCEDURE INIT;((np,nbw:integer))
..‘...“..........
Var I,J: In 1';
BEGIN 3°”
For I :2 1 to up Do

”$3111.00

ForJ :21to2Do

vvvvv vvv vvvvvvvvvvvvvvvvvvvv v v 'vvv vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv’

lamiJam 1am l:rea1;11e1:nleiz;xc,yc:d1siz;
1n: dfsizmmtl: nmsiz ,:dt,q 3111312)}
{..O....‘..y.....’........’
TYPE ESiz 2 ARRAYII. ..2,1.2] ofReal;
XYSiz 2 ARRAYII..2] of Real;
VAR ECM, ESM: ESiz;
EF, Y: XYSiz;
UTE, 8°:
IQE,1e}1§:1.u
.12! ”£116,333“ INTEGER;
consi'

ES: 3812 ' ((11'1)v('111));
LMP: ESil 2 ((1,0),(0,1));
BEGIN
II" ((IPTL >2 4) AND (kk2 1)) Then WRITELN (’ENTERING GLB M It eletnf);
{BI LINEAR RECTANGULAR ELEMENT WITHOUT TH
DERIVATIVE BOUNDARY CONDITION)
{RETRIEVAL OF NODAL COORDINATES AND NODE NUMBERS}

121%: 2 yclnellkhlll- yclnelfkk.211;

IF M) <2 0.00001) Then

LN(‘THEAREAOF ELEMENT’ ,:,’kk4 IS LESS THAN 0.00001;')
WRITELNCTHE NODE numbERS ARE IN THE WRONG ORDER OR THE NODES FORM A
STRAIGHT LINE');
mWRITELNC EXECUTION TERMINATED‘);

end; (INITIALIZATION: ET[ ] MATRIX FOR THE
AXISYMMETRIC - TIME DEPENDENT PROBLEM)

125

II :2 NMTLIkk];
DTE :2 D'ITII];
Q3 =2 Q1111;

go. 2 0.0;

for 1 :2 1 to 2 do

1] :2 ‘1 ;
for 4 :2 to 2
mm - oqrtau'kz) ‘oeﬁJl/len nh;gt
I-JI I dte‘sqlﬂlaml‘lam2)‘length‘lmp[iJ]/2;
if(i.21) and 62 21) then

bean
eam[i,j]. '2 k1 ° einJ Mon 3th;
.echiJ] :2 dte‘laml‘length‘lmpﬁJW;

1m} 2) and (i - 2) then

22116.51 2 k2 ‘ «tial/1m
.;gnlig]. :2 dte‘lam23‘lenit‘tlil‘lmpﬁj1/2;
and;
end;
{Direct Stiffness}

(The major diagonal 1e stored 1n column 1,
followed by the minor diagonals (to NBW’)
incolumne2-NBW. Zeroeareueedtofill

the column to NP, the matrix size ie NP 8 NB‘V)

gar I :2 1 To 2 Do
egin
J6: NELfkkJ];
FMIJB]. '2 FM[J6] + EH1]; (force vector, has ﬁxed nodal values)
ForJ .2 1 T02

IF(JJ > 0) AND (JJ <2 2) Then

J5.JJI :2 KIJ6JJ1 + ESMUJI:
CIIJ6,JJ] :- C[J6,JJ] + ECMHJ];

End;

 

ROCEDURE GLB M

.0.....‘...‘OOOOOOOOOOOOO0.0....00“0.$6....0'.OOOOOOOOOOOOOOOOOQCOOO}

PROCEDURE MODIFY;

.0...OOOOO..OQOOOOO
{THIS SUBROUTINE MODIFIES THE GLOBAL CAPACITANCE AND
STIFFNESS MATRICES WHEN THERE ARE NODAL VALUES THAT
REMAIN CONSTANT WITH TIME.)

{W}
VAR 1, J, m, M. N: mom;
BEGIN

mum. » 4) man

LN(OUTFILE,’EXECUTING MODIFY’);
WRITELN(OUTFILE,'NP 2’,NP:6

{MODIFY C AND K MATRICES BY DELETING ROWS AND COLUMNS)
ForI :2 1T0num_ sat_ ndeDo

Begin
J. '2 sat nd _arrayII];
N. '2 J-I';
{This zeros out the row for the ﬁxed node, except
JM 2 the ﬁrd column, which in actually the diagonal.)
M :.l 14111.1;
if (M <2 NP) then

m] :2 mmymamrpmm;
KIJ Jill:-

CIJJMI=2 0

End;

(I'hisaerosoutthecohmnfortheﬁxednodmbut

rememberthedi agonalemakeup thecolumneeo

calculations muetbedone at about a46 deaweeangle}
If(N>0)th81

[N] :2 NMHNJMI‘PHIIJI;

126

[F(IPTL < 4) Thai EXIT;

IF(IPTL 2 4) Then

WRITELN(OUTFILE,’ RETURNING FROM MODIFY’)
END; PROCEDURE MODIFY}

......................................................................... ,
)'ROCEDURE MATAB;(11 tl: integer;delta,thta:1eal);)

HIS Sn BROU NE GENERATES THE A AND P MATRICES
R THE SINGLE STEP METHODS USING THE EQUATIONS
A(,) 2C[, ]+(DELTA)"I‘HTA‘K[, ]

P( ). CI l-(DELTA)‘(1-THETA)’KI I)

{neueueeeeeeeeeeeee‘eeeeeee )
VAR I,J: INTEGER;
AA. BB, TC, 'I'K: REAL;
BEGIN
II" (IPTL > 4) Then WRITELN(OUTFILE,’ ENTERING MATAB');
AA :2 THTA‘delta;
B :2 (l-TIITN‘delta;
ForI :2 1 To NP Do

Begin
ForJ.:2 1 To 2 Do

IF (IPTL < 4) Then EXIT;
WRITELN(OUTFILE,’ LEAVING MATAB’)
“END; {PROCEDURE MATAB}

AAAAAAAAAAAAAAAAAAAAA‘AAAAAAAAAAAAA -A. -AAA. .A;A.A“-A- . - - . -AAAAAA-AAAAAA-

PROCEDURE DCMPBD;§(np,nbw,iptl: .integer);}

{”O”.”O”C“O“W
VAR
itz, .I,
'11.! MK, ND,
BgG'IN NR, NL. NP'I: INTEGER;
NEW?” Then WRITELN(OUTFILE,' ENTERING DCMPBD’);
Perl :2 1 To NP1 Do

BREE-1+1;

IN)” > NP) Thai MJ :2 NP;

NJ :2 1+1;

MK: 2, {aince 11111! 2 2}
IFSNP-I-I-l) < 2) Then MK. '2 NP 1+1;

0'
1’0er NJToMJDo

:2MK-1;
ND:2 ND+1;
NLle ND+1
ForK2z21'I‘oMKDo

NK: '2 ND+K2;
Ifcl'ikgl 0 thm writeln(outﬁle,’ C[’,i: 0, ’,1] 2 ’,c[i,1]);
C[J 2 CU, K2}C[I, NL]‘C[I, NK]/C[I, 1]

End;
ImPTI. < 4) 11m EXIT;
WRITELN 1.N,'(OUTI-'ILE 'I.EAVING DCMPBD’)
ND; {PROCEDURE DCMPBD}

AAAAAAAAA"AAA‘AAAAAAAAAAAA-AAAAAAAAAA- AA‘AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA-

PROCEDURE MULTBD;{(np,nhw,iptl. 'integer; phi: d1 az'var rf:d1liz);)

127

.“O”.~O~.“O”O”O”
VAR I, J, K2, M: INTEGER;
SUM: REAL;
BEGIN
IF(IPTL >2 4) Then WRITELN(OUTFILE,’ ENTERING MULTBD');
For I :2 1 To NP Do

m :2 0.0;
K2 :2 1-1;

IF(K2>0)TI'IEN

$6M :2 SUM+KIK2JPPHIIK213
K2 :2 K2 - 1

End;
ERFII]:28UM+KII,1]‘PH1[I];

IF( < 4) Then EXIT;
WRITELN(OUTFILE,’ LEAVING MULTBD’)

 

 

-EFPLM -1-“ 11111111111111 {R BQQEPHBE 1154;911:132); 11111111111111
..................................................................... I
ZERQPPHEE§£Y§D’ 11111111111
........................... t
VAR flag,
I, J, K2
M, 11.1,
1 NP1: N1
eki , elan: INTEGER,
8U REAL;
BEGIN
NP1 :2 NP-l;
{DECOMPOSITION OF THE COLUMN VECTOR RH ))
For I :2 1 To NP1 Do
:2 1+1
Imu > NP) Then 11.1.. NP;
13.1.1 :2 1+1;
ForJz2 NJ To MJ Do

L :2 L+1;
31311;” :2 mJl-(CILLl'RFﬂl/CUJI);

{BACKWARD SUBSTITUTION NR DETERMINATION OF PHI[ 1}

PHIINPl :2 Rl'TNPl/CINP, 11;
For K2 :2 1 To NPlDo

I :2 NP-K2;
“11:2 2;
.I’F((12+1)> NP) ThenMJ. 2 NP- 1+1;
NT- 11.1-1-
SUM:- cﬁm'PRIINJ.
E5121“). :2 (RHIl-SUMVCUJI;
{the following adds a saturated
node to the fixed node eet}

k2 :2num_ eat _nde+1;
forj :2k2tonpdo

(Phiﬁl > -10.33)thu1

:20;
2:1: '2 1tonum eat ndsdo
if(eat_ 11d _arrayfll2j)thenﬂag: 2 1;
if(ﬂag2 0)then

”:l-O 0;
num_eet_nde:2num_eat nde+1;

nd_array[num_ eat _nde] :2 j;

END; {PROCEDURE SLVBD)

-AAAAAAAAAAAAAAAAAA AA‘AAAAAAAAAAA- AAAAAAAAAA-AAAAAAA‘AAAAAAAAA-AAAAAAA-

L
L

128

P

P.“O”O~&O”O“ONO}

'"gim2.,.......:"'..'2

up);
Fcrhk :2 1 To NE Do

(Ir-hp him!
if (kh- 1) then (vel :2 k1;
LAMBDMkLphimelJamlJam2ﬂag);
if exit,
GMipthhkl,kl,la1n1,lam1,nel,yc,nmﬂ,dt,qg;
End; {calcul e element matnciee
{CALCULATION OF THE ELEMENT CAPACITANCE (C3).)
{S'I'IFFNESS(F), MATRICES AND THE ELEMENT}
MODIFY( {NRCEMVEC'i‘OR- (I?!) hi)'
np,num eat eat n_a1'ray,p ,
MATAB AB( .iptl, datgthta);

DCMPB
TIME 2Tﬂl) + DELTA;
MULTanpﬁptléophiJﬂ;
FOR 1. '2 1 to npdo
RI’Iilv2 RHiI+DELTA‘I"M[il;
End,SLVBD(n111n__ eat nde,np,1ptl,rf,ph1,eat nd _array);

-AAA‘A ‘-‘;AAAAAA‘AA- A- AA

vvv v‘vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv

1
Pmcedure flow calc;

eeeeeeeeeeeeel’eeeeeeeeeeeeee}
Var wc: dleiz;

3...":- "" :.'
mapiead tube, area .mrigteger'

‘II'

B?1l|.¢01'lI’-i£m 2 0) than

1 :2 0;

_old :2 0.0;
dh rev :2 0.0'

,1 =2 Phitll;
velocityl :2 0. 0;
writeln(outﬁlllei’aifial volume of water in the permeameter 2
,’vo D

writeln(outﬁle);
'3': “mi-gaggmga'eecpm
writaiﬁllfraph, ,iterationzo, ’ ',time/GO: 0 '4 ’ ', phi[1]: 10: 3,
Exit p 2]' 10 .3 ’ ’,phi13]: 10: 3,’ ’,phi[4]: 10: 3,’ ',phiI6]: 10: 3);

End; {initialeetup}

ﬂownim:20

area core '2 12611‘(core_ radi‘core _radi);
FchKn- 1ToneDo

BEGIN

WC _CALCGLpthcIkkl); call“procedune to calculate the

elm ( water content, ”(y C[kltt] eaaléknoldﬁ}

I 'clkklo In” ‘1‘..- 2 + ,

flow-21mm Wm elm_q; core c yc
E " ’
DH:2 Hﬂownm/AH’IUBE {ahtube - area of the head tube, cm3/cm2)
if ((111 < 0) than
writeln(w warnmwdh < 0, delta q total: ’,ﬂoweum:8:2);

:2 W+ ﬂoweum;

$11113"~ "1 -

IFU‘ 'te =
B-

vvvvvvvvvvvvv'-v;:;:vvvv’

129

Writeln(m1tﬁle,’1teration',1teration:3,’ at’,time/60:6: 3, '

minutes’);
Writeln(outfile,' head at each node (in centimeters): ’);
RESU LTS(np.Phi);

End,
IF (PHI[1] < 100.0) Then

1] 2
Writeleitﬁle,’ ’Reﬁll (’,mfillz3zo,’ cm) at ’,time/60:7:2,'

.2;: . ’
1f(1terat10n>1)then

area __head tube: Pi‘Sq1(dht11be/2);
larea core: 3 P1‘16 0; {this should probably be num_ x _core':c[2]}
en :2 6
IF( teration MOD rite2 2 0) THEN
writeln(graph,’ ’ ,iteration: 0, ’ ’ ,time/GO: 0' 4, ' ’,phiﬂ]: 10: 3, ’ ’,
phi[2]: 10: 3, ’ ’ ,phi[3]: 10: 3, ’ mghiﬁ]: .:10 3,’ phi[6]:10:3);
1F(Iteration MOD rite22 0)TH
writeln(head,’ ’,iteratioui: 0, ' ’ ,time/GO:0:4,'
’,headl-head2,’ ’

,ﬂowsum);
headl :2 head2;
aid;

”0.20”.”0”.”.”.”0”0”..OOOOO0.0.0.....‘...OOOOWOOOOOOO”O”O”O
”O”.”.”.”O“O”.“.“O”.“.OOO0.0.00....‘...‘OOOOOOOOOOOOOOOOOOOO.’
BEGIN {MAIN PROG
”0”.”‘OOOOOO”O”O”O“O”0.0....COOOOOOOOOOOOOOOOOOOOMO”OOOO”O”O}
{DATA INPUT SECTION OF THE PROGRAM
- OPEN DATA FILES)
{ WRITELNC What is your data ﬁle name? (include extension”;
READLN(NAME);

)

reeet(DATA, ’46ele1n.dat’);

WRITELNC data ﬁle name: ',;name)

{WRITELNC Isthistheinitialnin YorN (Yifstarttime: 0)');
QUES:2 READKEY;

writelniges: .0);

mﬁodmLE, 'clit .tad');

rvwrite(GRAPH, ’GPH .tad'); {file to input in subsquent runs)

rewrite(HEAD,'1-1ED.tad');

WRITELN(outi'ile,’ data ﬁle name: ’am,n ame);

Write(graph,'iteration time (min) phi[l] phi[2] phi[3]’);
Writeln<snphﬁ phi“) philél');

Writeln(head,1teration time(min) dh tot_flow’);

1F(Ques= 'Y')a-(Ques2 ’y')then

For.1:2 1 to 150 Do

P ' ]:2 0.0;
rﬂi] :2 0.0;

mmw :2 0.0;

12: nd ar19ayl'i]. '02
sat n3 arrayIi]. '2 0'
11312

End;
INDATA(np,ne,ncod‘ num & nds,iptl, kw,dhtube,
ah core ra
NthN el,&_ nd array,thta, delta,
Mumpmtfhi):
0 {call procedure to read original data, cm )
time .2

"611:2 phi[l];

{phi[2] :2 -30. 0;}

voll: 2 phi[1]‘ahtube;

FIDW Calc(Phi, ycmeLtiquelta, vollmeﬁll, ahtube,dhtube,cone _radi,

ﬁel ,dhheﬁ ,velocityl, length, totﬂow,np,0 ,netepe,

iwt, rgtmemumjx ___nds,num sat __nds,Kapp num,sat_ nd _array);
011:2 '

num sat_nds:21;

sat 11d _arrayll]. '2 f: 11d _;arrq[1]

END-

ELSE

begin
Redata(np,ne,nbw,ncoJ,num_&_nds,num_sat_nds,iptl,kw,dhtube,

130

r,ahtube,c0re_ radi,Dt,Q ,Yc,Nmtl,fx_ nd _.array,Nel sat_ nd _,array,thta
delta,time,nsteps,iwt,prnt, Phi, Rf ,wcp,t0tflow, reﬁll);
iteration. '2 neteps;
ndeps:2 Mops + iteration;

(procedure to read data from previous run}

writeln(outfile);
writeln(outﬁle);
WRITELN(OUTFILE,’ NUMERICAL SOLUTION OF THE SYSTEM OF DIFFERENTIAL
EQUATIONS');
writeln(outﬁle,’ (initial values)’);
WRITELN(OUTFI
RESUL'IS(np,Phi); rocedure t0 rint initial values}
REPEAT {SOLU ON OF THE IME DEPENDENT PROBLEM)

Iteration :2 Iteration + 1;

IF(1teration MOD iwt2 0) THEN

WRITELN (’Iteration ’,iteration: 3);

if(num sat nds> num & nds) then

Sat 013m(num_ sat _ndiphi);
Fnt _aem(np,ne,n00ef,kw,iptl,num_ sat nds,delta,thta,time,dye1,
Dt,Q, YC ,_Nmtl,sat nd Parray,Ne1, Phi ,rf);

FIDWC alc(Phi,wq1,yc,nel,t1me, delta, voll, refill ,_ahtube,dhtube,core radi,
(B'e1,dh,head1 ,velocityl, length, totﬂow, ,np,iterationms,
iwt,prnt,ne,num_ fx _nds,num sat nds,Ka pp_ num,sa _array);

UNTIL (Pressed) or (Iteration- - NSTEPS)

rewrite(DA V, ’SAV. tad’);

CREATE(np,ne,nbw,n.c;d', iptl,kw,num_ fx nds,nu1n_ sat nds, dhtuhe ,ahtube,core_ radi,
DtY,Q, c,N mndtle array,Nel,sat_ nd- array,thta,de1ta,t1me

END wt, pmt,Phi,Rf,ch,totﬂow, reﬁll);

APPENDIX D

Hydraulic Conductivity and Water Content Curves

mm

The values used for the determination of water content for a given water tential are
found in table D1. This data is from the US. Department of Agriculture, 8030 Conservation
Service. Logarithmic interpolation was used to calculate the water content. Procedure
WC_CALC performed these calculations. Procedure LAMBDA used these curves to

determine the slope of the ‘1’ — 6 line.

Table D1. Water potential - volumetric water content values

 

 

 

 

 

 

 

 

-3.4089

0.295

0.321

 

-6.198

0.262

0.296

 

water WC WC WC
potential Zigenfuss Capac Lenawee
(111) soil soil soil
0.0 0.453 0.372 0.435
-01033 0.452 0.37 0.433 II
-03099 0.44 0.366 0.429 I
-0.5165 0.42 0.362 0.426
-l.033 0.377 0.353 0.416 H
-2.066 0.329 0.338 0.4 I
ll

-10.33

0.236

0.269

 

-20.66

0.207

0.229

 

-61.65

0.175

0.181

 

-103.3

0.156

0.151

 

~154.95

0.146

0.137

 

-1033.0

 

0.1_1_

 

_ _ 0.093

 

 

Hydraulic conductivity

The values used for the determination of hydraulic conductivity for a given water
ﬁotential are found in table D2. Full logarithmic interpolation was used to calculate the
ydrauhc conductivity. Procedure KXY performed these calculations.

131

132

Table D2. Water potential - k values

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

water k (mm/hr V k k
potential (Pachappa (mm/hr) (mm/hr)
(m) sandy loam) (India loam) (Chino (31‘?)
-0.08 20.0 200.0 2.00
-0.10 20.0 200.0 1.90
-0.20 20.0 200.0 1.00
-0.30 18.2 190.0 0.90
-0.40 14.5 180.0 0.6012
-0.50 12.7 150.1 0.50
-0.65 9.1 100.0 0.30
-1.00 1.3 50.0 0.10
-1.50 3.6 8.2 10.5 3.0 E-2
~2.00 1.9 E-2 0.5 2.0 E-2
-2.50 1.3 E-2 0.2 1.05 E-2
-3.00 2.3 E-3 0.1 8.0 E-3
4.00 3.6 E-3 0.04 3.0 E-3
4.20 1.9 E-3 0.03 2.5 E-3
-5.00 1.4 E-3 0.01 1.051 E-3
-7.00 7.2 E4 7.0 E-3 1.008 E-3
-9.00 5.5 E-4 5.0 E-3 5.004 E4
o20.00 7.3 E-5 3.0 E4 9.972 E-5
-105.00 1.1 8.5 2.0 E-5 3.00 E-5

 

 

 

 

APPENDIX E

Procedures for Triangular Element Analysis.

The following procedures were incorporated into the axisymmetric rectangular element
VP model. The element matricies were calculated using these procedures for the triangular
elements with the results added to the global matricies. The ﬁnite element analysis was
then exactly as described for the rectangular VP model.

ALAAAALAA ALL ALAAALAmnkk‘AALAA

‘vvw—vvav—vvvvv v VT . Tw—rvv r

procedure fnt_ elem
this procedure sets up the global matricies and
directs the ﬁnite element analysis

....__..._......}
var i_fntj_fnt,elem_fnt,ﬂag: integer;
belaing, real;
gin
for i_fnt := 1 to 200 do

begm
fm[i_ fntlt=
forj_ fnt:- = 10t060do

k?nfnt4_ fnt]. = 0. 0;
cli _fntJ_ fnt]:=
end;
end,
for beeglmem _:fnt = 1 to ne do

belcl'y(elem_ fnt, kl ,klx);

sum_ k. = sum_ k + kl

lambda(elem_ fnt ,lam_l,ﬂa
sum_lambda := sum_lambSa + lam_l;
if (nel[elem_fnt,4] = 0) then
tri_mtx(elem_fnt,kl,klx,lam_l)

else
glb_mtx(elem_fnt,kl,klx,lam_l);

end;
[calculates the element matricies
calculation of the element capacitance (c),
stiffness“), matrices and the element
force vector (0}
modify;
matab;
dcmpbd;
run_time := run_time + delta;
multbd;
fori_ fnt: = l to l“ﬂingclo
rﬂi_ fntl:-- t]+delta‘fm[i_ fnt];
slvbd;
end;

AAAAAAAAAJ AAAAaAAA A; AAJAAAAA --‘ AAAmAAAAAAAAAAA;-LJAAAAAJA

‘1 freivv -77, -e-- VVVVTTVVTT.V-CVVV,V,,,f,}

procedure tri __mtx(elem: integer; kl, klx, lam _l: real);

 

this procedure calculates the element matricies

----—----l

type esiz = arrayll..3,1..3] of real;
xysiz = arrayll..3] of real;

var ecm, esm,
es, et ,lmp: esiz;
cf, ”8:, z: xysiz;

bi,b3,bk,
ci.<a.ck.

ge, rbar: real;
i _tri _m, 11, j__tl'1_m,
11'. 1‘5. 1'6: integer;

133

134

if (iptl >- 4) and (elem s I1)) then writeln (’en 'I‘ri mtx,1u elem. ');
{BILI NEAR RECTANGULAR LEMENT WITHOUT THE
DERIVATIVE BOUNDARY CONDITIO

N}
{RETRIEVAL OF NODAL COORDINATES AND NODE NUMBERS}
fori tri_n1:- 1to3do

j tri In :-nel[elem,i_tri_n1];
3.3-:3 3.2-:3
:IT:tn_:m]:-y1:.0; ’
and;

bi :-z[2] z[3];
bi . z[3l- 2[1];

b = Ill] 2[2l;

a2=x[3] xl2];
51:- xlll - 1(3);

=x[2]-1[ll;
rbar. - (x[1] + x[2] + xI3])/3

:- 1:de 5‘((:[1]’z[2l-x[2]‘z[1]) + (XI2J‘ZI3J- zl31‘2[2])

”“31.(zl3l‘zlll-xﬂl‘zl3l)»;
1! (area <= 0)tl1 en

wuteln(mrtﬁle,’area for element ’,elemzo,’ = 0. coordinates are: ’);
fori_tri 111:: 1to3do
writeln(outfileg[i_tri_m]:10:2,’ ',z[i_tri_m]:10:2);

and;
{AXISYMMETRIC - TIME DEPEN DENT PROBLEM}
{K - 'stiﬁ'ness' matrix}

{sz}

et[1,1] := ci‘ci;
et[1,2] := ci‘ '°
st[1 ,3] .= .533;
et[2, 1] :2 st[1,2];
atl2.2]: = «3':
et[2, 3]. =

“[311] 3"- 31113];
”[312] :3 «[213];
et[3,3] := ck‘ck;
{for kDr}
es[1,1].= bi‘bi;
es[1, 2]: = bi‘b;

es[3,3] ::bk’bk;
{C- -capacitance matrix, time dependent part}
fai tri 1n: - 1to3do
for) tag-11m . 1to3do

 

 

 

 

tr1_m,i tri m] :- 0. 0;
1, T := (3 31a?) 4» all]; {for C - lumped formulation)
3:: =33":
'8 + x

{CALCULATION or ran: snmsss AND
CAPACITANCE MATRICES}

ti := nmﬂleloml;

the := dtliil;
:-= qliil;

or i_tri_m := 1 to 3 do
for j_t1‘i_1n :a 1 to 3 do

esm[i_ tri_ _tri _m]. 2:213] ‘rbar‘klx
:_trim.) (3f urea»
+ (2‘p1_‘H‘etIi_2tri tri _ml/(4‘area));
sanIi_ tri 1n.) tri n1]:= '2‘pi" te‘area’lam _I
linpﬁ_ m _mJ_ tri _;ni]/12

«my-Ln] =- pi. ‘m‘m‘dmtiml/G;

end;
{DirectStiﬁness}
{lhemaiordiagonalisetoredincolumnh
followedby the minor diagonals (to NEW)
incolumns2- NEW. Zerosaneusedtoﬁll
thecolumntoNP,tI1ematrixsizeisNPxNBW}
for i_tri_1n:- 1 to 3 do

135

bjegin
=neI[elem,1
fm[i6] “ij511 3B t1'i_m];
forj__tri_m:dB=1to3

j :8 nel[ele1n.i_t1‘i_m];
1f(jj >160) ande <2 nbw)then

]+ esm[i__ tri _mj_ tri __m];
c 6631::cﬁ65'l]+ecm[i:tri_mj_ tn _m];
m 3
end;
end;
end; {PROCEDURE 'I‘ri_Mtx}

APPENDIX F

Pressure Head for Triangular Element Model

The model with rectangular elements and triangluar elements had the following
control parameters:

 

diameter head tube = 7.1 mm

1': k_x = 1.0 * k1;

initial MC = -20.0 m

core radius: 40 mm, core length 60 mm

k curve: 1 - Pachap a sandy loam (Ksat = 20 mm/hr)
wc curve: 1 - Zigengiss soil

 

 

 

The model may be seen in Figure F1.

136

137

 

 

20 4O 60 80 100 120 (MM)

I" axis

-40

-60

-100

 

Z GXiS

 

Figure F1. Axisymmetric model with triangular elements

 

138

The soil water nessure head distribution may be seen in the output data at 15 minutes and
30 minutes (Ta les F1 and F2). Figure F2 shows the isoheads for the file with triangular

elements at 15 minutes.

the total inflow

Table F1. Pressure head at each node - triangular
model
time: 15 minutes

at 15 minutes is 2.06 E5 mm’ 3

 

 

 

 

 

z \r 0.0 20.0 40.0 60.0 80.0 100.0 120.0
(mm) (m) (m) (m) (m) (m) (m) (m)
0.0 1.402 1.402 1.402 ~20.029 49.928 49.999 -20.000
-20.0 0.701 0.701 0.701 49.999 49.933 49.999 -20.000
-40.0 0.351 0.351 0.351 49.939 49.898 49.971 49.999
60.0 6.238 6.253 6.550 4.041 6.443 49.610 49.999
60.0 6.613 6.778 4.128 -2.253 6&7 49.620 49.999
400.0 4.057 4.189 4.600 6.485 48.841 49.654 49.999
420.0 -2.547 6.613 6.577 45.970 49.890 49.999 -20.000
440.0 48.829 49.193 49.648 49.929 49.992 -20.000 -20.000
460.0 49.996 49.997 49.998 49.999 -20.000 -20.000 -20.000

 

 

 

 

 

The nodal values shown above were at the (r,z) ition indicated. This model had many
nodes that were not on this id, however, so va use for those nodes follow, with their
coorilinates as indicated. e soil water pressure heads are given in meters, coordinates in
mm .

Table F2. Pressure head at each node, triangular model, more nodes

node 24 I 6.181 at ( 40,66) node 51 = 48.794 at ( 40,66)
node 25 I 6.312 at ( 28,60) node 52 :- 48.902 at ( 48,66)
node % = 6.346 at ( 32,60) node 53 c 6.625 at ( 44,60)
node 27 x 6.374 at ( 36,60) node 54 a 6.695 at ( 48,60)
node 28 = 6.550 at ( 40,60) node 55 I 6.781 at ( 52,60)
node 29 a 6.449 at ( 36,62) node 56 a 6.630 at ( 44,62)
node 30 c 6.546 at (40,62) nods 57 :- -0.637 at ( 44,64)
node 31 = 6.404 at ( 28,64) node 58 = 6.702 at ( 48,64)
node 32 a 6.429 at ( 32,64) node 59 = 6.778 at ( 52,64)
node 33 a 6.496 at ( 36,64) node 60 x 6.651 at ( 44,66)
node 34 = 6.566 at ( 40,64) node 61 I- 6.669 at ( 44,68)
node 35 = 6.543 at ( 36,66) node 62 a 6.729 at ( 48,68)
node 36 :- 6.587 at (40,66) node 63 = 6.787 at ( 52,68)
node 37 = 6.506 at ( 28,68) node 64 8 6.886 at ( 48,-72)
node 38 = 6.522 at ( 32,68) node 65 a -20.029 at ( 60,0)

node 39 I 6.568 at( 36,68) node 66 I 49 999 at (60,-20)
node 40 a 6.610 at ( 40,68) node 67 a 49.939 at (60,40)
node 41 = 6.652 at ( 32,-72) node 68 a 4.041 at ( 60,60)
node 42 I 6.747 at ( 40,-72) node 69 c 6.897 at ( 60,68)

 

139

EU 40 bill '81! 100 1.30 14") 16.0 ' mm ')
0 LI

4/

l“ 0 1 1’5

—1L’JU

420

 

- 5:130 _,_ ____________________ _J

{ m m ‘.

 

 

Figure F2. Isoheads at 15 minutes, triangular elements

140

 

model
time: 30 minutes

Table F3. Pressure head at each node - triangular

the total inﬂow
at 30 minutes is 3.316 E5 mm3

 

 

 

 

 

 

The nodal values shown above were at the (r,z)

nodes that were not on this

6'

mm

coor)dinates as indicated. (

id, however, so va

m)
49.843 49 .997
49.854 49.995
19.593 49.802
4.342 -7.467
4.663 6.158
6.630 43.085
41.945 49.640
48.233 49.946
49.940 49.995

 

 

Table F4. Pressure head at each node, triangular model, more nodes

 

 

node 24 I 6.142
node 25 = 6.%4
node x I 6.295
node 27 I 6.320
node 28 = 6.474
node 29 = 6.386
node 30 I 6.470
node 31 I 6.344
node 32 = 6.365
node 33 = 6.426
node 34 I 6.487
node 35 I 6.466
node 36 = 6.505
node 37 I 6.426
nods 38 I 6.443
node 39 I 6.486
node 40 I 6.522
node 41 I 6.545
node 42 I 6.624

at
at
at
at
at
at
at

51:2

at
at
at
at
at
at
at
at
at

at

 

40 -72)

 

node 51 I 47.246 at ( 40,66)
is -17 284 (48 66)

nods 5 . at ,
node 53 I 6.538 n ( 44,60)
node 54 I -0 5 at ( 48,60)

%
I
9. .6. . .
N) 05
8
--:::
AAAAAA
.5 .3
g I
e§e$s$$$$s

[0°
0&0

 

 

ition indicated. This model had many
use for those nodes _follow, with their
e soil water pressure heads are given in meters, coordinates in

 

141

Figure F3 shows the isoheads after 30 minutes. The extent of each has increased, particularly in the radial
direction.

 

 

C0 80 40 60 80 100 120 140 160 mm)
I, _
J l J l J 1 l \.
’ .71 r cums
-c’0 I .I J
’ l
g ’H, _I .’q . ”|
~4U ‘1' L“ 17”» i: ,/l
. x/1
-60 m ’/

-80 _l

—100 _.

-1E>O _
-IBU a

-800 __. _____________________ ._l

 

(mm )2

Figure F3. Isoheads at 30 minutes, triangular elements

List of References

Bosch, DD. 1991, "Error associated with point observations of matric
potential in heterogeneous soil proﬁles" Trans. of the ASAE,
Vol.34(6):Nov-Dec, 1991.

Bouma, J. 1981. "Soil morphology and preferential ﬂow along macropores"
Agric. Water Manage. 3,235-250.

Bouwer, H. and RD. Jackson, 1974. "Determining soil pro erties" Drainage
for Agriculture, ed. Jan vanSchilfgaarde. Madison, , Am Soc of

Agronomy.

Burdine, N .T. 1953. "Relative permeability calculations from pore size
distribution data" Trans. Am. Inst. Min. Metal]. Pet. Eng. 198:71-87.

Chiang, C .Y., K.R. Loos and R. A. Klo p, 1991. "Field determination of
eological/chemical properties 0 an aquifer by cone penetrometry and
eadspace analysis" Ground Water, Vol.30, No.3, May-June, 1992:

428-436.

Desai, CS. 1976. "Finite element residual schemes for unconﬁned ﬂow" Int.
J. Num. Meth. Engng., 10:1415-1418

Freeze, RA. and J .A. Cherry. 1979. Groundwater. Prentice-Hall, Inc.
Englewood Cliffs, New Jersey 07632.

Gardner, W. R. 1960. "Some steady-state solutions of the unsaturated
moisture ﬂow equation wiht application to evaporation from a water
table" Soil Sci. 85: 228-232.

Gardner, W.R. 1960. "Soil water relations in arid and semi-arid conditions",
UNESCO, Arid Aone Research, Vol 15, pp. 37-61.

J abro J. D. 1992. "Estimation of saturated hydraulic conductivity of soils
from particle size distribution and bulk density data", Am. Soc. Ag. Eng,
Vol 35(2). March-April, 1992.

J abro, J .D. and DD. Fritton, 1990. "Simulation of water ﬂow from a

percolation test hole in layered soil", Soil Science Soc. Am. J.
54:1214-1218.

142

143

Kanwar, R.S., H.A. Rizvi, M. Ahmed, R. Horton and SJ Marley, 1989.
"Measurement of ﬁels—saturated hydraulic conductivity by using Guelph
and velocity permeameters" Am. Soc. of Agric. Eng. Vol.32(6):Nov-Dec,
1989.

Kinzelbach, W. 1986. Groundwater Modelling. Elsevier Science Publishers
B.V., The Netherlands.

Klute, A, 1965. "Laboratory measurement of hydraulic conductivity of
saturated soil", MW, Part 1, ed. C.A. Black, Am.
Aoc. of Agronomy, Madison, Wis., pp. 210-221.

Kunze, R.J. and DR. Nielsen; 1982. "Finite-difference solutions of the
inﬁltration equation", Soil Science, Vol. 134, N o. 2.

Kunze, R.J. and DR. Nielsen; 1983. "Com arison of soil water infiltration
groﬁles obtained experimentally and y solution of Richard’s equation",
oil Science, Vol. 135, N o. 6.

Lambe, T.W. 1951. W, John wiley & Sons, New Yorp,
165 pp.

Luthin, J .N. 1966. W, John Wiley & Sons, Inc., New

York.

McIntyre, D.S., R.B. Cunninghanm, VI Vatanakul and GA. Stewart, 1979.
"Measuring lgydraulic conductivity in clay soils: methods, techniques
and errors oil Science, Vol.128, No.3:171-183.

Maulem, Y.I. 1976. "A new model for predicting the hyaraulic condctivity of
unsaturated porous media" Water Resour. Res. 12:513-522.

Merva, G. 1979. "Falling head ermeameter for ﬁeld investigation of
hydraulic conductivity" AE Paper # 79-2515, St. Joseph, MI: ASAE

Merva, G. 1987. "The velocity permeameter technique for rapid
determination of hydraulic conductivity in-situ' Proc. Third Inter.
Workshop on Land Drainage, Ohio State Univ., Dec, 1987:G56-G62.

Neuman, SF. and PA. Witherspoon; 1971. "Analysis of nonsteady flow with
a 3free surface using the ﬁnite element method" Water Resour. Res., 7,
n , 611-623.

N euman, SR 1973. "Saturated-unsaturated seepage by finite elements"
ASCE Hydraulics Division Journal, 99, 2233-2250.

Nielsen, D.R., J .W. Biggar and KT. Ehr; 1973. "Spatial variability of ﬁeld
measured soil-water properties" Hilgardia 42:215-259.

Philip, J .P. 1951. "The theory of inﬁltration: ; 4. Sorptivity and algebraic
inﬁltration equations" Soil Sci. 84, 257-264.

 

144

Philiﬁ, J. R. 1988. "Quasianalytic and analytic ap

roaches to unsaturated
ow' 'in '

  
 

 

Raats, P.A.C. 1983. "Implications of some analytical solutions for drainage of
soil water" Ag. Water Manag., 6(1983):161-175.

Raats, P.A. C. 1988. "Quasianalytic and analytic approaches to unsaturated
ﬂ0W= commentary" W

W, 47-50.

Reynolds, W.D. and DE. Elrick, 1985. "In situ measurement of
ﬁeld-saturated hydraulic conductivity, sorptivity, and the 111-parameter

using the Guelph permeameter, Soil Science, Vol.140 No.4, October,
1985:292-302.

Richards, LA. 1931. "Capillary conduction of liquid through porous media"
Physics 1: 318-333

Robertson, P. K, R. G. Campanella, D. Gillespie and J. Green. 1986. "Use of
iezometer cone data" Proc. of In-Situ, ASCE, Specialty Conference,
lacksburg, VA.

Rogers, J. S. 1986. "Finite element model of ﬂow to an auger hole in layered
soils" Trans. of the ASAE, Vol. 29(4). July-Aug, 1986: 1005- 1011.

Rogers, J. S. and C. E. Carter, 1987. "Auger hole hydraulic conductivity
determination in layered soils" Trans. of the ASAE 30(2). 374-37 8

Rogers, J. S., H. M. Selim, C. E. Carter and J. L. Fouss, 1991. "Variability of
auger hole hydraulic conductivity values for a commerce silt loam"
Trans. of the ASAE, Vol 34(3): May-June 1991.

Rose, K.J.1988. "Use of the velocity-head ermeameter to determine the
saturated h draulic conductivity ando orseptic system site evaluation",
Thesis, Mi igan State University.

Schuh, W. M. and R. L. Cline, 1990. "Effect of soil properties on unsaturated
1514d1rg1£ic5c1113ductivity pore-interaction factors" Soil Sci. Soc. Am. J.
1

Segeil'lind, L.J., 1984. Applied Finite Element Analysis, John Wiley & Sons,
nc.

Skagl‘gls, R. W., E. J. Monke and L. F. Huggins, 1973. "Experimental
valuation of a method for determining unsaturated hydraulic
conductivity" Trans. of the ASAE, 85- 88.

Soil Conservation Service, 1986. Computed soil water retention values,
Midwest National Technical Center, Lincoln, NE.

145

Stockton, J .G. and A.W. Warrick, 1971. "Spatial variability of unsaturated
hydraulic conductivity" Soil Sci. Soc. Am. J. 35:847-848.

Taherian, M.Q., J .W. Hummel and EC. Rebuck, 1976. "Comparison of
different techniques of measuring soil hydraulic conductivity" Third
National Drainage Symposium, 109-111,122.

Toledo, P.G., R.A. Novy, H.T. Davis and LE. Scriven. 1990. "Hydraulic
conductivity of porous media at low water content", Soil Sci. Soc. Am. J.
54:673-679.

 

United Nations Environment Program, 1992. "Saving Our Planet,
Challenges and Hopes" Printed by the United Nations Environment
Program, Nairobi, Kenya, chapters 5 & 6.

Unlu, K., D.R.Nielsen, J .W.Biggar, and F.Morkoc. 1990. "Statistical
arameters characterizin the spatial variability of selected soil
ydraulic properties" Soi Sci. Soc. Am. J. 54:1537-1547.

vanBavel, C.H.M. and D. Kirkham, 1949. "Field Measurement of soil
permeability using auger holes" Soil Sci. Soc. Am. J. 13:90-96.

Zienkiewicz, 0.0., P. Mayear and Y.K. Cheung, 1966. "Solution of anisotropic
seepage problems by ﬁnite elements", J. Eng. Mech. Div. Amer. Soc.
Civil Eng, 92, EMl, 111-120.