- Email: [email protected]
An analytical solution of nonsteady evaporation from bare soils with shallow ground water table
113
A n analytical solution of n o n s t e a d y e v a p o r a t i o n f r o m bare soils with shallow g r o u n d w a t e r table
Gh. Zarei 1, M. Homaee 2 and A.M. Liaghat 3 1 & 3. Dept. of Irrigation and Reclamation Eng., University of Tehran, Karaj-4111, Iran. 2. Dept. of Soil Science, University of Tarbiat Modarres, Tehran 14155-4838, Iran.
Evaporation from bare soils plays a significant role in hydrological cycle, particularly in arid and semi-arid regions. In such areas, a large amount of the water that enters into the soil returns to the atmosphere through evaporation. In these regions, there are some areas with shallow groundwater table, evaporating a considerable amount of water into the atmosphere and accumulating salts at the upper parts of the soil profile. One major difficulty for accurate estimate of evaporation in the field conditions arises from the lack of simple function with less needed input parameters to account for water losses in soil-water balance models. The main purpose of this study was to develop and verify an analytical solution for onedimensional non-steady upward flow from shallow ground water table with minimum input data. Consequently, an analytical solution was developed based on the Richards' equation with the initial and boundary conditions governing evaporation process. In this solution, the amount and instant of evaporation from the soil surface can be estimated as function of water table drawdown, impermeable layer depth, and soil hydraulic functions. The solution is based on the parametric closed form equation of van Genuchten for retention curve. Some lysimetric experiments with packed Silty Clay, Silty Clay Loam and Sandy Loam soils were carried out to evaluate the analytical solutions. The results indicate a reasonable agreement between the collected data and the theoretical scheme. The solution seems to be promising for different soil types, requiring only few accessible input parameters. However, the small discrepancies can be attributed to side gap of shrinkaged soil, evaporation as vapor phase, experimental errors and most importantly the collapse of macropores resulting from soil packing. 1. INTRODUCTION One major component of the water budget of soil-water system is evaporation from bare soils. This may exceed more than 50 percent of the evapotranspiration (ET) in arid and semi arid regions[7]. Even, when the soil surface is covered with vegetation, depends on type of irrigation system, about 30 to 60 percent of ET belong to evaporation alone. This indicates that evaporation plays an important role in the hydrological cycle. Having shallow groundwater table, evaporation may cause severe salinization resulting in a salt crust at the soil surface and hampering agricultural activities. Physically, evaporation occurs in three recognizable stages[7]. The first, happens when the soil is wet and conductive enough to evaporate water at a rate proportional to the evaporative demand. During this stage, the evaporation rate is dominated by external meteorological
114 conditions rather than soil physical properties. In the second stage, the evaporation rate falls progressively below the atmospheric evaporativity. The evaporation rate is then limited by the rate at which the gradually drying soil profile can deliver moisture toward the evaporation zone. The third stage is established after the surface zone becomes dry such that further conduction of liquid water through it effectively stops. While many researchers investigated steady-state upward flow, evaporation usually occurs under unsteady conditions[4,11,16,17]. Steady state solutions that are based on Darcy-Buckingham's law usually can not describe the reality. Because, even where a high water table exists, neither its depth nor external conditions can remain constant for long time. Unsteady evaporation studies can be divided into numerical and analytical solutions of Richards' equation with initial and boundary conditions governing evaporation process. Analytical and quasi-analytical solution of the flow equation for semi-infinite and finite soil for an initially uniform wetted soil column in absence of groundwater table have been presented by several researchers. These studies mainly focussed on the one dimensional diffusion equation, neglecting thermal and gravity effects. Some experimental and laboratory measurements showed good agreement with these solutions[5,6,9,10,12,13]. Lomen & Warrick (1987) used general solution of a linearized moisture flow equation for point and line sources whether buried or on the soil surface. Their solution can be applied to both steady and unsteady evaporation. Brandyk & Romanowicz (1989) give an analytical solution under transient state conditions to estimate the matrix flux potential distribution within a homogeneous soil profile, using an infinite time series. In this study, an analytical solution is developed for one dimensional unsteady upward flow from water table. This solution is based on the continuity and the Richards'equation with initial and boundary conditions goveming evaporation process. It is further assumed that the water status above the water table behaves statically at the beginning of evaporation process, so that can be described with the closed form equation of Van Genuchten [14]. Also, the following assumptions are considered: (a)The soil is homogeneous and water table is shallow; (b) evaporation rate is independent of external meteorological conditions and vapor transmission of water through the soil profile is negligible; and (c) water movement in liquid form during evaporation happens under isothermal conditions. 2. THE ANALYTICAL SOLUTION In respect to Fig. 1 the following boundary (B.C) and initial (I.C) conditions are chosen:
iO(Zc ,t
=
~q(O,t) = e [q(-L,t) = 0
B.C"
IZc, (1)
:
: Zc,
I.C" [q(Zc~ '0) = e
2)
where 0 is volumetric water content (L/L), Z is vertical coordinate taken positive upward (L), Zcy is water table depth at time t (L), t is time (T), 0 s is saturated water content (L/L), q is upward flow rate ( ~ , e is evaporation rate (L/~, L is impervious layer's depth (L) and Zc~ is initial water table depth (L). Changes in water content from zl to z2 in respect to time is then:
at
_ O ~ ~O(z,t)dz
-~
~,
By using Leibnitz's rule for differentiating the integral in Eq. 3 :
(3)
115
(4)
d [ z~(,, O ( z , t ) d z - [ z~:,, ~)O(z,t) dz + O(Zcy(t ) t)dZcy(t) -O(Zcy i t = O ) dZc€ ~ d t J zc~ - J zc~ i)t ' dt ' dt
The first right hand term of Eq. 4 describes the evaporation rate (e). This can be proved by integrating the Richards' equation in respect to Z and applying the boundary and initial condition (Eqs. 1 and 2). For this purpose, we assume that there is no Sink or Source of moisture during evaporation and therefore the upward flux from the capillary fringe region to unsaturated zone equals the evaporation rate. Having the boundary condition O(Zcl ,t) = 0 s , the water content in the second term (Eq. 4) equals to 0 s .Thus, the last term is reduced to zero: A
l
a ~.$///Y.,/,./,/'.,'5, o o z
- z~,(t)
Z--O
+.~,
C~jr rZ
~0
o•
B
h:ho
m"="
O0 01"
Imp ervious
Layer
Figure 1. Schematic representation of boundary and initial conditions given in Eqs. 1 and 2.
(5)
e =m O(z t ) d z - O ~ dt zc, ' " dt
Equation 5 is a partial differential equation and can be solved analytically if the 0 - h relation and initial and boundary conditions are known. In this study, we employed the closed form equation of Van Genuchten: O(h) = Or + (0 s - 0r)[1 + (~.h)n ] -m
(6)
where 0 r is residual water content, and, c~,n and m are shape parameters. Soil water pressure head in the unsaturated zone (Fig.l)
can be described by h = Z - Z c l - h
b where hbis
bubbling pressure.The corresponding water contents then can be obtained from: O(h) = O r + (0 s - O r ){ l + [a(Z - Zcf - h b )In }-m
(7)
By substituting Eq. 7 into Eq. 5 : d zcY{Or-k-(O -- O r ) ~ + [ 6 g ( Z c f - t - h b - Z ) ] n ~ m } d z - O ~ e-dt~z~r,
dZ e ~ " dt
(8)
Assuming y = o~(Zcl + h b - Z), and integrating Eq. 8: e
=d 1 ~(Zcl-Zc~+h~) -~(--d)~h~ {O r +(0 s
_
Reusing Leibnitz's rule in Eq. 9:
0r)[l+
yn
_
]-'n}dy
0s
d Z cI dt
(9)
116
e-
Or-~-(Os-Or),.lq-__(gcf-Zcfi-+-hb).n.i
-m
-0 s
} dZ~:
dt
(10)
The integration of Eq. 10 in respect to time ( [0, t ] ) and distance([Z~:,Zc: ])leads to:
E(t) = (O, -Or)
zc:
-
_
_
(11)
The integer of Eq. 11 with X = [a'(Z~: - Z~: + hb)]n :
~
1 f [=(z.-zr
E(t)
(0,
+ ( Z cf - Z cfi) - - -
--Or)
( l + X )-m x
a n ~ ( ~ )"
( 88
(12)
With C = (crhb )", the integer of Eq. 12 can be separated into two integrals as:
.
E(t)
(0 s
-
. + ( Z~y . Zc:. )
Or )
-
1
(<,-z~: +ho)]" X (•
.
OOvl
By assuming y = X / ( I + X )
(I+X) m
dX -
x
o ----~(1+
dX
(13)
and 7"=[oC(Zc: -Z~: +hb)]"/{l+[a'(Zc: - Z ~ : +hb)]" }, the
first integral in Eq. 13 can be rewritten as:
1"]O[~
~dX(l m+X (1-1) "X)
1 = I ory"l-1(1 - y)(m-l-1) " dy
(14)
and 1/ n = p and m - 1/ n = q , with the definition of Incomplete Beta Function:
o YP-I(1- y)q-ldy - Br(p,q) <7"<1 , p > 0
(15)
and q > O(if T" = l)
Similar operations as for the first part of Eq. 13, with K = C/(1 + C) will change the second integral into: Xe
n C
K
]o(1+ dX=] o y. e
(m-l)-I
(l-y)
"
dy=BK(p,q)
(16)
substituting Eq. 15 and 16 in Eq. 13, yields:
E(t)=(O~-Or){ I [ B r ( p ' -~n q)-BK(p'q)]-(Zc:"
-Zc:)}
(17)
Thus, from Eq. 17 one can easily estimate the amount of evaporation based on the water table drawdown and soil hydraulic properties. For solving the Incomplete Beta Functions the Continued Fraction Method[I] has been chosen in this study. To determine the instant of water table drawdown, the saturated zone is considered. Concerning the thickness of this zone, Darcy's law can be written as: q = -K
A ~ h __ - K
s 5Z
- Zc: - (-Zc:
s
L-Z~:
- hb)
= -Ks
hb L-Z~
(18)
117 For this derivation, we assumed that the influx from saturated zone into the unsaturated zone equals the evaporation rate at the soil surface. Therefore, by equating Eqs. 10 and 18:
b -(O~'-Or)" L- K-' hZ----~f
{
[1 + (oC(Zcf- - Zc•+ hb)) n ]
m --1
dt
(19)
Finally, with application of Variable Separation Method and integrating from both sides of Eq. 19 in respect to time ([0,t]) and distance([Zc,,Z~i ]) we have:
K~h b t = L(Zc I _ Zc~) _ 1 2 zr (L - Zcl )dZcl 1m "2 [(zcf )2 - (zcfi ) ]--~Z;~ {l + [c~(Zcy _ Zcfi +hb)ln
(20)
(Os __ ~
J Equation 20 can be easily used to estimate the instant of evaporation when the soil hydraulic parameters, water table drawdown and elevation of impermeable layer are known. Note that integration of this equation has no analytical solution and therefore should be solved by one of the numerical or graphical methods. The numerical Romberg's method [2] is chosen for this study. 3.
MATERIALS AND METHODS
Three air-dried soils were first sieved with a 2-mm sieve and then uniformly packed at 5-cm increments in cylindrical PVC columns (with 190 cm height and 40 cm diameter). The columns were packed with Sandy Loam ( 59.25 % Sand, 32.98 % Silt and 7.77 % Clay), Silty Clay Loam (18.08 % Sand, 45.27 % Silt and 36.65 % Clay) and Silty Clay (10.55 % Sand, 47.11% Silt and 42.34 % Clay ) soils in three replicates. The first 5-cm of the columns were filled with some gravels in two layers (3 cm with coarse and 2 cm with fine gravels, respectively). The columns were connected to a charge/discharge valves and a pipe to a moveable sink. Before starting any measurements, the columns were saturated twice and allowed to drain freely. Each soil column was provided with a piezometer to read the water table depth. The soil columns were placed in the open aired site of the Research Institute, Karaj, Iran (Fig. 1). The columns were fully isolated with glass wool in order to minimize the effect of lateral thermal advection. The columns were exposed freely to the atmosphere from the top. Some physical properties of the soils are given in Table 1. Insitu 0 - h measurements at the tensiometery range (0-30 kPa)were obtained, using some tensiometers (EE514-036, ELE Inter) to measure pressure heads and volumetric method to obtain the water contents, respectively. The samples were saturated from bottom by placing them in a few centimeters of water for 72 hours. For higher soil water pressure heads (30-1500 kPa) the pressure plate apparatus was used. Since the evaporation acts as a drying phenomena, the drying retention curves were measured and used. Thus, the hysteresis did not take place during the measurements. The soil samples were packed exactly at the same bulk density as for the
118 l:'ie~ol~aeter +
r
.,,.
~" .~:, ~ - i ...:,..,, .~.
:>: SA.,
i
.~ ~" .~.. r
.~ I
8,L. I ,e ,,',~-(~,s~ I
,~.C.L
g.C,I, o ,..
,~. ~.. ~..;,
s.,.+.L
~. ,~,~.
*' ~fl %*~
>,'.<,,i .~
,~
.,
' "," , o , I
I
L
(:oll1,~i
i.i
":; I - ~ ~
"11
|=tllll|alil|ll
Figure 1. A Cross-Section of the experimental set up. Table 1. Relevant Sample
Sand
Silt
properties of the three soils. Clay
%
1 2 3
10.55 18.08 59.25
47.11 45.27 32.98
42.34 36.65 7.77
Pb
Ks
(gr/cm 3)
(cm/hr)
1.39 1.25 1.30
0.059 0.349 O.669
experimental columns. The nonlinear least square optimizer RETC [ 15 ] was used to obtain the values for the parameters in the analytical expressions of Eq. 6 (Table 2).The saturated hydraulic conductivity of the soils (K s )were determined with a constant head method in three replicates. The reported K s data here are the geometric mean of three samples which have been packed with the same bulk density as for the experimental columns. The experiment was started by raising the water table gradually until it reached the soil surface. The evaporated water was measured two times a day with an electronic scale by lifting the columns up and weighing them. Water table drawdown was determined by piezometer readings. Evaporation was take place in uncontrolled environmental conditions and thus varied greatly with external changes of temperature and relative humidity. The data were collected for the period of 64 days (from Sep. 20 to Nov. 23,2001). Table 2. The Van Genuchten parametersderived from RETC. Parameters
Texture Sandy Loam Silty Clay Loam Silty Clay
0s
(%) 52.9 49.9 58.5
~~
Or
~
n
('~~ 17.5 20.4 24.9
(cm-1) 0.0123 0.1125 0.1005
(-) 1.379 1.12 1.116
m - 1 1/n (-)
0.275 0.107 0.1004
4. RESULTS AND DISCUSSIONS This study was conducted with the primary purpose of developing and verifying an analytical model to predict evaporation from the bare soils as function of water table drawdown, impermeable layer depth, and soil physical properties. To test the derived analytical equations 17 and 20, the experimentally collected data were compared with the predicted values. Figure 3 shows the cumulative evaporation (CE) versus water table drawdown (WTD) for both measured and estimated data in different soils. This Figure indicates that the CE-WTD
119 relationship is linear up to a certain water table depth and becomes nonlinear afterward. The reason can be explained by the fact that when the soil is initially saturated and the water table is shallow, evaporation is mainly controlled by external conditions. In such a case, the capillary rise can compensate the water lost from the soil surface and thus the evaporation rate remains nearly constant. When the water table falls down, the water content in the upper part of the soil is reduced and evaporation rate is then controlled mainly by the soil watertransmitting properties. Figure 3 also shows a reasonable agreement between the measured and that estimated with Eq.17 up to a certain range of water table drawdown. These results indicate that the analytical model can predict evaporation values well for shallow water table depths, which was the prime assumption of this study. Therefore, one can determine the range of shallow water table depth for different soils, looking at the curves presented in Fig.3. In this study, the shallow water table depths were determined to be in the range of 0-70 cm for the Sandy Loam, 0-95 cm for the Silty Clay Loam and 0-140 cm for the Silty Clay soils. Furthermore, Fig 3 shows that the analytical model (Eq. 17) underestimates evaporation values when compared to measured data. This could be due to following reasons; evaporation which is take placed from the side gap established in the soil columns due to soil shrinkage, evaporation from the vapor phase in soil that is not considered by the analytical model, and the experimental errors during data collection. Figure 4 shows variation of water table depth during evaporation period for both model (Eq. 20) and measured data in different soils. This figure indicates that the variation of water table elevation during the evaporation period is not linear. This is also confirmed by the mathematical form of Eq. 20. In this equation, the time is a polynomial function of water table drawdown with more than one order. Similar to the previous results, the predicted values estimated by Eq. 20 agreed well with the measured data in the defined range of water table drawdown. When the water table depth drawdown is more than the specified range, differences between observed data and model output increase. Figure 4 also shows that Eq. 20 underestimates the water table drawdown than the measured data at a specified time. The discrepancies can be attributed to the aforementioned reasons. REFERENCES 1. M. Abromowitz and I.A. Stegun, 1965. Handbook of mathematical function. Fourth printing. Applied Math. Ser. US GovernmentPrinting Office. Washington DC. pp:944. 2. A.W. AI-Khafaji and J.R. Tooley, 1986. Numerical methods in engineering practice. CBS College Publishing. pp:390-393. 3. T. Brandyk and R. Romanowicz, 1989. Some aspects of soil moisture control for soils with shallow groundwater levels. In Proc. "Symposium groundwater management : Quantity and quality". IAHS Publ. 188:19-28. 4. W.R. Gardner, 1958. Some steady-state solution of the unsaturated moisture flow equation with application to evaporation from a water tabel. Siol Sci. 85:228-232. 5. W.R. Gardner, 1959. Solution of the flow equation for the drying of soils and other porous media. Soil Sci. Soc. Am. Proc. 23:183-187. 6. W.R. Gardner and D.I. Hillel, 1962. The relatoin of external evaporation condition to the drying of soils. J. Geophys. Res. 67:4319-4325.
120
7. 8. 9. 10. 11.
12. 13. 14. 15. 16. 17.
D.I. Hillel, 1998. "Environmental soil physics". Chapter 18:Evaporation from bare -surface soil and wind erosion. Academic Press Inc. pp:508-522. D.O. Lomen and A.W. Warrick, 1978. Linearized moisture flow with loss at the soil surface. Soil Sci. Soc. Am. J. 42:396-400. M. Menziani, S. Pugnaghi, L. Pilan, R. Santangelo and S. Vincenzi, 1999. Field experiment to study evaporation from saturated bare soil. Phys. Chem. Earth(B). 24(7):813-818. M.D. Novak, 1988. Quasi-analytical solutions of the soil water flow equation for problems of evaporation. Soil Sci. Soc. Am. J. 52:916-924. C.D. Ripple, J. Rubin and T.E.A. Van Hylcame, 1972. Estimating steady-state evaporation rates from bare soils under conditions of high water table."Geological Survey Water Supply Paper 2019-A,US Geological Survey, Washington, DC. C.W. Rose, 1968. Evaporation from bare soil under high radiation conditions. Trans. Int. Congr. Soil Sci. 9 th. Adelaide 1:57-66. G.D. Salvucci, 1997. Soil and moisture independent estimation of stage-two evaporation from potential evaporation and albedo or surface temperature. Water Resources Res.33(1): 111-122. M.Th. Van Genuchten, 1980. A Closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Siol Sci. Soc. Am. Proc. 44(5):892-897. M.Th. Van Genuchten, F.J. Leij and S.R. Yates, 1991. The RETC code for quantifying the hydraulic functions of unsaturated soils. US Environmental Protection Agency. pp:85. W.A. Warrick, 1988. Aditional solution for steady -state evaporation from a shallow water table. Soil Sci. 146:63-66. W.O. Willis, 1960. Evaporation from layered soil in the presence of a water table. Siol Sci. Soc. Am. Proc. 24:239-242.
(Sandy Loam) 250
250 ~
I ............... I
Model output
~
100
~,
1oo
0
Mode/output
9
40
so
20
o
..: 40
6O
80
lOO
120
140
160
0
180 40
water Table Depth ( c m )
60
80
lOO
12o
14o
16o
18o
..
80
'.
. . . .
1oo
120
140
160
180
200
W,,,r r = ~ o~oth (am)
water Table Depth ( c m )
Figure 3. Comparison between the experimental and modeled evaporation versus water table depth for different soils.
(Sandy Loam)
(Silty Clay)
(Silty Clay Loam) 1400 T
1600 T ......................................................................................
400
.........................................
300
12o~ I
I
.............
~
150
200
60
80
100
t
.... ~.... ;.... . .... ! 120
140
160
Water Table Depth ( c m )
180
200
"0'
....
lOO
400
40
i ]
Observation Data
2OO
~ 600
o k ~ : . - . . . .B~:
9
i
'~ ..................................
0
40
60
80
100
120
*; ....
140
Water Table Depth ( c m )
160
o,.,.:
; ....
180
2O0
80
.;...;... 100
120
140
180
Water Table Depth ( c m )
Figure 4. Comparison between the experimental and modeled water table d r a w d o w n for different soils.
180
200
A n analytical solution of n o n s t e a d y e v a p o r a t i o n f r o m bare soils with shallow g r o u n d w a t e r table
Gh. Zarei 1, M. Homaee 2 and A.M. Liaghat 3 1 & 3. Dept. of Irrigation and Reclamation Eng., University of Tehran, Karaj-4111, Iran. 2. Dept. of Soil Science, University of Tarbiat Modarres, Tehran 14155-4838, Iran.
Evaporation from bare soils plays a significant role in hydrological cycle, particularly in arid and semi-arid regions. In such areas, a large amount of the water that enters into the soil returns to the atmosphere through evaporation. In these regions, there are some areas with shallow groundwater table, evaporating a considerable amount of water into the atmosphere and accumulating salts at the upper parts of the soil profile. One major difficulty for accurate estimate of evaporation in the field conditions arises from the lack of simple function with less needed input parameters to account for water losses in soil-water balance models. The main purpose of this study was to develop and verify an analytical solution for onedimensional non-steady upward flow from shallow ground water table with minimum input data. Consequently, an analytical solution was developed based on the Richards' equation with the initial and boundary conditions governing evaporation process. In this solution, the amount and instant of evaporation from the soil surface can be estimated as function of water table drawdown, impermeable layer depth, and soil hydraulic functions. The solution is based on the parametric closed form equation of van Genuchten for retention curve. Some lysimetric experiments with packed Silty Clay, Silty Clay Loam and Sandy Loam soils were carried out to evaluate the analytical solutions. The results indicate a reasonable agreement between the collected data and the theoretical scheme. The solution seems to be promising for different soil types, requiring only few accessible input parameters. However, the small discrepancies can be attributed to side gap of shrinkaged soil, evaporation as vapor phase, experimental errors and most importantly the collapse of macropores resulting from soil packing. 1. INTRODUCTION One major component of the water budget of soil-water system is evaporation from bare soils. This may exceed more than 50 percent of the evapotranspiration (ET) in arid and semi arid regions[7]. Even, when the soil surface is covered with vegetation, depends on type of irrigation system, about 30 to 60 percent of ET belong to evaporation alone. This indicates that evaporation plays an important role in the hydrological cycle. Having shallow groundwater table, evaporation may cause severe salinization resulting in a salt crust at the soil surface and hampering agricultural activities. Physically, evaporation occurs in three recognizable stages[7]. The first, happens when the soil is wet and conductive enough to evaporate water at a rate proportional to the evaporative demand. During this stage, the evaporation rate is dominated by external meteorological
114 conditions rather than soil physical properties. In the second stage, the evaporation rate falls progressively below the atmospheric evaporativity. The evaporation rate is then limited by the rate at which the gradually drying soil profile can deliver moisture toward the evaporation zone. The third stage is established after the surface zone becomes dry such that further conduction of liquid water through it effectively stops. While many researchers investigated steady-state upward flow, evaporation usually occurs under unsteady conditions[4,11,16,17]. Steady state solutions that are based on Darcy-Buckingham's law usually can not describe the reality. Because, even where a high water table exists, neither its depth nor external conditions can remain constant for long time. Unsteady evaporation studies can be divided into numerical and analytical solutions of Richards' equation with initial and boundary conditions governing evaporation process. Analytical and quasi-analytical solution of the flow equation for semi-infinite and finite soil for an initially uniform wetted soil column in absence of groundwater table have been presented by several researchers. These studies mainly focussed on the one dimensional diffusion equation, neglecting thermal and gravity effects. Some experimental and laboratory measurements showed good agreement with these solutions[5,6,9,10,12,13]. Lomen & Warrick (1987) used general solution of a linearized moisture flow equation for point and line sources whether buried or on the soil surface. Their solution can be applied to both steady and unsteady evaporation. Brandyk & Romanowicz (1989) give an analytical solution under transient state conditions to estimate the matrix flux potential distribution within a homogeneous soil profile, using an infinite time series. In this study, an analytical solution is developed for one dimensional unsteady upward flow from water table. This solution is based on the continuity and the Richards'equation with initial and boundary conditions goveming evaporation process. It is further assumed that the water status above the water table behaves statically at the beginning of evaporation process, so that can be described with the closed form equation of Van Genuchten [14]. Also, the following assumptions are considered: (a)The soil is homogeneous and water table is shallow; (b) evaporation rate is independent of external meteorological conditions and vapor transmission of water through the soil profile is negligible; and (c) water movement in liquid form during evaporation happens under isothermal conditions. 2. THE ANALYTICAL SOLUTION In respect to Fig. 1 the following boundary (B.C) and initial (I.C) conditions are chosen:
iO(Zc ,t
=
~q(O,t) = e [q(-L,t) = 0
B.C"
IZc, (1)
:
: Zc,
I.C" [q(Zc~ '0) = e
2)
where 0 is volumetric water content (L/L), Z is vertical coordinate taken positive upward (L), Zcy is water table depth at time t (L), t is time (T), 0 s is saturated water content (L/L), q is upward flow rate ( ~ , e is evaporation rate (L/~, L is impervious layer's depth (L) and Zc~ is initial water table depth (L). Changes in water content from zl to z2 in respect to time is then:
at
_ O ~ ~O(z,t)dz
-~
~,
By using Leibnitz's rule for differentiating the integral in Eq. 3 :
(3)
115
(4)
d [ z~(,, O ( z , t ) d z - [ z~:,, ~)O(z,t) dz + O(Zcy(t ) t)dZcy(t) -O(Zcy i t = O ) dZc€ ~ d t J zc~ - J zc~ i)t ' dt ' dt
The first right hand term of Eq. 4 describes the evaporation rate (e). This can be proved by integrating the Richards' equation in respect to Z and applying the boundary and initial condition (Eqs. 1 and 2). For this purpose, we assume that there is no Sink or Source of moisture during evaporation and therefore the upward flux from the capillary fringe region to unsaturated zone equals the evaporation rate. Having the boundary condition O(Zcl ,t) = 0 s , the water content in the second term (Eq. 4) equals to 0 s .Thus, the last term is reduced to zero: A
l
a ~.$///Y.,/,./,/'.,'5, o o z
- z~,(t)
Z--O
+.~,
C~jr rZ
~0
o•
B
h:ho
m"="
O0 01"
Imp ervious
Layer
Figure 1. Schematic representation of boundary and initial conditions given in Eqs. 1 and 2.
(5)
e =m O(z t ) d z - O ~ dt zc, ' " dt
Equation 5 is a partial differential equation and can be solved analytically if the 0 - h relation and initial and boundary conditions are known. In this study, we employed the closed form equation of Van Genuchten: O(h) = Or + (0 s - 0r)[1 + (~.h)n ] -m
(6)
where 0 r is residual water content, and, c~,n and m are shape parameters. Soil water pressure head in the unsaturated zone (Fig.l)
can be described by h = Z - Z c l - h
b where hbis
bubbling pressure.The corresponding water contents then can be obtained from: O(h) = O r + (0 s - O r ){ l + [a(Z - Zcf - h b )In }-m
(7)
By substituting Eq. 7 into Eq. 5 : d zcY{Or-k-(O -- O r ) ~ + [ 6 g ( Z c f - t - h b - Z ) ] n ~ m } d z - O ~ e-dt~z~r,
dZ e ~ " dt
(8)
Assuming y = o~(Zcl + h b - Z), and integrating Eq. 8: e
=d 1 ~(Zcl-Zc~+h~) -~(--d)~h~ {O r +(0 s
_
Reusing Leibnitz's rule in Eq. 9:
0r)[l+
yn
_
]-'n}dy
0s
d Z cI dt
(9)
116
e-
Or-~-(Os-Or),.lq-__(gcf-Zcfi-+-hb).n.i
-m
-0 s
} dZ~:
dt
(10)
The integration of Eq. 10 in respect to time ( [0, t ] ) and distance([Z~:,Zc: ])leads to:
E(t) = (O, -Or)
zc:
-
_
_
(11)
The integer of Eq. 11 with X = [a'(Z~: - Z~: + hb)]n :
~
1 f [=(z.-zr
E(t)
(0,
+ ( Z cf - Z cfi) - - -
--Or)
( l + X )-m x
a n ~ ( ~ )"
( 88
(12)
With C = (crhb )", the integer of Eq. 12 can be separated into two integrals as:
.
E(t)
(0 s
-
. + ( Z~y . Zc:. )
Or )
-
1
(<,-z~: +ho)]" X (•
.
OOvl
By assuming y = X / ( I + X )
(I+X) m
dX -
x
o ----~(1+
dX
(13)
and 7"=[oC(Zc: -Z~: +hb)]"/{l+[a'(Zc: - Z ~ : +hb)]" }, the
first integral in Eq. 13 can be rewritten as:
1"]O[~
~dX(l m+X (1-1) "X)
1 = I ory"l-1(1 - y)(m-l-1) " dy
(14)
and 1/ n = p and m - 1/ n = q , with the definition of Incomplete Beta Function:
o YP-I(1- y)q-ldy - Br(p,q) <7"<1 , p > 0
(15)
and q > O(if T" = l)
Similar operations as for the first part of Eq. 13, with K = C/(1 + C) will change the second integral into: Xe
n C
K
]o(1+ dX=] o y. e
(m-l)-I
(l-y)
"
dy=BK(p,q)
(16)
substituting Eq. 15 and 16 in Eq. 13, yields:
E(t)=(O~-Or){ I [ B r ( p ' -~n q)-BK(p'q)]-(Zc:"
-Zc:)}
(17)
Thus, from Eq. 17 one can easily estimate the amount of evaporation based on the water table drawdown and soil hydraulic properties. For solving the Incomplete Beta Functions the Continued Fraction Method[I] has been chosen in this study. To determine the instant of water table drawdown, the saturated zone is considered. Concerning the thickness of this zone, Darcy's law can be written as: q = -K
A ~ h __ - K
s 5Z
- Zc: - (-Zc:
s
L-Z~:
- hb)
= -Ks
hb L-Z~
(18)
117 For this derivation, we assumed that the influx from saturated zone into the unsaturated zone equals the evaporation rate at the soil surface. Therefore, by equating Eqs. 10 and 18:
b -(O~'-Or)" L- K-' hZ----~f
{
[1 + (oC(Zcf- - Zc•+ hb)) n ]
m --1
dt
(19)
Finally, with application of Variable Separation Method and integrating from both sides of Eq. 19 in respect to time ([0,t]) and distance([Zc,,Z~i ]) we have:
K~h b t = L(Zc I _ Zc~) _ 1 2 zr (L - Zcl )dZcl 1m "2 [(zcf )2 - (zcfi ) ]--~Z;~ {l + [c~(Zcy _ Zcfi +hb)ln
(20)
(Os __ ~
J Equation 20 can be easily used to estimate the instant of evaporation when the soil hydraulic parameters, water table drawdown and elevation of impermeable layer are known. Note that integration of this equation has no analytical solution and therefore should be solved by one of the numerical or graphical methods. The numerical Romberg's method [2] is chosen for this study. 3.
MATERIALS AND METHODS
Three air-dried soils were first sieved with a 2-mm sieve and then uniformly packed at 5-cm increments in cylindrical PVC columns (with 190 cm height and 40 cm diameter). The columns were packed with Sandy Loam ( 59.25 % Sand, 32.98 % Silt and 7.77 % Clay), Silty Clay Loam (18.08 % Sand, 45.27 % Silt and 36.65 % Clay) and Silty Clay (10.55 % Sand, 47.11% Silt and 42.34 % Clay ) soils in three replicates. The first 5-cm of the columns were filled with some gravels in two layers (3 cm with coarse and 2 cm with fine gravels, respectively). The columns were connected to a charge/discharge valves and a pipe to a moveable sink. Before starting any measurements, the columns were saturated twice and allowed to drain freely. Each soil column was provided with a piezometer to read the water table depth. The soil columns were placed in the open aired site of the Research Institute, Karaj, Iran (Fig. 1). The columns were fully isolated with glass wool in order to minimize the effect of lateral thermal advection. The columns were exposed freely to the atmosphere from the top. Some physical properties of the soils are given in Table 1. Insitu 0 - h measurements at the tensiometery range (0-30 kPa)were obtained, using some tensiometers (EE514-036, ELE Inter) to measure pressure heads and volumetric method to obtain the water contents, respectively. The samples were saturated from bottom by placing them in a few centimeters of water for 72 hours. For higher soil water pressure heads (30-1500 kPa) the pressure plate apparatus was used. Since the evaporation acts as a drying phenomena, the drying retention curves were measured and used. Thus, the hysteresis did not take place during the measurements. The soil samples were packed exactly at the same bulk density as for the
118 l:'ie~ol~aeter +
r
.,,.
~" .~:, ~ - i ...:,..,, .~.
:>: SA.,
i
.~ ~" .~.. r
.~ I
8,L. I ,e ,,',~-(~,s~ I
,~.C.L
g.C,I, o ,..
,~. ~.. ~..;,
s.,.+.L
~. ,~,~.
*' ~fl %*~
>,'.<,,i .~
,~
.,
' "," , o , I
I
L
(:oll1,~i
i.i
":; I - ~ ~
"11
|=tllll|alil|ll
Figure 1. A Cross-Section of the experimental set up. Table 1. Relevant Sample
Sand
Silt
properties of the three soils. Clay
%
1 2 3
10.55 18.08 59.25
47.11 45.27 32.98
42.34 36.65 7.77
Pb
Ks
(gr/cm 3)
(cm/hr)
1.39 1.25 1.30
0.059 0.349 O.669
experimental columns. The nonlinear least square optimizer RETC [ 15 ] was used to obtain the values for the parameters in the analytical expressions of Eq. 6 (Table 2).The saturated hydraulic conductivity of the soils (K s )were determined with a constant head method in three replicates. The reported K s data here are the geometric mean of three samples which have been packed with the same bulk density as for the experimental columns. The experiment was started by raising the water table gradually until it reached the soil surface. The evaporated water was measured two times a day with an electronic scale by lifting the columns up and weighing them. Water table drawdown was determined by piezometer readings. Evaporation was take place in uncontrolled environmental conditions and thus varied greatly with external changes of temperature and relative humidity. The data were collected for the period of 64 days (from Sep. 20 to Nov. 23,2001). Table 2. The Van Genuchten parametersderived from RETC. Parameters
Texture Sandy Loam Silty Clay Loam Silty Clay
0s
(%) 52.9 49.9 58.5
~~
Or
~
n
('~~ 17.5 20.4 24.9
(cm-1) 0.0123 0.1125 0.1005
(-) 1.379 1.12 1.116
m - 1 1/n (-)
0.275 0.107 0.1004
4. RESULTS AND DISCUSSIONS This study was conducted with the primary purpose of developing and verifying an analytical model to predict evaporation from the bare soils as function of water table drawdown, impermeable layer depth, and soil physical properties. To test the derived analytical equations 17 and 20, the experimentally collected data were compared with the predicted values. Figure 3 shows the cumulative evaporation (CE) versus water table drawdown (WTD) for both measured and estimated data in different soils. This Figure indicates that the CE-WTD
119 relationship is linear up to a certain water table depth and becomes nonlinear afterward. The reason can be explained by the fact that when the soil is initially saturated and the water table is shallow, evaporation is mainly controlled by external conditions. In such a case, the capillary rise can compensate the water lost from the soil surface and thus the evaporation rate remains nearly constant. When the water table falls down, the water content in the upper part of the soil is reduced and evaporation rate is then controlled mainly by the soil watertransmitting properties. Figure 3 also shows a reasonable agreement between the measured and that estimated with Eq.17 up to a certain range of water table drawdown. These results indicate that the analytical model can predict evaporation values well for shallow water table depths, which was the prime assumption of this study. Therefore, one can determine the range of shallow water table depth for different soils, looking at the curves presented in Fig.3. In this study, the shallow water table depths were determined to be in the range of 0-70 cm for the Sandy Loam, 0-95 cm for the Silty Clay Loam and 0-140 cm for the Silty Clay soils. Furthermore, Fig 3 shows that the analytical model (Eq. 17) underestimates evaporation values when compared to measured data. This could be due to following reasons; evaporation which is take placed from the side gap established in the soil columns due to soil shrinkage, evaporation from the vapor phase in soil that is not considered by the analytical model, and the experimental errors during data collection. Figure 4 shows variation of water table depth during evaporation period for both model (Eq. 20) and measured data in different soils. This figure indicates that the variation of water table elevation during the evaporation period is not linear. This is also confirmed by the mathematical form of Eq. 20. In this equation, the time is a polynomial function of water table drawdown with more than one order. Similar to the previous results, the predicted values estimated by Eq. 20 agreed well with the measured data in the defined range of water table drawdown. When the water table depth drawdown is more than the specified range, differences between observed data and model output increase. Figure 4 also shows that Eq. 20 underestimates the water table drawdown than the measured data at a specified time. The discrepancies can be attributed to the aforementioned reasons. REFERENCES 1. M. Abromowitz and I.A. Stegun, 1965. Handbook of mathematical function. Fourth printing. Applied Math. Ser. US GovernmentPrinting Office. Washington DC. pp:944. 2. A.W. AI-Khafaji and J.R. Tooley, 1986. Numerical methods in engineering practice. CBS College Publishing. pp:390-393. 3. T. Brandyk and R. Romanowicz, 1989. Some aspects of soil moisture control for soils with shallow groundwater levels. In Proc. "Symposium groundwater management : Quantity and quality". IAHS Publ. 188:19-28. 4. W.R. Gardner, 1958. Some steady-state solution of the unsaturated moisture flow equation with application to evaporation from a water tabel. Siol Sci. 85:228-232. 5. W.R. Gardner, 1959. Solution of the flow equation for the drying of soils and other porous media. Soil Sci. Soc. Am. Proc. 23:183-187. 6. W.R. Gardner and D.I. Hillel, 1962. The relatoin of external evaporation condition to the drying of soils. J. Geophys. Res. 67:4319-4325.
120
7. 8. 9. 10. 11.
12. 13. 14. 15. 16. 17.
D.I. Hillel, 1998. "Environmental soil physics". Chapter 18:Evaporation from bare -surface soil and wind erosion. Academic Press Inc. pp:508-522. D.O. Lomen and A.W. Warrick, 1978. Linearized moisture flow with loss at the soil surface. Soil Sci. Soc. Am. J. 42:396-400. M. Menziani, S. Pugnaghi, L. Pilan, R. Santangelo and S. Vincenzi, 1999. Field experiment to study evaporation from saturated bare soil. Phys. Chem. Earth(B). 24(7):813-818. M.D. Novak, 1988. Quasi-analytical solutions of the soil water flow equation for problems of evaporation. Soil Sci. Soc. Am. J. 52:916-924. C.D. Ripple, J. Rubin and T.E.A. Van Hylcame, 1972. Estimating steady-state evaporation rates from bare soils under conditions of high water table."Geological Survey Water Supply Paper 2019-A,US Geological Survey, Washington, DC. C.W. Rose, 1968. Evaporation from bare soil under high radiation conditions. Trans. Int. Congr. Soil Sci. 9 th. Adelaide 1:57-66. G.D. Salvucci, 1997. Soil and moisture independent estimation of stage-two evaporation from potential evaporation and albedo or surface temperature. Water Resources Res.33(1): 111-122. M.Th. Van Genuchten, 1980. A Closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Siol Sci. Soc. Am. Proc. 44(5):892-897. M.Th. Van Genuchten, F.J. Leij and S.R. Yates, 1991. The RETC code for quantifying the hydraulic functions of unsaturated soils. US Environmental Protection Agency. pp:85. W.A. Warrick, 1988. Aditional solution for steady -state evaporation from a shallow water table. Soil Sci. 146:63-66. W.O. Willis, 1960. Evaporation from layered soil in the presence of a water table. Siol Sci. Soc. Am. Proc. 24:239-242.
(Sandy Loam) 250
250 ~
I ............... I
Model output
~
100
~,
1oo
0
Mode/output
9
40
so
20
o
..: 40
6O
80
lOO
120
140
160
0
180 40
water Table Depth ( c m )
60
80
lOO
12o
14o
16o
18o
..
80
'.
. . . .
1oo
120
140
160
180
200
W,,,r r = ~ o~oth (am)
water Table Depth ( c m )
Figure 3. Comparison between the experimental and modeled evaporation versus water table depth for different soils.
(Sandy Loam)
(Silty Clay)
(Silty Clay Loam) 1400 T
1600 T ......................................................................................
400
.........................................
300
12o~ I
I
.............
~
150
200
60
80
100
t
.... ~.... ;.... . .... ! 120
140
160
Water Table Depth ( c m )
180
200
"0'
....
lOO
400
40
i ]
Observation Data
2OO
~ 600
o k ~ : . - . . . .B~:
9
i
'~ ..................................
0
40
60
80
100
120
*; ....
140
Water Table Depth ( c m )
160
o,.,.:
; ....
180
2O0
80
.;...;... 100
120
140
180
Water Table Depth ( c m )
Figure 4. Comparison between the experimental and modeled water table d r a w d o w n for different soils.
180
200