Statens Planteavlsforsøgv\T

t

Landbrugsministeriet

Statens Planteavlsforsøg v\T

J 4 IAH 1993

Numerical modelling and simulation of bare soil evaporation from micro-lysimeters

Numerisk modellering og beregning a f barjordsfordampning fra mikro-lysimetre

Finn Plauborg

The Danish Institute of Plant and Soil Science Department of Agrometeorology

Research Centre Foulum DK-8830 Tjele

Tidsskrift for Planteavls Specialserie

Landbrugsministeriet

Statens Planteavlsforsøg

Numerical modelling and simulation of bare soil evaporation from micro-lysimeters

Numerisk modellering og beregning a f barjordsfordampning fra mikro-lysimetre

Finn Plauborg

The Danish Institute of Plant and Soil Science Department of Agrometeorology

Research Centre Foulum DK-8830 Tjele

Tidsskrift for Planteavls Specialserie

Beretning nr. S 2227 -1992

Preface

This report describes modelling and simulations of bare soil evaporation involving numerical approximations of the partial differential equation governing unsaturated water flow. The report finalize a study in mathematics included in my Ph.D. programme at The Royal Veterinary and Agricultural University in Copenhagen.

The Ph.D. programme is funded jointly by The Danish Research Academy and The Danish Institute of Plant and Soil Science.

I want to thank my supervisor Professor M. Flensted-Jensen, Department of Mathematics and Physics, Section of Mathematics at The Royal Veterinary and Agricultural University in Copenhagen for his instructions. Further I want to thank collegues at my department for the inspiring discussions.

Dept, of Agrometeorology December 1992

Finn Plauborg

Contents

Preface... 3

Contents... 5

Summary ... 6

Danish su m m a ry ... 6

1. Introduction... 7

2. M icro-lysim eters... 8

3. Bare soil evaporation ... 10

3.1. Physical processes 10 3.2. Modelling 12 3.3. Numerical approximations 16 4. Simulations ... 28

4.1. Duration of micro-lysimeters 28 5. D iscussion... 38

References ... 39

A. List of used symbols ... 40

Summary

The duration of micro-lysimeters for measuring bare soil evaporation is not known a priori, where duration is defined as the number of days during which measured evaporation is similar to the evaporation from the surrounding soil surface. A physically based mathematical model on bare soil evaporation is described and used to simulate evaporation from a bare silt loam soil and evaporation from a 15 cm long micro-lysimeter. The results indicate for an initially wet soil the duration of 3 days and 6 days under high and low evaporative demand, respectively.

Keywords: evaporation, micro-lysimeters, modelling, unsaturated water flow.

Danish summary

Varigheden af mikro-lysimetre til måling af fordampning fra bar jord er ikke kendt på forhånd. Varighed skal i denne sammenhæng forstås som antal dage, hvor den målte fordampning er lig fordampningen fra den omkringliggende og uforstyrrede jordoverflade. En fysisk baseret matematisk model er beskrevet og anvendt til simulering af fordampning fra en bar leijord og fordampning fra et 15 cm langt mikro-lysimeter. Resultaterne indikerer, at varigheden er 3 og 6 døgn under henholdsvis høj og lav potentiel fordampning, når jorden ved begyndelsestidspunktet er ved markkapacitet.

Nøgleord: fordampning, mikro-lysimetre, modellering, umættet vandtransport.

1. Introduction

Bare soil evaporation and the soil water content in the upper soil layers is important in relation to the germination of seeds, the activities of micro-organisms — e.g. fungi — and the timing and effect of soil cultivation.

Different methods are known for measuring bare soil evaporation. One is the use of micro- lysimeters, which are small tubes vertically installed in the soil and sealed at the bottom (Boast and Robertson, 1982). A drawback with this method is that the duration is not known a priori.

The duration is in this context defined as the number of days during which evaporation is similar to the evaporation from the surrounding soil surface.

The duration can be expected to depend on the hydraulic properties of the soil, the initial water content, the length of the micro-lysimeter and evaporative demand. Boast and Robertson (1982) compared different length of micro-lysimeters and concluded from laboratory experiments the method to be valid for 1 or 2 days. Shawcroft and Gardner (1983) stated that similar behaviour for a micro-lysimeter and the surrounding soil can be expected only if the soil water content and the water distribution are the same in and outside the micro-lysimeter.

The main objective in this study is to evaluate the duration of a 15 cm long micro-lysimeter within a drying period from simulations with a physically based mathematical model on bare soil evaporation. Also the model and the applied numerical approximations are described.

2. Micro-lysimeters

Micro-lysimeters similar to the design shown in Figure 1 were used in a field experiment at Foulum during 1991.

Tube of stainless steel

| 4 - 8,5 cm

J

15 cm

Rubber stopper

Figure 1. Outline of a micro-lysimeter.

Daily values of soil evaporation were achieved from the following procedure. The tube was installed in the soil with the upper rim levelled at the soil surface. Then the soil filled tube was removed and a rubber stopper was inserted level with the bottom rim of the tube. The weight of the micro-lysimeter was determined and the lysimeter was replaced in the hole left in the soil. Next day the micro-lysimeter was removed for weighing. An estimate of daily evaporation was then calculated from the daily weight loss of the micro-lysimeter.

For several reasons new micro-lysimeters were established every day. The most important reason was that the duration of a micro-lysimeter was not known a priori. Some of the "old"

micro-lysimeters were used in parallel with the new ones in periods of 1-15 days.

Comparisons of the evaporation data from the new and "old" micro-lysimeters could give an

r

Shawcroft and Gardner (1983) compared evaporation from micro-lysimeters and evaporation from the surrounding soil surface. The latter was calculated from a soil water balance method, where the soil-water content was measured gravimetrically from composite 1-cm-increment cores sampled daily. Another approach is to assess the duration from model simulations, cf.

Sections 3.2. and 4.

3. Bare soil evaporation

Water exchange from bare soil to the atmosphere is a physical complex phenomena overall controlled by the radiation from the sun.

3.1. Physical processes

The one-dimensional steady state energy balance for the soil surface is given by

Rn = kE + H + G W

where Rn is net radiation, A.E and H are the latent and sensible heat flux density and G is the soil heat flux density, all in Win'2. By convention Rn is positive when directed towards the soil surface and A.E, H and G are positive when directed away from the soil surface.

The latent heat flux density is equivalent to the rate of evaporation. The evaporation process can be described as the change of the state of water, where water molecules by use of energy slip from the liquid phase into the air phase.

There is a strong feedback between the fluxes in equation 1 and the soil surface conditions and properties. Net radiation is dependent on the surface albedo and the surface temperature. The albedo and the surface temperature are affected by the specific humidity at the surface. Latent and sensible heat fluxes are governed by the surface gradient in specific humidity and temperature, respectively. The turbulent exchange coefficients for latent and sensible heat are dependent on windspeed which on the other hand is affected by the roughness of the surface.

These examples also depict a mutual interrelation between the fluxes and the surface conditions. A relational diagram is shown in Figure 2.

In wet soils evaporation takes place at the soil surface. In drying soils where the level of water filled pores slowly moves downwards evaporation takes place into the air phase of the soil.

mass en ergy

s c r e e n humidity

---- 1--- g lo b a l rodioM o n

soil woter content

soil temperature

1 soil water

\

content temperature

Figure 2. Relational diagram of the soil-atmosphere system. Rectangles for main state variables, valves for rates and fluxes, circles and ellipses indicate auxiliary soil and atmospheric boundary layer variables. Solid lines represent flows of entities, broken lines flows of information. Boundary conditions are underlined. Within each column, feedback mechanisms between flux and state variable (gradient) — usually indicated by broken lines — have been omitted (Berge, 1990).

The flow processes in soil are mainly water flow in the unsaturated zone to the evaporating surface but also diffusion and convection of water vapour to the soil surface. The further transport of water vapour from the soil surface is governed by laminar and turbulent exchange processes. The soil-atmosphere pathway for water is shown in Figure 3.

Modelling this complex system to describe bare soil evaporation is not the aim of this study, but merely to adapt a well-known physically based model including unsaturated water flow and evaporation, cf. Section 3.2.

Ra=1/(C^Ua) Atmosphere

Figure 3. A schematic description of the soil-atmosphere pathway for water. q„ qg and q, are specific humidity. R(0i, and R„ are diffusion and aerodynamic resistance to the flow of water vapour (adapted from Mahfouf and Noilhan, 1991).

3.2. Modelling

The mathematical model used in this study will include only few of the physical processes involved in bare soil evaporation. More comprehensive modelling will be undertaken in the thesis within this Ph.D. study.

The model is adapted from Hansen et al. (1990). It includes a physically based description of isothermal water flow in the unsaturated soil zone (explanatory sub model) and a simple description of the evaporation process (descriptive sub model). It is believed, however, that this model can explain the major part of the water exchange especially when the soil is wet.

Therefore it is also believed that simulations can give a reasonable estimate of the duration of a micro-lysimeter, which is the main objective of this study.

In the model Darcy flow of water is assumed. The one-dimensional flow is then described by

q = - K r ø + K t ø dz

(2⁾

where z is the soil depth (cm) positive downwards, i|r is the matric potential (cm), K(ijr) is the unsaturated hydraulic conductivity (cm3 cm'2 hour') and q is the water flux density (cm3 cm'2 hour'1) positive in the positive z-direction. Flow of water satisfies the law of conservation of matter which for incompressible water and non-deformable soil is given by

30 = _ dq (3)

dt dz

where 0 is the volumetric water content (cm3 cm"3) and t is time (hours). Introduction of C(i|r) and combination of (2) and (3) leads to the governing partial differential equation (PDE) named the Richards equation

c r o * - 4 - ( K ( D p ) - m

dt dz dz dz

where

dty

is the specific water capacity (cm 1) defined as the derivative of a single valued water release function. Thus hysteresis is not included.

The PDE can be solved for known initial and boundary conditions. In this study numerical methods are applied because of the non-linearity of the PDE, cf. Section 3.3.

The initial conditions are given by

i|r = Tjr0 for 0 < z < z, and t = 0 (5) where i|f0 defines the matric potential at a well defined water content given by the iJr-0 relationship and z, is the soil depth at the lower boundary.

Only specific boundary conditions will be considered in this study:

Upper boundary

q = - E„ for z = 0 and t > 0 Lower boundary

(

)

q = qg for z = z, and t > 0 or

(

7

q = 0 for z = zm and t > 0 (8)

where E„ is the actual evaporation (cm3 cm'2 hour'1) positive in the negative z-direction, qg is the gravity flux density (cm3 cm'2 hour'1) and zm is the depth of the micro-lysimeter.

Assuming the flow at the lower boundary to be governed by only the gradient in gravitational potential is appropriate when the ground water table is far away from the modelled depth- domain (Jensen, 1983).

Actual evaporation depends on the hydraulic conditions in the soil and the evaporative demand of the atmosphere and is not known a priori. A solution is then given by

is the maximum possible flux density. Ep (cm hour'1) is the potential evaporative demand of the atmosphere calculated from equation 10 (Hansen et al., 1990).

where is the daily potential evaporation (cm day') and t is the hour of the day, t = 1 ,2 , ... 24. The parameter E^ is in this study assumed constant, E^ e [1,5].

E _P E i , 0_P m in tk lv E ] £ > 0

(

9

where

(10)

The maximum possible flux density q, is calculated in accordance with Hansen et al. (1990) by assuming the water content at a minimum level 0O at z = 0.

The hydraulic parameter function K(i|r) is developed using the model of van Genuchten (1980).

The water release function 0(h) is in this model given by

6(h) = (0 - 0j[l + (a h ) f + 0 O 1) where h is the tension (cm) (= -ijr), 0, and 0f are the saturated and residual (e.g., at wilting point) water contents (cm3 cm'3), respectively. Assuming that m = 1 - n'1 and 0, and 0, are known, then n and a (cm 1) can be found by non-linear regression applied to measured values of (h,0(h)).

The inverse function h(0) is found from equation 11 as

h(6) = — a

0 - 0 0 - 0 The water retention function is then given by

t|<0) = - h(6)

The unsaturated hydraulic conductivity is calculated from equation 14.

. [l - (a h)"~' • [l + (a h fX m = K }

[l + (a h y

(12)

(13)

(14)

where K, is the saturated hydraulic conductivity (cm3 cm'2 hour'1). The mathematical derivation of equation 14 is given by van Genuchten (1980).

The capacity function is found by differentiation of equation 11 with respect to h.

C(h) = - m (0s - 0r) [l + (a h)"]""' n (a h y ' a The specific water capacity expressed as a function of i|r is found by

C(i|r) = - C(h) (16) 3.3. Numerical approximations

The Richards equation (4) can be solved analytically when special initial and boundary conditions are known or in fact only "quasi-analytical", because some coefficients in the solutions are obtained by numerical methods.

Solution of equation 4 for known initial and boundary conditions can in general only be obtained by using numerical methods.

A widely used method is the finite difference scheme where derivatives are approximated with their difference form.

The difference notation j, i, Azj and Ati+I can be visualized as a grid superimposed on the depth-time domain as shown in Figure 4.

In this notation used by Jensen (1983) and Hansen et al. (1990) i|r- represent iJr(Zj,t1) or, from the retention curve, the water content 6(Zj,t1) in the depth interval or layer Azj to time t‘.

Soil surface z0 j = l z, j=2 z2 j=3 z3

j-1 Zj.

j Z|

j+1 _{‘ i+ l}

j=N zN

i=0 i=l i=2 t° t 1 t2

-LZy ■

-AZ:

-Azm

t'+1 i Upper boundary --- ►Time

<— At1'

♦L,

*i

*i+. — I *i:i

T Lower boundary

The derivative with respect to time in (4) at depth zs is approximated by a forward difference (17) at a t M

and the space derivative in (3) can be approximated by a centered difference at depth zi to time ti+1,

_ dq = _ C - g “

dz A Zj

Or the approximation can be centered in both space and time.

(18)

^ - - 'A

dz A ^{Z ,}

(19)

Using centered values of K(ijr), the flow variables on the right sides of (18) and (19) are given by

- ü rM V - 1 7J - z Mj- 1

K

z j<

+ k M

+

* (c)

(20)

(21)

(22)

(23)

z j - z j->

With C(i|r) centered in time substitution of equation 17, 18, 20 and 21 into (4) leads to the

"fully" implicit finite difference approximation where evaluation of the space derivatives are put forward in time, cf. equation 24.

(24)

With C(tJt) centered in time substitution of equation 17, 19, 20, 21, 22 and 23 into (4) leads to the Crank-Nicolson finite difference approximation where the space derivatives are evaluated at time ti+'\

The "fully" implicit formulation 24 is first and second-order accurate in time and space, respectively and thereby less accurate than the Crank-Nicolson formulation 25 which is second-order accurate in both time and space (Chapra and Canale, 1988). Both methods can be characterized by being unconditional stable and convergent for A t-0 and A z-0 when applied to linear parabolic differential equations.

In the case of the non-linear Richards equation the approximations lead to non-linear finite difference equations, which are not straight forward to solve.

Equations 24 and 25 can be rearranged to explicitly express the unknown dependent variables (25)

expressed as linear algebraic equations valid for each nodal point, cf. (26) respectively,

where

öür'*1 + = r

aj =

m 1;*) (Az)2

(A z)2

(A z)2

or

where

A u r 1 + _{j * j-i} B l ir 1 + _j D Ur'*1 = R _j

’ V i j

A - - K ( O

' (A z)2

B, _ 2 c ( t r ) ^ k f e j + k (* ;:)

A t M (A z)2

D - - * M (A z)2

and (31)

(26)

(27)

(28)

(29)

(30)

(

31

(32)

(33)

(34)

R,

(Az)2

K f e )

(A z f

m-i

_ *(<d „,, _ * M - k (»;:)

( i z ) ‘

*(«d - « ta ll

T;*i

(35)

As mentioned the finite difference equations are difficult to solve because of the non-linearity.

The hydraulic parameter functions C(ijr) and K(i|r) are assumed constant within a time interval and depend on the unknown at time ti+* at one or two depth levels, respectively, cf. equation 26 and 31. Different approximations have been proposed to overcome these problems.

Feddes et al. (1978) approximated i|ri+w by extrapolations of some functional relationship involving »Jr11, i|r! and time step size At and Ati+I, cf. equations 36 and 37.

i + A L L 2 A t ' A t‘*'

(< • + *') - (36)

°-25 fe + *'")

= 0.5 1 + A r 1 0.25

2 A t ' A t w

(*;+ <■) -

A t ‘

(♦r *

From these approximations, the unsaturated hydraulic conductivity can be found from the functional relationship K(i|r), where K(-i|r) = K(h) cf. (14).

To assess the capacity C(i|rj+W) Feddes et al. (1978) applied the following approximation

= 0.5 (*;;* + * ; ; j (38)

With these approximations Feddes et al. used the Crank-Nicolson method to solve the Richards equation. To overcome convergence and stability problems Feddes et al. restricted the time step size to the following conditions given by Zaradny (1978).

r

Ar'*1 < 6 A z (39)

k r

where q is the actual flux at the lower or upper boundary at time t1 and e is a factor where 0.015 < e < 0.035. For high value of q the lower value should be assigned to e and vice versa.

Jensen (1983) used an iterative algorithm to calculate the hydraulic properties at time ti+l\

which he found could eliminate stability and convergence problems when used in combination with the "fully" implicit or the Crank-Nicolson method. This approach will be described in the following and be implemented in the simulation model. Nodal points involved in calculation of the space averages K(iJrj±1J are shown in Figure 5.

♦ j - i t

j-l --- •--- 1--- Az ---® ---

t j A

J --- •--- 1-- Az

K(*j+Jj) ▼

---® ---

*3+1 t

3+1 --- •--- 1-- Az

▼

Figure 5. Nodal points ( • ) involved for calculation of K(i|rj±Vi) (®).

Variations in soil properties and water content among calculation layers Azj mean that estimation of K(i|rj±,J is not simple especially because the K(i|r) relation is highly non-linear.

Often used weighting methods are the arithmetic and the geometric means cf. equation 40 and 41, respectively.

y _ (K M (40)

K K . ) - /(* W K (♦ J

Jensen (1983) tested numerical solutions using both weighting methods and found that the geometric mean in general was superior to the arithmetic mean. Equation 41 will therefore be implemented in the simulation model.

Jensen (1983) found that an important criterion for achieving stability in the numerical solutions was how the time averaging was performed when calculating K(t|ri+W) and C(i|r,+''4).

He found that the following procedure had a stabilizing effect without introducing problems with convergence.

c (*;♦*) = 0.5 c ( r ; (A*)) + 1 £ c (*;*’ (m)) (42)

k te* M + i t K te« (»»))

where M is the number of iterations within the current time step.

The iterative procedure involved for calculation of is terminated when the following convergence criteria are fulfilled for all j (Hansen et al., 1990),

nr;-' (M) - nr1(M-1) r ; (M -1) ^{< 6.}

(44)

or

(AO - (M-l)| < 62 (45) where t|ri+l (M) represents the solution to the iterative step M. The time step is decreased by a factor two between the iterative steps. If convergence is not achieved at a certain number of iterations, the calculations continue without convergence. The iteration procedure is

Substitution of the approximations (42) and (43) into equation 26 and applying this equation at all node points in the depth-time domain, a solution can be found by solving the following

matrix equation for known boundary and initial conditions,

L - 1 = f (46)

which describes a system of N linear equations with N unknowns represented by the vector ifc the vector r and tri-diagonal matrix L, cf. Figure 6.

<^N-1

b„ 1C1

Figure 6. Illustration of the matrix equation 46.

A simular matrix equation can be formulated to solve (31).

M ■ ip = R (47)

The equations at the first and last node points j = l and j= N , respectively, depend on the boundary conditions. If the boundary is known as a specified potential, no equation has to be solved at this node and the total number of equations to be solved is reduced by one. If the boundary is known as a flux, the equation to be solved can be deduced by involving an imaginary node point outside the depth-time domain, cf. Figure 7.

The known fluxes at the upper and lower boundary qj^j and qj^*, respectively can be approximated by

i - - X * « ( * ; . : ) (4 8 )

* K ( ♦ - ) (4 9 )

Solving for TjC'in (48) and in (49) gives A z

♦o

---4

f 2

a.

Az Az

— upper boundary ---

Az Az

lower boundary —

Az Az

b .

♦m-i

♦n

Figure 7. a. Imaginary node point i|r0 at the upper boundary, b. Imaginary node point i(rN+1 at the lower boundary.

+ (50)

tfrM = ürM -Y AM T N (51)

Substitution of (50) into equation 26 for the note j = 1 leads after rearrangement to the first equation in the matrix equation 46 which is valid for the "fully" implicit method.

where

b i|iM + I * 1 d ilr'*1 = 1 *2 rI (52)

(53)

d, = K

M

(Az)2

(54)

_ £ k l, „ K M -K (♦::)

At M 1 A z

K (< j

(A z)2

k : -g (* d ]Az

K ( K l )

(55)

simplified to

A f w 1 A z

Substitution of (51) into equation 26 for the node j= N , leads after rearrangement to the last equation in the matrix equation 46.

a ilr'*1 + b tir'*‘ = r (57)

“ W ^ y-1 "n V N N

where

a = _ K (* C ) (58)

(A z)2

b m C M + (59)

* A t M (A z)2

A / M " A z

_

4 c ) [ c

< iz)’

simplified to

(61)

For obtaining the first and last equation in the matrix equation 47 valid for the Crank-Nicolson method, a simular procedure should be used. Here the known fluxes at the boundaries should be expressed at time ti+l and t1 (cf. equations 20-23) and solved for ijrj, and i|rj,+l, and i|rj<+l and respectively. Substitution of these expressions into equation 31 forj = l and j = N would after rearrangement lead to the equations needed.

These derivations are not shown because only the "fully" implicit method will be implemented in the simulation model. The "fully" implicit method was also used by Jensen (1983) who found that less iterations were needed with this method compared to the Crank-Nicolson method.

In this study the upper boundary flux is equal to the actual evaporation, cf. Section 3.2.

equations 6, 9 and 10. Then in equation 56

The lower boundary flux is set to zero or to gravity flow, cf. Section 3.2, equations 7 and 8.

Then in equation 61

E M E ws 0

P P (62)

min[|<?.M1,2?M] £ w > 0 where qi+l is approximated by

(63)

o " 1 = 0*N*'A (64)

or

(65)

F

(66⁾

4. Simulations

The "fully" implicit numerical model described in Section 3.3 has been programmed in Turbo Pascal. Simulations of water flow and actual evaporation can then be performed for defined systems of interest.

In this study simulations are performed for a silt loam soil which is assumed homogeneous and characterized by the following parameters:

K, = 1.2 cm hour'1 0O — 0.061 cm3 cm'3 8, = 0.48 cm3 cm'3 0, = 0.061 cm3 cm 3 a = 2.452-10'2 cm'1 n = 1.568

m = i - n-' = 3.622-10-'

The van Genuchten parameters are calculated from tabulated retention data listed by Hanks (1991).

The soil is initially assumed at field capacity equivalent to a water content of 30 vol %. The only driving input variable is (cm day ') which is assumed constant.

4.1. Duration of micro-lysimeters

To describe and evaluate the duration of a micro-lysimeter simulations are performed for two different systems.

Actual evaporation is calculated from a 100 cm soil profile. In this system, free drainage is assumed from the profile and approximated by gravity flow, cf. Section 3.2.

In the other system, actual evaporation is calculated from a 15 cm long micro-lysimeter. At the bottom no water exchange is possible because of the stopper. This is introduced in the simulation as a zero-flux boundary.

The geometry of the systems is shown in Figure 8.

t

t

soil surface 15 cm

soil

profile zero flux

100 cm

free drainage

Figure 8. Outline of the geometry of the two simulated systems.

Conservation of mass in the simulations was achieved reasonable well on daily basis (cf.

Figure 9 and 12) with Az = 1 cm, ö2 = 0.01 cm and an allowable number of iterations equal to 12 within every hourly output step. This implies an initial time step of 1 hour decreasing to a minimum of 0.9 sec. The accumulated mass balance error to time t Err1 (cf. Figures 9b and 12b) was calculated from equation 67.

In figures 9-11 and 12-14 outputs from the simulations are shown for the 1 m soil profile and the micro-lysimeter, respectively, when E^ was set to 0.5 cm day'1.

Total evaporation after 10 days was 24 mm and 20 mm for the soil profile and the micro- lysimeter, respectively, cf. Figure 9c and Figure 12c.

(67)

0.

o. - o.

r 1.0 0.5 £

c. c.

(U

0.0 o

-0.5 “

oL>

<

Water content Total flux

Figure 9. Simulations for the 1 m soil profile with = 0.5 cm day'1, a. Accumulated water loss from the soil profile was calculated from integrated soil water profiles or from the sum of the upper and lower actual fluxes, b. Accumulated mass balance error, c. Accumulated actual evaporation.

roi

roE (J CCD

-i~>

CO o

C_0)

0.5 cm 3.5 cm 34 .5 cm

1 .5 cm B . 5 cm 6 5 . 5 cm

2 . 5 cm 14 .5 cm 9 9 . 5 cm

Hours

Figure 10. Simulated soil-water contents for the soil profile at different depths. = 0.5 cm day'1.

Soildepth,cm

---- o ---- 12 Water content, c m3 c m ' 3

--- 3 6 --- 8 4 --- 1 3 2

---- 180 .... 240 hours

Figure 11. Simulated soil-water profiles for the soil profile at different times. = 0.5 cm day1.

r

1 . 0

0.5 £ C.C- Q) 0.0 .

U3 -0.5 “

EE O)O) C(O c o c.

-p0)

<0 2 3 UO

<

1— i— I— '— I— '— I— '— i— 1— i— 1— i— 1— i 1 I 1 I ' r

0 24 48 72 96 120 144 168 192 216 240 Hours

Water content Total flux

Figure 12. Simulations for the micro-lysimeter with = 0.5 cm day'1, a. Accumulated water loss from the soil profile was calculated from integrated soil water profiles or from the sum of the upper and lower actual fluxes, b. Accumulated mass balance error, c. Accumulated actual evaporation.

---- 0.5 cm ---- 1.5 cm ---- 2.5 cm H o u r s

---- 3.5 cm ---- 4 . 5 c m ---- 6.5 cm

Figure 13. Simulated soil-water contents for the micro-lysimeter at different depths. = 0.5 cm day'1.

O ---- 12 Water content, crn^ cm 3

---- 36 ---- 84 --- 132

---- 180 .... 240 hours

Figure 14. Simulated soil-water profiles for the micro-lysimeter at different times. = 0.5 cm day1.

In Figure 15 is shown, at high evaporative demand the hourly difference (15a) and the accumulated difference (15b) between simulated actual evaporation from the soil profile and the micro-lysimeter. From Figure 15b it is seen that the duration for the micro-lysimeter is around 3 days.

0.3-i

0

-1

^ L 1^ lU L L L I

"i 1 i 1 i 1 i 1 i ' i 1 r

96 120 144 168 192 216 240 Hours

Figure 15. a. Calculated difference between actual evaporation from the soil profile Ea_s and the micro-lysimeter Ea_m. b. Accumulated difference. = 0.5 cm day'1.

In Figure 16 is shown at low evaporative demand the hourly difference (16a) and the accumulated difference (16b) between actual evaporation from the soil profile and the micro- lysimeter. From Figure 16b it is seen that the duration for the micro-lysimeter is around 6 days.

Hours

Figure 16. a. Calculated difference between actual evaporation from the soil profile Ea_s and the micro-lysimeter Ea_m. b. Accumulated difference. = 0.2 cm day'.

5. Discussion

From the simulations it was found that the duration of a micro-lysimeter depends on the evaporative demand. When the soil was initially wet the duration was found to be around 3 days under an evaporative demand of 5 mm day' and 6 days under an evaporative demand of 2 mm day'1. These findings illustrate the usefulness of model simulations, allowing for quick and low cost answers to questions on evaporation from micro-lysimeters — but precaution must be taken if the results are applied in reality. This is because a model is a picture of the reality and in particular because this model has not been validated against measured data on bare soil evaporation. As mentioned, only few processes are included in the model — conditions which possibly could introduce a certain bias in the calculations.

The empirical elements in the model describing the conditions at the upper boundary (cf.

Section 3.2) probably introduce some errors in the calculation of actual evaporation.

Calculation of hourly potential evaporation from equation 10 leads to negative values (i.e.

dewfall) at nighttime, which is unrealistic under some climatic conditions. This causes, in the model, an infiltration of water into the soil at nighttime, which increases the duration of the micro-lysimeter.

To assess a more accurate and reliable estimation of the duration of a micro-lysimeter the model has to be applied using measured input variables and the calculated output has to be validated against field measurements.

The results are, however, comparable to measurements of Sadras et al. (1991), who found a lower evaporation from micro-lysimeters compared to the surrounding soil when the micro- lysimeters were removed for weighing after a week. Shawcroft and Gardner (1983) found similar cumulative water loss from micro-lysimeters and the surrounding soil within 5 and 12 days in a field experiment in 1975 and 1976, respectively. They argued these findings to be caused possibly by several compensating errors due to unknown drainage, upward flow and plant water uptake in the surrounding soil and the restricted flow in the micro-lysimeters.

References

Berge, H .F.M. ten. (1990). Heat and water transfer in bare topsoil and the lower atmosphere.

Simulation Monographs 33. Wageningen. 207 pp.

Boast, C.W. and Robertson, T.M. (1982). A "micro-lysimeter" method for determining evaporation from bare soil: Description and laboratory evaluation. Soil Sei. Soc. Am.

J. 46, 689-696.

Chapra, S.T. and Canale, P.C. (1988). Numerical methods for engineers. 2. ed. McGraw-Hill Book Company. New York. 812 pp.

Feddes, R.A., Kowalik, P.J. and Zaradny, H. (1978). Simulation of field water use and crop yield. Simulation Monographs. Wageningen. 189 pp.

Hanks, R.J. (1991). Infiltration and redistribution. In: Modelling Plant and Soil Systems, Hanks, J. and Ritchie, J.T. (Eds.). Agronomy Monograph no. 31. p 181-204.

Hansen, S., Jensen, H.E., Nielsen, N.E. and Svendsen, H. (1990). DAISY — Soil plant atmosphere system model. No. A10. The National Agency of Environmental Protection.

Ministry of the Environment. Copenhagen. 272 pp.

Jensen, K.H. (1983). Simulation of water flow in the unsaturated zone including the root zone. Series Paper No. 33. Institute of Hydrodynamics and Hydraulic Engineering, Technical University of Denmark. 259 pp.

Mahfouf, J.F. and Noilhan, J. (1991). Comparative study of various formulations of evaporation from bare soil using in situ data. J. Appl. Meteor. 30, 1354-1365.

Sadras, V.O., Whitfield, D.M. and Connor, D.J. (1991). Regulation of evapotranspiration, and its partitioning between transpiration and soil evaporation by sunflower crops: a comparison between hybrids of different stature. Field Crops Research 28, 17-37.

Shawcroft, R.W. and Gardner, H.R. (1983). Direct evaporation from soil under a row crop canopy. Agricultural Meteorology 28, 229-238.

van Genuchten, M. T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sei. Soc. Amer. J. 44, 892-898.

Zaradny, H. (1978). Boundary conditions in modelling water flow in unsaturated soils. Soil Sei. 125(2), 75-82.

A. List of used symbols

a Soil specific constant, cm'1 C(ijr) Specific water capacity, cm'1 6, Convergence limit

52 Convergence limit, cm

E, Actual evaporation, cm3 cm'2 hour'1

E!,+' Actual evaporation at time ti+l, cm3 cm'2 hour'1 Ep Potential evaporation, cm hour'1

Ej,+I Potential evaporation at time ti+l, cm hour'1 Fv Potential evaporation, cm day'1

Err' Accumulated mass balance error to time t, mm G Soil heat flux density, Wm'2

H Sensible heat flux density, Wm'2 h Water tension, cm

K, Saturated hydraulic conductivity, cm3 cm'2 hour'1 K(i|r) Unsaturated hydraulic conductivity, cm3 cm'2 hour'1 A.E Latent heat flux density, Wm'2

M Iterations within a time step m,n Soil specific constants

q Water flux density, cm3 cm'2 hour'1

qj±Vi Water flux density at time t1 and depth zj±vl, cm3 cm'2 hour'1 qE Gravity flux density, cm3 cm'2 hour'1

q^+1 Gravity flux density at time tm , cm3 cm'2 hour'1 q, Maximum flux density, cm3 cm'2 hour'1

i|r Matric potential, cm

>|rj Matric potential at time t1 and depth zjt cm Rn Net radiation, Wm'2

t Time, hour t1 i'th time, hour At1 i'th time step, hour 9 Water content, cm3 cm'3

0- Water content at time t: and depth zj( cm3 cm'3 0O Water content at the soil surface, cm3 cm'3 03,0, Saturated and residual water content, cm3 cm'3 z Soil depth, cm

Zj j'th soil depth, cm Az Soil layer, cm

f

'å

fl

Afdelinger under Statens Planteavlsforsøg