Journal of Hydrology 373 (2009) 545–561

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Effects of meteorological models on the solution of the surface energy balanceand soil temperature variations in bare soils
Hirotaka Saito a,*, Jiri Šimůnek bab

Department of Ecoregion Science, Tokyo University of Agriculture and Technology, 3-5-8 Saiwaicho, Fuchu, Tokyo 183-8509, Japan
Department of Environmental Sciences, University of California, Riverside, California 92521, USAa r t i c l e

i n f oArticle history:
Received 3 July 2008
Received in revised form 14 March 2009
Accepted 23 May 2009This manuscript was handled by P. Baveye,
Editor-in-Chief, with the assistance ofStephanie Roulier, Associate Editor
Keywords:
HYDRUS1D
Water and vapor transport
Thermal conductivity
Net radiation
Evaporations u m m a r y
A complete evaluation of the soil thermal regime can be obtained by evaluating the movement of liquidwater, water vapor, and thermal energy in the subsurface. Such an evaluation requires the simultaneous
solution of the system of equations for the surface water and energy balance, and subsurface heat transport and water ﬂow. When only daily climatic data is available, one needs not only to estimate diurnal
cycles of climatic data, but to calculate the continuous values of various components in the energy balance equation, using different parameterization methods. The objective of this study is to quantify the
impact of the choice of different estimation and parameterization methods, referred together to as meteorological models in this paper, on soil temperature predictions in bare soils. A variety of widely accepted
meteorological models were tested on the dataset collected at a proposed low-level radioactive-waste
disposal site in the Chihuahua Desert in West Texas. As the soil surface was kept bare during the study,
no vegetation effects were evaluated. A coupled liquid water, water vapor, and heat transport model,
implemented in the HYDRUS-1D program, was used to simulate diurnal and seasonal soil temperaturechanges in the engineered cover installed at the site. The modiﬁed version of HYDRUS provides a ﬂexible
means for using various types of information and different models to evaluate surface mass and energybalance. Different meteorological models were compared in terms of their prediction errors for soil temperatures at seven observation depths. The results obtained indicate that although many available meteorological models can be used to solve the energy balance equation at the soil–atmosphere interface incoupled water, vapor, and heat transport models, their impact on overall simulation results varies. For
example, using daily average climatic data led to greater prediction errors, while relatively simple meteorological models may signiﬁcantly improve soil temperature predictions. On the other hand, while models for the albedo and soil emissivity had little impact on soil temperature predictions, the choice of theatmospheric emissivity models had a greater impact. A comparison of all the different models indicatesthat the error introduced at the soil atmosphere interface propagates to deeper layers. Therefore, attention needs to be paid not only to the precise determination of the soil hydraulic and thermal properties,
but also to the selection of proper meteorological models for the components involved in the surface
energy balance calculations.
Ó 2009 Elsevier B.V. All rights reserved.

Introduction
A complete evaluation of the movement of liquid water, water
vapor, and heat in the subsurface can be obtained by simultaneously solving the system of equations describing the surface
water and energy balance, and subsurface heat transport and water
ﬂow. The use of a coupled liquid water, water vapor, and heattransport model to simulate continuous changes in water contents,
soil temperatures, and a variety of ﬂuxes has been presented by,among many others, Nassar and Horton (1989), Noborio et al.
(1996a), Fayer (2000), and Saito et al. (2006). When simulations* Corresponding author. Tel./fax: +81 42 367 5584.
E-mail address: hiros@cc.tuat.ac.jp (H. Saito).
0022-1694/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved.
doi:10.1016/j.jhydrol.2009.05.019

are conducted at the ﬁeld scale, boundary conditions at the soil–
atmosphere interface for water, vapor, and heat transport are usually determined using the surface water and energy balance (e.g.,van Bavel and Hillel, 1976; Boulet et al., 1997). However, direct,
continuous measurements of all the components needed to fully
evaluate such surface mass and energy balances rarely exist. Usually, only standard daily climatic data from nearby weather stations and/or daily irrigation schedule are available. If detailed
predictions of water and heat ﬂuxes are needed, components ofthe surface mass and energy balance at much smaller time steps
than those taken at daily intervals will need to be evaluated using
daily standard climatic data.
The energy balance at the soil–atmosphere interface is expressed as

546

Rn  H  LE  G ¼ 0

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

ð1Þ

where Rn is the net radiation (W m2), H is the sensible heat ﬂux
density (W m2), LE is the latent heat ﬂux density (W m2), L is
the latent heat (J kg1), E is the evaporation (kg m2 s1), and G is
the surface heat ﬂux density (W m2). While Rn and G are positive
downward, H and LE are positive upward. To solve Eq. (1) for thesurface heat ﬂux G, which is needed as the upper boundary condition in the solution of the heat transport equation, continuous variations in Rn, H, and LE must be calculated or measured. When only
daily information is available, continuous values for various components in the energy balance equation must be obtained using existing estimation and parameterization methods.
Continuous diurnal cycles in climatic variables, such as air temperature, are usually generated from their mean daily values using
the analogy between their cycles during the day and trigonometricfunctions (e.g., Jury and Horton, 2003). Once the values for air temperature, relative humidity, precipitation, and wind speed at any
given time are obtained, they can be further used in parameterization formulas to calculate the components of the energy and waterbalance equations (Saito et al., 2006). It was not our intention to
separately consider above functions and formulas. Therefore, in
the remainder of this manuscript, the functions for generation ofcontinuous diurnal cycles of climatic variables and the parameterization formulas for components of the energy balance equation
are both referred to as ‘‘meteorological models”. A number of comparative studies, in which measured components were comparedto calculated ones, have appeared in the literature. For example,
Ortega-Farias et al. (2000) compared measured and predicted air
emissivities. As there are a number of available and accepted meteorological models to calculate the atmospheric variables, it is hardto determine which model is most suitable for a particular application. In addition, an extensive amount of work has been carried out
over the last few decades to develop meteorological models that
accurately predict evaporation rates from the soil and vegetation(e.g., Brutsaert, 1982; Monteith and Unsworth, 1990). Because of
the high complexity of both air and subsurface conditions, it is extremely difﬁcult to choose the best model to estimate evaporation
rates from a particular soil. To our knowledge, the impact of the
choice of particular meteorological models on simulated waterﬂow and heat transport in the vadose zone has, so far, not been discussed or investigated. Since all the components of the surface
mass and energy balance affect each other more or less in the coupled water ﬂow and heat transport model, investigating their mutual interactions is not a straightforward task. A change in one
variable can easily alter all the others.
Thus, the main objective of this study is to quantify the impactof the choice of particular meteorological models on the predictionof bare soil temperatures. We use various meteorological models inthe surface energy balance equation, and then evaluate how thesemethods affect soil temperatures calculated with the coupled liquid water, water vapor, and heat transport model. We compare
predicted soil temperatures at different soil depths with measured
values. The results will allow investigators and/or practitioners toevaluate their choice of meteorological models, and provide aquantiﬁed assessment of the effects these models have on predictions of soil temperatures.
A variety of meteorological models are reviewed and tested inthis study, using a dataset collected at a proposed low-level radioactive-waste disposal site in the Chihuahua Desert in West Texas,
10 km east of Sierra Blanca, where prototype engineered covers
were installed (Scanlon et al., 2005). The energy balance assessment in the engineered covers is as important for evaluation of
their performance as the mass (water) balance analysis. Whilethe long-term water balance of the site was evaluated by Scanlon
et al. (2005), who showed that a capillary barrier can signiﬁcantly

reduce drainage in arid and semi-arid regions, the energy balanceof the site has not yet been fully assessed. The coupled liquidwater, water vapor, and heat transport model, based on the modiﬁed HYDRUS-1D software package (Saito et al., 2006), is used inthis study to simulate soil temperatures in the engineering cover,the surface of which was kept bare during the analyzed time period. The modiﬁed version of HYDRUS provides a ﬂexible way to use
various types of climatic information to evaluate surface mass andenergy balance when continuous changes in water contents, temperatures, and ﬂuxes are simulated.
Method description
Generating diurnal cycle of climatic data
The solution of the energy balance equation (Eq. (1)) at a timeinterval of interest requires knowledge of the values of climatic
variables such as air temperatures, atmosphere relative humidities,
and wind speeds at the same or similar time intervals. However,weather stations do not always provide standard data at timeintervals of interest. Thus, diurnal changes in these variables need
to be calculated from available daily average values using meteorological models (e.g., Ephrath et al., 1996). In this study, we compared relatively simple approaches for generating the diurnal
cycles of climatic variables from available daily information.
Air temperatureContinuous values for air temperature, Ta, can be obtained fromthe daily maximum and minimum air temperatures usually available from the weather station using a trigonometric function witha period of 24 h as follows (Kirkham and Powers, 1972):

 t  13
T a ¼ T þ At  cos 2p24

ð2Þwhere T is the average daily temperature (°C), At is the amplitude of
the cosine wave (°C) calculated from the difference between thedaily maximum and minimum temperatures, and t is a local time
during the day (h). The argument of the cosine function shows thatthe highest temperature is assumed to occur at 1 p.m. and the lowest at 1 a.m.
Wind speed
It is well known that wind transports heat and effectively mixes
the soil–atmospheric boundary layer. Wind is generally highly variable in speed and direction, since it involves mostly turbulent ﬂow
characterized by random ﬂuctuations (Campbell, 1977). Because
several parameterization formulas for components in the energybalance equation require continuous inputs of the wind speed, continuous diurnal changes of the wind speed must somehow be calculated even when daily average values are all that is available. At
present we have two simple approaches, in addition to using a constant daily value for the entire day. Both approaches use the following maximum to minimum wind speed ratio Ur, which is deﬁned as

Ur ¼

U max
U min

ð3Þ

where Umax and Umin (m s1) are the unknown minimum and maximum wind speeds of the day, respectively. This ratio may be determined from prior knowledge or calibration using available data. The
maximum and minimum wind speeds can then be calculated fromthe daily average wind speed as follows:

2U r
U
1 þ Ur2
U min ¼U
1 þ Ur
U max ¼

ð4Þ
ð5Þ

547

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

where U (m s1) is the daily average wind speed.
FAO (1990) simply used the maximum and minimum values forday-time (0700–1900) and night-time (1900–0700) wind speeds.
On the other hand, Gregory (1989) modeled the cyclic behaviorof the wind speed Uh during the day as follows: t  tmax
U h ¼ U þ ðU max  UÞ  cos 2p24

ð6Þ

where t is the time of the day and tmax is the time when the maximum wind speed occurs. In Gregory’s study, the maximum windspeed occurred at 3:30 pm.
Relative humidity
Since typical diurnal air temperature variations show a cyclic
behavior throughout the day, it is reasonable to assume that therelative humidity also shows such a cyclic pattern. As for the airtemperature (Eq. (2)), a trigonometric function can be used to calculate continuous values of the relative humidity from the daily
information (e.g., Gregory et al., 1994). When both the daily maximum and minimum relative humidity values are available, diurnal
cycles of relative humidity can be generated using a cosine function with a period of 24 h as follows:  þ Ar  cos 2p t  t max
Hr ¼ H24

ð7Þ

 is the average daily relative humidity (–), Ar is the ampliwhere H
tude of the cosine wave calculated from the difference between the
maximum and minimum relative humidity values, t is the local
time within the day, and tmax is the hour of the day when the relative humidity is at its maximum. The relative humidity is generally
at its minimum during day-time when the air temperature and thewind speed are at their maxima. Gregory et al. (1994) showed thatthe relative humidity was highest between 5 and 6 a.m. in Lubbock,
Texas. In this study, we assume that the relative humidity peaks at 5a.m.
In many cases, only daily average relative humidity values areavailable. The maximum and minimum relative humidity values
can then be speciﬁed similar to those for wind speed, i.e., using
the ratio of maximum and minimum values (Eq. (3)). After this ratio is determined, the maximum and minimum relative humidity
values can be calculated, as with wind speed, using Eqs. (4) and(5). However, the maximum relative humidity value needs to beconstrained so as not to exceed 100%. All values exceeding 100%
need to be reset to 100%. The continuous diurnal cyclic change inthe relative humidity can then be obtained using Eq. (7), in which
the amplitude of the cosine wave can be obtained from the difference between the calculated maximum and minimum relative
humidity values.
Parameterization of components in surface energy balance equation
Surface precipitation, irrigation, evaporation, and heat ﬂuxes
are used as boundary conditions for simulating soil water contents
and temperatures in bare ﬁeld soils. Evaporation and heat ﬂuxes
can be calculated from the surface energy balance, Eq. (1). Net radiation, Rn, is deﬁned as (e.g., Brutsaert, 1982)

Rn ¼ Rns þ Rnl ¼ ð1  aÞSt þ ðes Rld#  Rlu" Þð8Þ

where Rns is the net shortwave radiation (W m2), Rnl is the net
longwave radiation (W m2), a is the surface albedo (–), St is theincoming (global) shortwave solar radiation (W m2), es is the soil
surface emissivity (–) representing the reﬂection of longwave radiation at the soil surface, Rld; is the incoming (thermal) longwave
radiation at the soil surface (downward ﬂux) (W m2) as emitted
by the atmosphere and cloud cover, and Rlu" is the sum of the out-

going (thermal) longwave radiation emitted from the surface (vegetation and soil) into the atmosphere (W m2).
Shortwave radiation
The value of the incoming shortwave solar radiation St (W m2),
at any given time and location, can be calculated by taking into account the position of the sun (e.g., van Bavel and Hillel, 1976) using
the following equation (Campbell, 1985):

St ¼ maxðGsc T t sin e; 0Þð9Þ2where Gsc is the solar constant (1360 W m ) and Tt (–) is the atmospheric transmission coefﬁcient, which is deﬁned as the ratio of themeasured daily global solar radiation Stm (W m2) and the daily potential global (extraterrestrial) radiation Ra (W m2):

Tt ¼

Stm
Ra

ð10Þ

The last term of Eq. (9), e (rad), is the solar elevation angle givenby (Monteith and Unsworth, 1990):

sin e ¼ sin u sin d þ cos u cos d cos

2pðt  t 0 Þ
24

ð11Þ

where d is the solar declination (rad), u is the latitude (rad), t is the
local time within a day, and t0 is the time of solar noon.
AlbedoThere are several accepted models to determine the surface albedo, a. A constant value can be applied for given surface conditions (e.g., 0.23 for grass). A simple linear relation relating albedo
with the surface water content may also be used, since albedo depends, especially for bare soils, on the soil surface wetness. The following formula was proposed by van Bavel and Hillel (1976) toestimate the albedo as a function of the water content:

a ¼ 0:25

h0 < 0:1

a ¼ 0:35  h0

0:1 6 h0 < 0:25

a ¼ 0:10

h0 P 0:25

ð12Þ

where h0 is the water content at the ground surface. In this function,
the albedo is constant for low and high water contents, and varies
linearly for intermediate water contents. On the other hand, in theirstudy Kondo et al. (1992) used loam with a slightly different formula, where the albedo depends on the surface water content, evenfor dry conditions:a ¼ 0:24  0:21h0h0 < 0:14

a ¼ 0:35  h0

0:14 6 h0 < 0:22

a ¼ 0:13

h0 P 0:22

ð13Þ

Longwave radiation
Both incoming and outgoing thermal longwave radiations Rld
and Rlu in Eq. (8) need to be evaluated in order to obtain continuous
net longwave radiation values. At a given temperature T (K), thermal longwave radiation is given by the Stefan–Boltzmann law (e.g.,Jury and Horton, 2003):

X¼ erT 4

ð14Þ
8

where r is the Stefan–Boltzmann constant (5.67  10 W m K4
or 4.89  109 MJ m2 K4 d1), and e is the emissivity (–), which isequal to 1 for a blackbody and has values between 0 and 1 for othersurfaces. The net longwave radiation at any given time t can then be
rewritten as (van Bavel and Hillel, 1976):

Rnl ðtÞ ¼ es ea rT 4a  es rT 4s

2

ð15Þ

where subscripts a and s are used for the atmosphere and the soil,
respectively.

548

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

The emissivity of the soil surface, es, generally depends on thewater content and vegetation of that surface. Although different
empirical models are available for the emissivity of a soil surface
with vegetation (i.e., surface emittance), soil surface emissivity
usually has a value close to 1 (e.g., Ortega-Farias et al., 2000; Fuchs
and Tanner, 1966). For bare soils without vegetation, the following
two empirical models have been commonly used in the literature.
The ﬁrst formula (van Bavel and Hillel, 1976) gives the surface
emittance of 0.9 for a dry surface and 0.98 for a saturated soil when
the saturated water content is 0.45:

es ¼ 0:90 þ 0:18hs

ð16Þ

Noborio et al. (1996b) used the following formula for the loamy
sand, which gives nearly constant values for the surface emissivity
for a range of water contents:2

es ¼ 0:898 þ 2:17  10 hs

ð17Þ

Since water vapor and carbon dioxide in the air are the main
emitting gases, the atmospheric emissivity ea depends on manyother factors, such as the air temperature and the air humidity.
In general, more than half of the longwave radiation received at
the ground comes from gases in the lower 100 m of the atmosphere, and almost 90% comes from the ﬁrst 1 km height (Monteith
and Unsworth, 1990). Therefore, the longwave radiation can be calculated accurately if the proﬁles of the air temperature and the air
humidity are known. However, that is rarely the case. To cope withthis problem, a large number of parameterization methods havebeen developed empirically to estimate the atmospheric emissivity
from the near-surface temperature and/or humidity, which areboth commonly available and easily measured. While one set of
formulas depends on both the surface air temperature and the surface air humidity, the other set is formulated as a function of the airtemperature only. Some models of the former type are expressed
as a function of the vapor pressure, ea (kPa), or the vapor density,qva (kg m3), both of which can be expressed conveniently as afunction of the air temperature and the relative humidity of theair. The relative humidity is deﬁned as either the ratio of the actual
vapor pressure to the saturation vapor pressure, or the ratio of theactual vapor density to the saturation vapor density at a given
temperature.
One of the ﬁrst empirical models for calculating atmosphericemissivity was suggested by Brunt (1932). In this model, ea depends on the near-surface vapor pressure ea (millibars) as follows:

pﬃﬃﬃﬃﬃ

ea ¼ 0:51 þ 0:066 ea

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

ð19Þ

One of the more-commonly used models is a power function
originally derived by Brutsaert (1975):

ea ¼ 1:24  1=7eaTa

ð20Þ

where Ta is the near-surface air temperature (K). While Brutsaert’s
equation was derived from the physics of the process, the Brunt’s
equation is purely empirical. Because Brutsaert’s model is fairly
insensitive to changes in the air temperature, it can be simpliﬁed
as follows (Brutsaert, 1975):

ea ¼ 0:553  ðea Þ1=7

ea ¼ 0:179ea1=7 expð350=T a Þ

ð21Þ

ð22Þ

Idso (1981) also proposed the following formula to avoid having
a zero emissivity at the zero intercept in the above formula. The revised model has a linear dependency of ea on ea:ea ¼ 0:70 þ 5:95  105 ea expð1500=T a Þ

ð23Þ

Satterlund (1979) proposed the following equation to improvethe estimation of atmospheric emissivity when the air temperature
is below 0 °C:

hnea ¼ 1:08 1  exp ðea ÞT a =2016

oi

ð24Þ

Other formulas that depend solely on the near-surface air temperature, Ta, were derived by Swinbank (1963) and Idso and Jackson (1969), respectively:ea ¼ 9:2  106 T 2an4

ea ¼ 1  0:26 exp 7:77  10 ð273  T a Þ

2

o

ð25Þ
ð26Þ

Evaporation and sensible heat ﬂux
Since surface evaporation is controlled by atmospheric conditions, surface moisture, and water ﬂow in the soil, a formulaaccounting for all of these factors is needed. Evaporation, E, in arid
or semi-arid regions is often calculated as (Camillo and Gurney,1986):q  qv a
E ¼ vsra þ rs

ð27Þ

where qvs (kg m3) is the water vapor density at the soil surface, qva(kg m3) is the atmospheric vapor density, ra (s m1) is the aerodynamic resistance to water vapor ﬂow from the soil surface to theatmosphere, and rs (s m1) is the soil surface resistance to water vapor ﬂow in soil pores acting as an additional (source of?) resistancecoupled with the aerodynamic resistance (Camillo and Gurney,1986).
Four variables need to be determined to calculate the evaporation rate using Eq. (27): qva, qvs, ra, and rs. The atmosphere vapor
density qva (kg m3) is obtained as:

ð18Þ

Empirical coefﬁcients in Brunt’s formula were originally derived
for a speciﬁc location. However, this relationship can be used forother locations, and can be thought of as a more-or-less universal
equation (Hatﬁeld et al., 1983). More recently, Chung and Horton(1987) used a similar formula that is a function of the vapor density, qva (kg m3):

ea ¼ 0:605 þ 0:048 1370  qv awhich corresponds to a typical air temperature of 285 K. Idso (1981)
proposed a similar equation, in which the atmospheric emissivity
depends on both Ta (K) and ea (millibars), while also keeping thedependence:e1=7a

qv a ¼

RH
 q
100 v að28Þ

where RH is the relative humidity of the air (%) and qv a is the saturated vapor density (kg m3), which can be expressed conveniently
as a function of the air temperature (see for example, Campbell,1977). The water vapor density at the soil surface is calculated following the equilibrium model of Philip and de Vries (1957) as:

qv s ¼ exp




hMg
 qv sRT sð29Þwhere qv a is the saturated vapor density at the soil surface (kg m3),
h is the pressure head (m), M is the molecular weight of water(M mol1), g is the gravitational acceleration (m s2), Ts is the soil
surface temperature (K), and R is the universal gas constant(J mol1 K1) (8.314 J mol1 K1).
The aerodynamic resistance, ra, at the soil–atmosphere boundary is expected to vary depending upon the wind speed and the level of turbulence, which are determined by the roughness of thesoil surface, the distance from the surface, and the condition of
the atmosphere above the surface. The aerodynamic boundary
layer resistance is commonly calculated as (Campbell, 1985):

549

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

ra ¼







zref  d þ zH
1
ln
þ wH
2
zH
uk
 


zref  d þ zm
 ln
þ wm
zm



ð30Þ

where k is von Karman’s constant (=0.41), zref is the reference heightof measurements (for both temperature and wind speed) (m), u is
the mean wind speed (m s1) at height zref, d is the zero plane displacement (m), zH is the surface roughness for the heat ﬂux (m), zm
is the surface roughness for the momentum ﬂux (m), wH is theatmospheric stability correction factor for the heat ﬂux (–), andwm is the atmospheric stability correction factor for the momentum
ﬂux (–). For bare soils, the zero plane displacement, d, is equal tozero, while typical surface roughness values of 0.001 m are used
for both zH and zm (Oke, 1978). Eq. (30) can be rewritten as:

ra ¼



1
2


ln

uk
zref  zH ; zm


  


zref  dzref  dþ wH  ln
þ wm ;zHzm
ð31Þ

when zref is much greater than the surface roughness (Brutsaert,
1982). In the remainder of the study, calculations are made using
Eq. (31), because zref (1.5 m) is much greater than the typical roughness lengths for bare soils (e.g., 0.001 m).
The aerodynamic resistance depends on the so-called stability
condition of the atmosphere, which can be assessed using Monin–Obukhov’s stability parameter (or the MO length). The MO
length is calculated as follows (e.g., Camillo and Gurney, 1986; Aluwihare and Watanabe, 2003):

MO ¼ 

qa ca T a U 3

ð32Þ

kgH

where H is the sensible heat ﬂux at the soil surface, Ta is the atmosphere temperature (K) at zref, and U* is the frictional velocity deﬁned based upon the logarithmic wind proﬁle law (e.g., Aluwihare
and Watanabe, 2003; Camillo and Gurney, 1986; van de Griend
and Owe, 1994; Brutsaert, 1982):

 

1
zref  d
U  ¼ uk ln
þ wm
zm

ð33Þ

The evaluations of wm and wh are then determined using theatmospheric stability parameter, f, deﬁned as (Brutsaert, 1982):

f¼

zref  dMO

ð34Þ

For the neutral atmosphere, deﬁned by the condition|Ta  Ts| 6 0.1, the stability correction factors wm and wh become
equal to 0 (e.g., Camillo and Gurney, 1986; Choudhury et al.,1986). For the unstable atmosphere (Ta < Ts or MO < 0), the stability
correction factors are empirically expressed as follows (e.g., Paulson, 1970; Brutsaert, 1982):

wh ¼ 2 ln

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1 þ 1  16f
2

!

ð36Þ

and for the stable atmosphere (Ta > Ts or MO > 0), the stability factors for the heat and momentum ﬂuxes can be simply expressed as(for moderately stable, f < 1, and extremely stable conditions, f > 1):
wh ¼ wm ¼

5f 0<f<1
5

f>1

MO ¼ 

qa ca T a U 3
kgH
T a u2

gðT a  T s Þ ln¼

T a u2
gðT a  T s Þ ln



zref
z0



1 þ lnðzwh=z0 Þrefð1 þ lnðzwm=z0 Þ Þ2
ref



ð38Þ

zref
z0

where zh and zm are assumed to be equal to the roughness length, z0.
They claimed that the right-hand-side portion of the third term of
the above equation is always between 0.97 and 1.0 for their data,
and may thus be set equal to 1.0. This results in there being no needfor the iterative procedure.
A number of studies use the Richardson number, Ri, which is analternative stability parameter (e.g., van Bavel and Hillel, 1976), in
calculations of the aerodynamic resistance. The Richardson number can be estimated from (e.g., Brutsaert, 1982; Beljaars and Holtslag, 1991):

Ri ¼

g  ðzref  dÞ  ðT a  T s Þ
T a  u2

ð39ÞLike the aforementioned deﬁnition of the MO length, the atmosphere is stable when the Richardson number is positive, andunstable when the Richardson number is negative. The aerodynamic resistance is then deﬁned using the neutral resistance, raa(when wm and wh are 0), and the stability correction w as:

ra ¼ r aa  w

ð40Þ

Koivusalo et al. (2001) used the following stability correction
parameter for stable and unstable conditions:w ¼ ð1  10Ri Þ1

Ri < 0ð41Þ

w ¼ ð1 þ 10Ri Þ Ri > 0

Unlike cases where the MO length is used, the Richardson number does not require an iterative approach, but allows direct estimation of the aerodynamic resistance. As can be seen from theirdeﬁnitions, the MO length (Eq. (32)) and the Richardson number(Eq. (39)) are strongly related, and there are a number of studies
showing their relationship (e.g., Monteith and Unsworth, 1990).
On the other hand, the soil surface resistance, rs, depends
greatly on soil structure and texture. The relationship between rs
and the surface water content, h0, has been empirically formulated,typically using an exponential function (e.g., van de Griend andOwe, 1994), indicating that surface resistance increases dramatically as the soil dries out. A general formula for soil surface resistance was proposed by van de Griend and Owe (1994), forh0 6 hmin:

rs ¼ r 0 expðFðhmin  h0 ÞÞ
ð35Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1 þ ð1  16fÞ1 þ 1  16f
 ln
wm ¼ 2 ln
2
2
n
o p
þ 2 arctan ð1  16fÞ0:25 2
0:25

In general, calculations of ra require an iterative procedure, as
the stability factors are evaluated from variables that themselves
depend on the stability factors. To avoid the iterative approach, Camillo and Gurney (1986) approximated the MO length as follows:

ð37Þ

ð42Þ
3

3

where h0 is the soil water content (m m ) in the top 1 cm, r0 is theminimum surface resistance, which is theoretically equal to theresistance of the molecular diffusion across the water surface itself,
hmin is the empirical minimum water content above which the soil
is able to deliver vapor at a potential rate, and F is an empiricalparameter. The following equation was ﬁtted to the soil resistancedata with the minimum resistance occurring at a soil water content
of approximately 15% and an r0 value of 10 (s m1):


rs ¼

10:0  expð35:63ð0:15  h0 ÞÞ h0 6 0:15
10:0

h0 > 0:15

ð43Þ

Camillo and Gurney (1986) used the formula for the surface
resistance (s m1) originally proposed by Sun (1982):

550

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

r s ¼ 3:5

 2:3hsat
þ 33:5h0

ð44Þ

where h0 and hsat are the volumetric water contents in the top
0.5 cm of the soil, and at full saturation, respectively. One set of thisdata can also be represented very well by a linear model (correlation coefﬁcient = 0.98):

r s ¼ 805 þ 4140ðhsat  h0 Þð45ÞAdditional soil resistance formulas are cited by Mahfouf andNoilhan (1991). The Passerat de Silans’ model has an exponential
form as:r s ¼ A expðBh0 Þ

ð46Þ
3

3

where h0 is the soil water content (m m ) of the surface layer of
speciﬁed depth (0–20 mm or 0–50 mm) and A and B are empiricalparameters speciﬁc for the soil. In their study, the A and B parameters had values of 3.8113  104 and 13.515/hfc, respectively,
where hfc is the volumetric water content at ﬁeld capacity. On theother hand, Kondo et al. (1990) proposed a different formula forsand, which is also a function of the water content:

F1
rs ¼
ðhsat  h0 ÞF 2
Dm

ð47Þ

where h0 is the soil water content (m3 m3) of the top 20-mm layer,
F1 and F2 are ﬁtting parameters with values of 8.32  105 and 16.6,
respectively, and Dm is the molecular diffusion of water vapor,which has been formulated by Camillo et al. (1983) as:Dm ¼ D0



Ts
273:16

1:75¼ 0:229  104



1:75
Ts
273:16

ð48Þ

These parameter values were given by Kondo et al. (1990) for
sand, which has a saturated water content of 0.392.
It is widely accepted that the aerodynamic resistance to vapor
ﬂow is usually very close to the aerodynamic resistance to heatﬂow (e.g., Aluwihare and Watanabe, 2003; van Bavel and Hillel,1976). In this study, the sensible heat ﬂux is thus calculated using
the aerodynamic resistance (e.g., van Bavel and Hillel, 1976) as:

H ¼ CaTs  Tarað49Þ

where Ta is the air temperature (K), Ts is the soil surface temperature (K), and Ca is the volumetric heat capacity of air (J m3 K1).
Numerical simulation
The modiﬁed HYDRUS-1D code used in this study solves thevariably-saturated water ﬂow and heat transport equations simultaneously (See Appendix), using the method of ﬁnite elements forthe spatial discretization and ﬁnite differences for the temporal
discretization (Saito et al., 2006). When the surface energy balance(Eq. (1)) is solved, the code can use either directly observed climatic data or the continuous diurnal values generated from available daily information. In this study, continuous diurnal values ofclimatic data were generated from daily information collected at
the study site.

At each time step, the net radiation (Rn), the sensible heat ﬂux(H), and the latent heat ﬂux (LE) are calculated to obtain the surface
heat ﬂux (G) from the surface energy balance equation (Eq. (1)).
The surface heat ﬂux, G, is subsequently used as a known heat ﬂux
boundary condition. Calculations of Rn, H, and LE require knowledge of the temperature and pressure head at the soil surface. At
a given time step, the energy balance equation is solved ﬁrst withthe heat ﬂux at the boundary calculated using the soil moistureand temperature data obtained in a previous time step (old timestep). Soil water contents are then updated using the resulting
boundary water (evaporation) ﬂux by solving Richards equation.
The energy balance equation is then solved again, using updated
soil water contents to obtain updated heat and water ﬂuxes atthe boundary. Finally, with updated boundary conditions, waterﬂow and heat transport equations are solved. Convergence ischecked after new soil moisture and temperature data are
obtained.
Study site
The study site is a proposed low-level radioactive-waste disposal site in the Chihuahua Desert in West Texas, 10 km east of
Sierra Blanca (31°8.7730 N, 105°16.2370 E; elevation 1337 m),
where prototype engineering covers were installed (Scanlon
et al., 2005). Engineered cover designs, with a conductive or capillary barrier of sand at the 2-m depth, were installed at the site in
the summer of 1997. Water movement is generally restricted toonly those zones above the capillary barrier in this type of cover.
These covers also function as evapotranspiration (ET) covers, relying on surface vegetation to increase the water storage capacity of
the covers by removing water through ET. As a result, drainage todeeper depths can be minimized and negligible. There was, however, no vegetation during the time period analyzed in this study.
The cover consists of 0.3 m of topsoil (sandy clay loam, bulk
density of 1.5 Mg m3) underlain by 1.7 m of compacted soil (sandy clay loam, bulk density of 1.8 Mg m3). Gravel was added to theupper 0.3 m of the topsoil to reduce erosion. There is a 1.05 m capillary barrier underneath the compacted sandy clay loam layer,
consisting of four layers: sandy loam, muddy gravel, gravel, andsand. The thicknesses of each layer are 0.3, 0.3, 0.3, and 0.15 m,
respectively (Table 1).
Various hydrological and climatic data have been monitored at
the site since October 1997. The monitoring system consists of an
onsite weather station that monitors daily and hourly precipitation, wind speed, relative humidity, temperature, and solar andnet radiations at 1.5 m above the surface. However, to test the performance of various meteorological models, continuous diurnalchanges of climatic data were generated from daily informationin this study.
Soil temperatures and pressure heads were measured at seven
depths: 0.15, 0.3, 0.6, 0.9, 1.2, 1.5, and 2.0 m, using thermistors
and heat dissipation sensors (Scanlon et al., 2002). Soil hydraulicproperties of each soil layer were obtained from previous studies(Scanlon et al., 2002). Water retention data for soils in the upper
three soil layers were based on laboratory measurements using ahanging water column (pressure heads from 0 to 2 m) and a

Table 1
Soil textures and soil hydraulic parameters for each layer of the ET cover at the study site (Scanlon et al., 2002).
Layer

Thickness (m)TextureKs (cm/day)hshra (cm1)

n

1
2
3
4
5
6

0.3
1.7
0.3
0.3
0.3
0.15

Sandy clay loam
Sandy clay loam
Sand loam
Muddy gravel
Gravel
Sand

41
20
639
10
159840
587

0.45
0.35
0.40
0.14
0.51
0.38

0.00
0.00
0.00
0.00
0.00
0.00

0.027
0.010
0.020
0.007
10.95
0.050

1.276
1.167
1.464
1.188
1.722
1.774

551

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561
Table 2
Parameters of Chung and Horton’s (1987) soil thermal conductivity function for eachlayer of the ET cover at the study site.
Layer

b1

b2

b3

1–3
4–6

0.244
0.228

0.393
2.406

1.534
4.909

pressure plate apparatus (pressure heads from 1 to 50 m).
Water retention data and saturated hydraulic conductivities forunderlying soil layers were, on the other hand, obtained from various databases (e.g., UNSODA, Leij et al., 1996). Parameter values
from Scanlon et al. (2002) for all six layers are summarized in Table
1. HYDRUS-1D default values for thermal conductivities were used
for each layer of the soil proﬁle based on its textural class: the thermal conductivity of loam for layers 1–3 and that of sand for layers
4–6 (Table 2).
Soil proﬁle and initial conditions
The simulated soil proﬁle is 3.05 m deep, with six soil layers. Seven observation nodes were deﬁned in HYDRUS-1D, at depths
where soil temperatures and pressure heads were measured. The
soil proﬁle was discretized using 103 nodes, with nodal spacing
ranging from 0.2 cm near the soil surface and 2 cm around materialinterfaces to 15 cm within soil horizons. Simulations of the coupledtransport of liquid water, water vapor, and heat were conducted to
predict soil temperatures at each observation node for 365
consecutive days beginning October 1, 1997. There was novegetation at the soil surface (i.e., bare soil) during this period (inAugust 1998, ﬁve perennial warm-season bunch grass species
were planted). Initial conditions for temperatures and pressure
heads are based on a linear interpolation of measured values atthe site.

Boundary conditions
Boundary conditions at the soil surface for heat transport are
determined from the surface energy balance equation. At any given
time, the surface energy balance equation (Eq. (1)) is solved to obtain the surface heat ﬂux, G, which is used as a known heat ﬂux
boundary condition. Precipitation and evaporation rates were used
as upper boundary conditions for water ﬂow. Zero pressure andtemperature gradients were used as the bottom boundary conditions. These conditions assume that the water table is located far
below the domain of interest, and that heat transfer across the lower boundary occurs only by convection of liquid water and water
vapor.
Prediction performance
A set of meteorological models used in a previous study (Saito
et al., 2006) was maintained as a combination of models used forreference. Models used for the reference run are listed in Tables
3–5 as Scenario 0. The impact of the choice of different meteorological models on predicted soil temperatures was then evaluatedusing the following procedure.
1. Select models to generate continuous diurnal values of the airtemperature, the wind speed, and the relative humidity fromdaily climatic information.
2. Select models to calculate the soil albedo, the soil and atmospheric emissivities, and the aerodynamic and soil resistances.
3. For each selected model, specify parameter values.
4. Predict hourly soil temperatures numerically at seven observation nodes for the period of 365 days (from October 1, 1997 to
September 30, 1998) using the modiﬁed HYDRUS-1D code withselected meteorological models with speciﬁed parametervalues.
 Observation nodes are referred to as ON1 through ON7, from top(0.15 m depth) to bottom (2.0 m depth).

Table 3Models used to generate continuous diurnal variations in air temperature, wind speed, and relative humidity. Bold font indicates the method’s difference from the reference
scenario.
Scenario

Air temperatureWind speed

RatioaRelative humidity

Ratiob0 (Reference)12
3
4

Eq. (2)
Eq. (2)
Eq. (2)
Eq. (2)
Eq. (2)

Eq. (6)
Eqs. (4) and (5)Daily average
Eq. (6)
Eq. (6)

3.0
3.0
–
3.0
3.0

Eq. (7)
Eq. (7)
Eq. (7)
Eq. (7)
Daily average

3.0
3.0
3.0
2.0
–

a
b

Ratio of maximum and minimum wind speeds.
Ratio of maximum and minimum relative humidity.

Table 4Models used to calculate the albedo and the atmospheric and soil emissivities. For all scenarios listed in this table, the aerodynamic resistance model, ra, given by Campbell (1985)
is used. Bold font indicates the method’s difference from the reference scenario. (vB&H – van Bavel and Hillel (1976), C&H – Chung and Horton (1987), I&J – Idso and Jackson(1969)).
Scenario

Albedoeaes

0 (Reference)
5
6
7
8
9
10
11
12
13
14
15

Eq. (12) (vB&H)
Eq. (13) (Kondo et al., 1992)
Constant (0.23)
vB&H
vB&H
vB&H
vB&H
vB&H
vB&H
vB&H
vB&H
vB&H

Eq. (23) (Idso, 1981)
Idso
Idso
Eq. (18) (Brunt, 1932)
Eq. (19) (C&H)
Eq. (20) (Brutsaert, 1975)
Eq. (21) (Brutsaert, 1975)
Eq. (22) (Idso, 1981)
Eq. (24) (Satterlund, 1979)
Eq. (25) (Swinbank, 1963)
Eq. (26) (I&J)
Idso

Eq. (16) (vB&H)
vB&H
vB&H
vB&H
vB&H
vB&H
vB&H
vB&H
vB&H
vB&H
vB&H
Eq. (17) (Noborio et al., 1996b)

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

Table 5Models used to calculate the aerodynamic and soil surface resistances. Bold fontindicates the method’s difference from the reference scenario. (C&G – Camillo andGurney (1986), vG&O – van de Griend and Owe (1994), M&N – Mahfouf and Noilhan(1991)).
Scenariorars

0 (Reference)
16
17
18
19
20
21

MO, Eq. (38) (C&G)
MO, Eq. (32) (Iterative approach)
Ri, Eq. (41) (Koivusalo et al., 2001)
C&G
C&G
C&G
C&G

Eq. (45) (C&G)
C&G
C&G
Eq. (44) (C&G)
Eq. (42) (vG&O)
Eq. (46) (M&N)
Eq. (47) (Kondo et al., 1992)

35
30

Temperature [oC]

552

25
20
15
10
Generated
Observed

5
0

0

2

4

6

8

10

12

Day

 Day 1 corresponds to October 1, 1997 in the remainder of themanuscript.
5. At each observation node, xj, compute the correlation coefﬁcient
between predicted and measured temperatures, and prediction
errors: the mean square error (MSE) and the mean absolute
error (MAE), which are deﬁned, respectively, as follows:

Fig. 1. Hourly variations of observed air temperatures at the study site (opencircles) and those approximated from daily maximum and minimum values (solidline).

8N

j
1 X
jT  ðxj ; t i Þ  Tðxj ; t i Þj
MAEj ¼
Nj i¼1

ð51Þ

where Nj is the number of measurements at the xj node, T*(xj, ti)represents the temperature measured at time ti at the xj node,
and T(xj, ti) is the corresponding predicted temperature.
 MSE is used to assess the overall prediction performance, while
MAE is used to check discrepancies between observed and predicted temperature values.
6. Change one model at a time, and repeat from 1.
Since it is impractical to test all possible combinations of available meteorological models, we focus mainly on investigating what
impact a change in a particular model has on soil temperature predictions, and not on ﬁnding the best combination of these models.
We thus limit the number of combinations by starting with the reference combination (Table 3), and allowing only one model tochange at a time. Table 3 shows models used to generate continuous diurnal changes in the wind speed and relative humidity. Models used in the calculation of the net radiation are listed in Table 4,
while those used to obtain the surface resistances (soil and atmospheric) are summarized in Table 5. The remaining parameters,such as soil hydraulic and thermal properties, and initial andboundary conditions, are the same for the different scenarios.
Results and discussion
Generation of climatic data
Continuous diurnal changes in the air temperature can be generated using a sinusoidal function with a period of 24 h from daily
maximum and minimum data (Kirkham and Powers, 1972). The
formula (Eq. (2)) includes a term to allow the highest temperature
to occur at 1 p.m. and the lowest at 1 a.m. Fig. 1 depicts the ﬁrst 12
days of hourly measured air temperatures and those approximatedfrom daily maximum and minimum values. Calculated air temperatures match observed values well, except for day 11. Because daily values were measured separately, they may not represent thoseof day 11 due to some problems in measurement.
Fig. 2 shows the ﬁrst 12 days of diurnal wind speed hourly variations measured at the site, along with those calculated from theWind speed [m/sec]

ð50Þ

Daily
FAO
Gregory
Observed

7

N

j
1 X
MSEj ¼½T  ðxj ; t i Þ  Tðxj ; ti Þ 2
Nj i¼1

6
5
4
3
2
1
0

0

2

4

6
Day

8

10

12

Fig. 2. Hourly variations of observed wind speeds at the study site (open circles)
and those calculated from daily mean values using different models. In addition todaily mean values, the FAO model with day-time to night-time ratios of 2, and theGregory model with the maximum to minimum ratio equal to 3 are compared. Timeof the maximum wind speed is set to 3 pm.

daily mean wind speed using different models. The simplest approach is to use the constant daily mean value. This approach obviously understates temporal variations of the wind speed, and may
not be appropriate to simulate hourly variations of the soil temperature. When the FAO model (Eqs. (4) and (5)) is used, the day-time
to night-time ratio needs to be speciﬁed. FAO recommends a ratio
of about 2 for average conditions, which leads to a day-time windspeed 1.33 times greater than the mean value. Following the previous study of Saito et al. (2006), a ratio of 3 was used in this study.
When Ur = 3, the day-time wind speed is 1.5 times greater than themean value. Generated wind speeds approximate the observed values reasonably well, although detailed ﬂuctuations were not modeled (Fig. 2). An obvious disadvantage of the FAO model is that
there is always a sudden change in wind speed at 7 a.m. and 7
p.m., which may cause undesired consequences when the obtained
value is used to calculate other climatic variables.
On the other hand, the Gregory’s model allows for smoothly
varying wind speed variations by applying a trigonometric function (Eq. (6)). The maximum to minimum wind speed ratio, equal
to 3, was applied as in the FAO model above. Calculated windspeeds indeed change smoothly and approximate better observed
wind speeds. Detailed ﬂuctuations were again not modeled. Reproduction of such ﬂuctuations can never be obtained using these simpliﬁed methods. Overall, both the FAO and Gregory’s models
approximate the general trend of the wind speed throughout the
day from the mean daily value well.

553

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

Atmospheric emissivity [-]

1

80

60

40
Daily
Cyclic: R atio=3
Cyclic: R atio=2
Observed

20

0.9
0.8
0.7
0.6
0.5
0.4

Relative humidity = 10%
0

2

4

6
Day

8

10

0.3

12

Fig. 3. Hourly variations of observed relative humidities at the study site (opencircles) and those calculated from daily mean values. The maximum to minimum
ratios of 2 and 3 are used with the peak time set to 5 a.m. Interpolated daily meanvalues are also displayed.

The generation of hourly variations of the atmospheric relative
humidity from daily mean values requires a maximum to minimum relative humidity ratio that is similar to the one required
for wind speed generation (Eq. (3)). Following the analysis of windspeed presented above, relative humidities calculated using ratiosof 2 and 3 were compared with daily mean values and hourly measured values (Fig. 3). Relative humidities calculated using a ratio of
3 match the measured values well, especially around the measured
maximums (i.e., 5 a.m.), while those calculated with the ratio of 2ﬁt the observed values well around the measured minimums (i.e., 5p.m.). Daily mean values clearly underestimate temporal cyclic
patterns of relative humidity. As can be observed in Fig. 3, relative
humidity values can reach 100%, which leads to calculations ofevaporation rates and heat ﬂuxes that are numerically unstable.
This is because having 100% relative humidity reverses the gradient in Eq. (27), even when the soil surface is relatively wet and predicted evaporation is equal to zero. However, in reality, theatmosphere is always mixing, and evaporation from the wet soil
surface rarely stops completely. To avoid such conditions, the maximum atmospheric relative humidity in the program was set to99.5% to allow soil water to evaporate even when the atmospherichumidity is near saturation.
Climatic variables
Net radiation calculations require values of the albedo, the soil
emissivity, and the atmospheric emissivity at any given time.
While albedo and soil emissivity can be computed using relatively
simple equations, there is a great number of available and accepted
formulas to calculate the atmospheric emissivity, as it is physicallymuch more complicated than other climatic variables. Fig. 4 showsthat different models provide quite different atmospheric emissivity values, especially for the low atmosphere humidity (i.e., 10%).
As water molecules in the air are primary energy emittancesources, the atmospheric emissivity values increase with the relative humidity. When the relative humidity is 10%, the maximum
difference among different formulas is about 0.4 (or 40%). Models
that depend on both the air temperature and the relative humidity
have, in general, similar trends; atmospheric emissivity values increase almost linearly as the air temperature increases. Slopes ofthese quasi-linear lines are quite similar. The two models that depend only on the air temperature provide much higher atmospheric emissivity values than the other formulas when therelative humidity is small (Idso and Jackson, Eq. (25) and Swinbank, Eq. (26)). As the relative humidity increases, the difference
between different models decreases. The relative humidity in arid

275

280

285

290
295
300
Temperature [K]

305

310

315

305

310

315

1
Atmospheric emissivity [-]

0

0.9
0.8
0.7
0.6
0.5
0.4

Relative humidity = 50%
0.3

275

280

285

290
295
300
Temperature [K]1Atmospheric emissivity [-]Relaitve humidity [%]

100

0.9
0.8
0.7
Brunt (1932)
Chung & Horton (1987)
Brutsaert (1975)
Brutsaert 2 (1975)
Idso (1981)
Idso 2 (1981)
Satterlund (1979)
Idso & Jackson (1969)
Swinbank (1969)

0.6
0.5
0.4

Relative humidity = 90%
0.3

275

280

285

290
295
300
Temperature [K]

305

310

315

Fig. 4. The atmospheric emissivity obtained using 10 different models as a functionof the air temperature (273–315 K) for three different relative humidity values: 10(top), 50 (middle), and 90 (bottom)%.

or semi-arid regions can be very low and thus the choice of theatmospheric emissivity model may have an important impact on
predictions of soil water contents and temperatures.
Soil temperature predictions
Different meteorological models mainly affect the various energy ﬂuxes (e.g., reﬂection, emissivity, etc.) and their distribution
during the day. Consequently, these models signiﬁcantly affectpredicted temperatures, while having only a limited affect on predicted pressure heads or water contents (especially at deeper
depths), which showed only very small differences between different runs. Since simulated pressure heads and water ﬂuxes at the
study site were already extensively evaluated by Scanlon et al.
(2002, 2005), in this section we will focus primarily on measured
and predicted temperatures at different depths.
Predicted and observed soil temperatures at four observations
nodes (ON1, 3, 5, and 7) are depicted in Fig. 5 for the Reference Scenario. Soil temperatures were predicted very well for the ﬁrst 250

554

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

oTemperature [ C]40

Reference30
20
ON1: Predicted
ON1: Observed
ON3: Predicted
ON3: Observed
ON5: Predicted
ON5: Observed
ON7: Predicted
ON7: Observed

10
0
0

100

Day

200

300

Fig. 5. Hourly temperature variations observed (open circles) at four observation nodes (ON1, 3, 5, and 7) during 365 days, along with those predicted using the Reference
Scenario.

days, with the largest deviations within few temperature degrees.
Despite some small discrepancies between the measured and predicted soil temperatures after 250 days, the overall agreement is rea-

sonably good, considering that no calibration was involved, only
default values of thermal properties were used, and only daily climatic information was used at input. Simulated soil temperatures

4

12
Reference
Wind
RH
Albedo
Ea
Es
Ra
Rs

ON1

10
8

MSE

MAE

3

2

Reference
Wind
RH
Albedo
Ea
Es
Ra
Rs

ON1

6
4

1
2
0

0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

Scenario

Scenario

3

8

ON3

ON3
6

MSE

MAE

2
4

1
2

0

0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

Scenario

Scenario
5

ON7

ON7

2

4

MSE

MAE

3

1

2

1

0

0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

Scenario

Scenario

Fig. 6. Mean absolute errors (MAE) for soil temperatures at observation nodes 1, 3,
and 7 for different scenarios.

Fig. 7. Mean square errors (MSE) for soil temperatures at observation nodes 1, 3,
and 7 for different scenarios.

555

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

1
Reference
Wind
RH
Albedo
Ea
Es
Ra
Rs

ON1
0.98

R

2

0.96

0.94

0.92

0.9
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

Scenario
Fig. 8. Correlation coefﬁcients (R2) between observed and predicted temperaturesat the observation node 1 (ON1) for different scenarios.

overestimate the observed values after 300 days, although discrepancies between observed and predicted temperatures are held within 5 °C, even at the 2-m depth, most of the time. The MAE at ON7 after
300 days increased to 3.54 °C, from the overall MAE of 1.43 °C. The

reason for this overestimation after 300 days is unknown, but it
may be associated with soil thermal properties that were not measured, but selected from the database based on the textural class.
Diurnal ﬂuctuations, as well as seasonal variations, are reproduced
well at each observation node. Seasonal variations, or in other words
annual ﬂuctuations, are more important in deeper layers where
there are almost no observed diurnal temperature changes. Thesevariations are also more critical for the long-term simulations, which
are important for the engineering cover performance assessment.
Mean absolute errors (MAE) and mean square errors (MSE) at
observation nodes 1, 3, and 7 (i.e., at depths of 15, 60, and200 cm, respectively) for all considered simulation scenarios areplotted in Figs. 6 and 7. Additionally, Fig. 8 shows the correlation
coefﬁcients between observed and predicted temperatures at allobservation nodes for all simulation scenarios. For most scenarios,
both MAE and MSE decrease from the upper node (ON1) to thelower node (ON7). Prediction errors decrease with depth, becausetemperature variations are damped in deeper layers, which lead tosmaller variability in prediction errors. There are, however, several
scenarios where MAE is greater at ON7 than at ON3. There is little
difference in the MAEs at ON7, probably because the overestima-

40

Temperature [oC]

35
30
25
20
15
10

7

8

9

10

11

12

13

14

44

45

46

47

Day

Temperature [oC]

20

15

10

5
40

41

42

43

Day
30

Temperature [oC]

25
20
15
10

Observed
Reference
Scenario 1
Scenario 2
Scenario 3
Scenario 4

5
0
181

182

183

184

185

186

187

188

Day
Fig. 9. Observed and predicted hourly temperature variations at ON1 for one week at three different periods (days 1–8, 51–58, 151–158). Different simulation scenarios(Reference, 1, 2, 3, and 4) are compared.

556

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

tion of the soil temperature after 300 days screens out small errors
before day 300. This may also explain increases in MAEs for somescenarios at ON7 when compared to ON5. Some unexpectedbehavior at ON7 may also be due to the boundary effect, sinceON7 is located on the boundary between layers 2 and 3, where
the soil type changes from loam to gravel. Overall, however, theranking of different scenarios based on errors is almost the sameat different depths, except for the MSE at ON7. In other words,
the propagation of errors at greater depths is independent of the
choice of meteorological models. As the choice of models affects
only the calculation of boundary conditions, errors induced at thesurface boundary for a particular scenario simply propagate to deeper layers without being magniﬁed or increased by errors fromother sources.
Prediction errors calculated using the daily average values of
the wind speed and relative humidity (Scenarios 2 & 4) are greaterthan those obtained with simple functions (Reference Scenario and
Scenarios 1 & 3). Although improvements in predictions are relatively small in the long-term simulation, as depicted in Figs. 6
and 7, the results suggest that simple functions (e.g., Eq. (6)) used
to generate these two climatic variables should be utilized when
only daily average values are available. Although improvements

in predictions seem to be much larger for short-term simulations
that focus primarily on short-term variations of temperatures, aslong as continuous values of air temperatures are obtained froma simple function (Eq. (2)), the improvements are relatively small(Fig. 9).
The impact of the choice of the albedo formula on temperature
predictions is relatively small, as reﬂected by both MAE and MSE at
all considered depths (Scenarios 5 & 6 in Figs. 6 and 7). It is obvious
from Eqs. (12) and (13) that calculated albedo values are almost
the same for different models when the soil surface is dry(h < 0.1), and such conditions were prevalent during most of thesimulation time at the study site. Most accepted albedo models
can thus be used for the bare soils of arid or semi-arid regions. Thisconclusion may, however, may not apply for humid conditions orfor soil proﬁles with vegetation.
When an appropriate atmospheric emissivity model is selected,
the prediction error can easily be reduced by 50% or more (compare Scenarios 7–14 with the Reference Scenario in Figs. 6 and
7). Among the various climatic variables considered in this study,
models computing the atmospheric emissivity produce the widest
range of values (Fig. 4). The Brunt equation (Eq. (18)), used in Scenario 7, leads to the smallest prediction errors and produces the

Latent heat flux [W m-2]

15

10

5

0
25

30

35

40

45

215

220

Day

Latent heat flux [Wm-2]

2

1.5

1

0.5

0
200

205

210

Day
80

Latent heat flux [Wm-2]

70

Reference
Scenario 16
Scenario 17

60
50
40
30
20
10
0
300

305

310

315

320

Day
Fig. 10. Calculated hourly latent heat ﬂuxes at the surface for 20-day time intervals at three different periods (days 25–45, 200–220, and 300–320). Different simulation
scenarios (Reference, 16, and 17) are compared.

557

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

smallest atmospheric emissivity values for the large interval of relative humidity and temperature values, except when both relative
humidity and temperature are high (Fig. 4). The higher the atmospheric emissivity, the more incoming longwave radiation fromthe atmosphere reaches the soil surface, and vice versa. There isno obvious explanation why the Brunt equation works better thanthe other formulas for this site. Although models that depend only
on the temperature (Eqs. (25) and (26)) should not necessarily be
excluded because of their relatively poor prediction performancein our study, those that depend either on the atmospheric pressureor the vapor density (Eqs. (18)–(24)) are likely to produce better
predictions.
Two models that calculate the soil emissivity (Eqs. (16) and(17)) were compared. As expected, because of their similarity,
there is no apparent difference in prediction errors between the
Reference Scenario and Scenario 15 in Figs. 6 and 7.
Estimation of evaporation rates requires that more parameters
be determined in comparison to other climatic variables. There isno distinct difference in prediction errors between the approximate approach (the Reference Scenario) and the iterative approach(Scenario 16) when the MO length is used in aerodynamic resistance calculations. This implies that there is no need to use theiterative approach, which demands more computational effort,
for the long-term simulations. Both MAE and MSE for the scenario
using the Richardson number (Scenario 17) are consistently smaller than those of the Reference Scenario, at all depths. Fig. 10 shows
the calculated latent heat ﬂuxes for three scenarios at the soil surface for three 20-days periods. Except for some minor ﬂuctuations
observed when latent heat ﬂuxes are relatively small (from day

200 to day 220), the three scenarios are comparable regardless ofthe magnitude of latent heat ﬂuxes. This demonstrates that the
MO stability correction is no more accurate than the simpler method. Because of its simplicity, the approach that uses the Richardson
number may be preferable to the one using the MO length. Thisconclusion should be valid not only for this particular dataset,but for other sites as well.
Among different models calculating the soil surface resistance,rs (Scenarios 18–21), the equation cited by Mahfouf and Noilhan(1991) outperforms the others with regard to prediction errors(Figs. 6 and 7). The empirical parameters used by Mahfouf andNoilhan (1991) also worked reasonably well for our study site.
Although different in details, all the formulas considered in thisstudy are in the exponential form, which partly explains relatively
small differences in prediction errors.
The correlation coefﬁcients, r2, between the observed and predicted temperatures, all very high (greater than 0.95), can be used
to assess agreement in the overall trend of model predictions. Even
with relatively large errors, r2 can be close to one if the generaltrend is predicted well. High r2s for all scenarios indicate that
although differences in meteorological models do not signiﬁcantly
affect the overall trend of predicted temperatures, they may shift
the predicted values away from observed temperatures. This is
conﬁrmed in Figs. 11–13. One can see that the predicted temperatures using different meteorological models are more-or-less parallel. As previously mentioned, the ranking of different methods
with depth remains almost the same as errors introduced at thesoil–atmosphere interface are propagated to deeper layers. It
seems that no additional errors are introduced as the simulations

Temperature [ oC]

20
15
10
5
0

ON1
50

75

100

125

150

225

250

Day

Temperature [ oC]

40
30

ON1

20
10
0
150

175

200

Day

Temperature [ oC]

45
40

ON1

4
7
17
Reference
ON1

35
30
25
20
250

275

300

325

350

Day
Fig. 11. Observed and predicted hourly temperature variations at ON1 from day 50 to day 365. Different simulation scenarios (Reference, 4, 7, and 17) are compared (notice
different scales in different ﬁgures).

558

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

proceed. Although Figs. 11–13 show results for only three observation nodes, similar trends were observed at all depths.
Our results indicate that using more complicated models does
not necessarily lead to better predictions of the soil temperature
with depth than using simpler models. This conclusion should beapplied to other sites as well. The majority of available meteorological models resulted in relatively good predictions of the soil temperature. However, we need to stress that the choice ofmeteorological models may lead to unwanted effects on temperature predictions, and thus needs to be done carefully. As importantas it is to exert efforts to parameterize the soil hydraulic and thermal properties, determining the proper formulas for the components required to solve the energy balance equation is alsoessential.
Summary and conclusions
A complete evaluation of the soil thermal regime can be obtained by simulating the coupled movement of liquid water, water
vapor, and heat energy in the subsurface. In our study, we modiﬁed
the HYDRUS-1D software package to increase its ﬂexibility inaccommodating various types of climatic information, and to solve
the coupled equations describing liquid water, water vapor, andheat transport in soils. Evaluation of the thermal regime of ﬁeld
soils requires, in addition to the solution of the subsurface waterﬂow and heat transport equations, the solution of equations
describing the surface water and energy balance. When only daily
climatic information is available, available meteorological models

need to be used to obtain continuous diurnal values of climatic
data when solving the energy balance equation. The impact of
the choice of such models on soil temperature predictions in thebare soil was investigated in this study. In addition, many different
models for components used in the energy balance equation werecompared. A variety of meteorological models were tested at the
proposed low-level radioactive-waste disposal site in Chihuahua
Desert in West Texas, where a thorough assessment of thermal
conditions needs to be conducted. Different meteorological models
were compared in terms of prediction errors of soil temperaturesat seven observation depths.
Models used to generate diurnal cycles of the wind speed andthe relative humidity were compared ﬁrst. As expected, using daily average values led to greater prediction errors, while relatively
simple extension methods could signiﬁcantly improve soil temperature predictions. Models needed to obtain the net radiation
were investigated next. While different models for the albedo
and the soil emissivity had little impact on soil temperature predictions, the choice of the atmospheric emissivity model had alarger impact. Some methods could decrease prediction errors
by 50% or more when compared to other methods. Soil and aerodynamic resistance models used in the calculation of the evaporation rate were tested last. As for the aerodynamic resistance, themodel based on the Richardson number outperformed that using
the Monin–Obukhov length. Comparable results were obtainedfor different soil resistance models, as all use similar exponential
forms. A comparison of all the different meteorological models
indicates that the error introduced at the soil–atmosphere inter-

Temperature [ oC]

20
15
10
5
0

ON2
50

75

100

125

150

225

250

Day

Temperature [ oC]

40
30

ON2

20
10
0
150

175

200

Day

Temperature [ oC]

45
40

ON2

4
7
17
Reference
ON2

35
30
25
20
250

275

300

325

350

Day
Fig. 12. Observed and predicted hourly temperature variations at ON2 from day 50 to day 365. Different simulation scenarios (Reference, 4, 7, and 17) are compared (notice
different scales in different ﬁgures).

559

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

Temperature [ oC]

30

20

10

0

ON5
50

75

100

125

150

225

250

325

350

Day

Temperature [ oC]

30

ON5
20

10

0
150

175

200

Day

Temperature [ oC]

30

ON5
25
4
7
17
Reference
ON5

20
250

275

300

Day
Fig. 13. Observed and predicted hourly temperature variations at ON5 from day 50 to day 365. Different simulation scenarios (Reference, 4, 7, and 17) are compared (notice
different scales in different ﬁgures).

face propagates to deeper layers. In other words, the ranking ofdifferent methods does not change much for different depths.
We found that there is a best combination of meteorologicalmodels for this particular study site. However, the objective of
our study was not to recommend the best combination of models
to predict soil temperatures. It was rather to show quantitatively
how the choice of meteorological models affects the soil temperature predictions. The results indicate that most of the meteorological models compared in this study can be used with conﬁdence tosolve the energy balance equation at the soil–atmosphere interface
in coupled water, vapor, and heat transport models. Improvements
in soil temperature predictions can be achieved not only by determining precisely the soil hydraulic and thermal properties, but also
by selecting proper meteorological models to solve the surface energy balance equation. Most conclusions should hold not only forthe site analyzed in this study, but for other sites where soil temperatures are to be predicted.

Acknowledgements
We would like to express our sincere gratitude to Bridget R.
Scanlon for providing us with her soil temperature data. This paper
is based on work supported in part by the BARD (Binational Agricultural Research and Development Fund) Project IS-3823-06 andby the National Science Foundation Biocomplexity programs#04-10055 and NSF DEB 04-21530. Financial support for the travel

costs of the second author by the Japan Society for the Promotionof Science is greatly appreciated.
Appendix
The governing equations of variably-saturated water ﬂow andheat transport in the vadose zone, as well as the soil hydraulicand thermal properties used in HYDRUS-1D are presented here.
Details can also be found in Saito et al. (2006).
Coupled liquid water, water vapor, and heat transport
The governing equation for one-dimensional ﬂow of liquidwater and water vapor in a variably saturated rigid porous medium
is given as follows:



@h
@
@h
@T
@h
@T
S
¼
K Lh
þ K Lh þ K LT
þ K vhþ K vT
@t @z
@z
@z
@z
@z


@
@h
@T
S
¼
K Th
þ K Lh þ K TT
@z
@z
@z

ð52Þ

where h is the total volumetric water content (m3 m3), h is the pressure head (m), T is the temperature (K), z is the spatial coordinate positive upward (m), and S is a sink term usually accounting for the root
water uptake (s1). KLh (m s1) and KLT (m2 K1 s1) are the (isothermal and thermal) hydraulic conductivities for liquid phase ﬂuxes due
to gradients in h and T, respectively, Kvh (m s1) and KvT (m2 K1 s1)

560

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561

are the isothermal and thermal vapor hydraulic conductivities,
respectively, and KTh (m s1) and KTT (m2 K1 s1) are the isothermal
and thermal total hydraulic conductivities, respectively.
The governing equation for the movement of energy in a variably saturated rigid porous medium, on the other hand, is given
as follows:@C p T
@hv
@
@T
@q T
@q
@q T
þ L0
kðhÞ
 C w L  L0 v  C v v  C w ST
¼
@t
@z
@z
@z
@z
@t
@z
ð53Þ
3

1

where Cw, Cv and Cp are volumetric heat capacities (J m K ) of theliquid water, water vapor, and moist soil (de Vries, 1963), respectively, L0 is the volumetric latent heat of vaporization of water(J m3), and the last term on the right side represents a sink of energy associated with root water uptake. k(h) is the apparent thermal
conductivity of the soil (J m1 s1 K11), and qL and qv are the ﬂux
densities of liquid water and water vapor (m s1), respectively.
Model parameters
The pore size distribution model of Mualem (1976) was used to
predict the isothermal unsaturated hydraulic conductivity function, KLh, from the saturated hydraulic conductivity Ks and thevan Genuchten’s (1980) model of the soil water retention curve:(hl ðhÞ ¼

r
hr þ ½1þjhsah
hjn m

h<0

hs

hP0

K Lh ¼ K S Sle ½1  ð1  S1=m
Þm 2
e

ð54Þ
ð55Þ

where hs and hr are the saturated and residual water contents(m3 m3), respectively, Ks is the saturated hydraulic conductivity(m s1), Se is the effective saturation (–), and a (m1), n (–), m(=1  1/n), and l (–) are empirical shape parameters. The parameterl was given a value of 0.5 as suggested by Mualem (1976).
The thermal hydraulic conductivity function, KLT, in Eq. (52) isdeﬁned as follows (e.g., Noborio et al., 1996b):



1 dc
K LT ¼ K Lh hGwT
c0 dT

ð56Þ

where GwT is the gain factor, which quantiﬁes temperature dependence of the soil water retention curve (Nimmo and Miller, 1986),
c is the surface tension of soil water (J m2), and c0 is the surface
tension at 25 °C (=71.89 g s2).
The isothermal, Kvh, and thermal, KvT, vapor hydraulic conductivities are described as (e.g., Nassar and Horton, 1989):

Mg
Hr
RT
D
dq
K vT ¼
gH v sqw r dT

K vh ¼

D

qw

qv s

ð57Þ
ð58Þ

where D is the vapor diffusivity in soil (m2 s1), qvs is the saturated
vapor density (kg m3), M is the molecular weight of water(M mol1) (=0.018015 kg mol1), g is the gravitational acceleration(m s2) (=9.81 m s2), R is the universal gas constant (J mol1 K1)(=8.314 J mol1 K1), g is the enhancement factor (–) (Cass et al.,1984), and Hr is the relative humidity (–) (=RH/100).
The apparent thermal conductivity of the soil, k(h), in Eq. (53)combines the thermal conductivity of the porous media in the absence of ﬂow, and the macrodispersivity which is assumed to be a
linear function of the velocity (de Marsily, 1987). The apparent
thermal conductivity, k(h), may then be expressed as (e.g., Šimůnek
and Suarez, 1993):

kðhÞ ¼ k0 ðhÞ þ bC w jqjl

ð59Þ

where b is the thermal dispersivity (m). The thermal conductivity,k0(h), accounts for the tortuosity of the porous medium, and canbe described with a simple equation given by Chung and Horton(1987):

k0 ðhÞ ¼ b1 þ b2 h þ b3 h0:5

ð60Þ

where b1, b2, and b3 are empirical regression parameters(W m1 K1). Chung and Horton (1987) also provided average values for the b coefﬁcients for three textural classes (i.e., clay, loam,
and sand), which are implemented in the HYDRUS-1D program(Šimůnek et al., 1998).
References
Aluwihare, S., Watanabe, K., 2003. Measurement of evaporation on bare soil andestimating surface resistance. Journal of Environmental Engineering 129,
1157–1168.
Beljaars, A.C.M., Holtslag, A.A.M., 1991. Flux parameterization over land surfaces foratmospheric models. Journal of Applied Meteorology 30, 327–341.
Boulet, G., Braud, I., Vauclin, M., 1997. Study of the mechanisms of evaporation
under arid conditions using a detailed model of the soil–atmosphere
continuum. Application to the EFEDA I experiment. Journal of Hydrology 193,
114–141.
Brunt, D., 1932. Notes on radiation in the atmosphere. Quarterly Journal of the Royal
Meteorological Society 58, 389–418.
Brutsaert, W., 1975. On a derivable formula for long-wave radiation from clear
skies. Water Resources Research 11, 742–744.
Brutsaert, W., 1982. Evaporation into the Atmosphere: Theory, History, andApplications. D. Reidel Publishing Company, Dordrecht, The Netherlands.
Camillo, P.J., Gurney, R.J., 1986. A resistance parameter for bare-soil evaporation
models. Soil Science 141, 95–105.
Camillo, P.J., Gurney, R.J., Schmugge, T.J., 1983. A soil and atmospheric boundary
layer model for evapotranspiration and soil moisture studies. Water Resources
Research 19, 371–380.
Campbell, G.S., 1977. An Introduction to Environmental Biophysics. Springer-Verlag,
New York.
Campbell, G.S., 1985. Soil Physics with Basic, Transport models for soil–plantsystems. Elsevier, New York.
Cass, A., Campbell, G.S., Jones, T.L., 1984. Enhancement of thermal water vapor
diffusion in soil. Soil Science Society of America Journal 48, 25–32.
Choudhury, B.J., Reginato, R.J., Idso, S.B., 1986. An analysis of infrared temperatureobservations over wheat and calculation of latent heat ﬂux. Agricultural andForest Meteorology 37, 75–88.
Chung, S.-O., Horton, R., 1987. Soil heat and water ﬂow with a partial surface mulch.
Water Resources Research 23, 2175–2186.
de Marsily, G., 1987. Quantitative Hydrogeology, Groundwater Hydrology forEngineers. Academic Press, Inc., Orlando, FL.
de Vries, D.A., 1963. The thermal properties of soils. In: van Wijk, R.W. (Ed.), Physics
of Plant Environment. North Holland, Amsterdam, pp. 210–235.
Ephrath, J.E., Goudriaan, J., Marani, A., 1996. Modelling diurnal patterns of air
temperature, radiation, wind speed and relative humidity by equations fromdaily characteristics. Agricultural Systems 51 (4), 377–393.
FAO, 1990. FAO Penman-Monteith Formula. In Report: Expert consultation onrevision of FAO methodologies for crop water requirement. May 28–31, Rome,
Italy.
Fayer, M.J., 2000. UNSAT-H version 3.0: unsaturated soil water and heat ﬂow model,
theory, user manual, and examples. Rep. 13249. Paciﬁc Northwest National
Laboratory, Richland, Washington.
Fuchs, M., Tanner, C.B., 1966. Infrared thermometry of vegetation. Agronomy
Journal 58, 597–601.
Gregory, J.M., 1989. Wind data generation for Great Plains locations. In: Proceedingsof American Society of Agricultural Engineers, International Winter Meeting,
December 12–15, New Orleans, LA.
Gregory, J.M., Peterson, R.E., Lee, J.A., Wilson, G.R., 1994. Modeling wind and relative
humidity effects on air quality. In: International Conferences on Air Pollution
from Agricultural Operations. Midwest Plan Service, Ames, IA, p. 183–190.
Hatﬁeld, J.L., Reginato, R.J., Idso, S.B., 1983. Comparison of long-wave radiation
calculation methods over the United States. Water Resources Research 19,
285–288.
Idso, S.B., 1981. A set of equations for full spectrum and 8- to 14-lm and 10.5- to
12.5-lm thermal radiation from cloudless skies. Water Resources Research 17,
295–304.
Idso, S.B., Jackson, R.D., 1969. Thermal radiation from the atmosphere. Journal ofGeophysical Research 74, 5397–5403.
Jury, W.A., Horton, R., 2003. Soil Physics, sixth ed. John Wiley & Sons, Inc., New York.
Kirkham, D., Powers, W.L., 1972. Advanced Soil Physics. John Wiley & Sons, NewYork.
Koivusalo, H., Heikinheimo, M., Karvonen, T., 2001. Test of a simple two-layer
parameterisation to simulate the energy balance and temperature of a snow
pack. Theoretical and Applied Climatology 70, 65–79.
Kondo, J., Saigusa, N., Sato, T., 1990. A parameterization of evaporation from baresoil surface. Journal of Applied Meteorology 29, 385–389.
Kondo, J., Saigusa, N., Sato, T., 1992. A model and experimental study of evaporationfrom bare-soil surfaces. Journal of Applied Meteorology 31, 304–312.

H. Saito, J. Šimůnek / Journal of Hydrology 373 (2009) 545–561
Leij, F.J., Alves, W.J., van Genuchten, M.T., Williams, J.R., 1996. The UNSODA
unsaturated soil hydraulic database, user’s manual, version 1.0, EPA/600/R-96/
095, Natl. Risk Manage. Lab., Off of Res. and Dev., Cincinnati, Ohio.
Mahfouf, J.F., Noilhan, J., 1991. Comparative study of various formulations ofevaporation from bare soil using in situ data. Journal of Applied Meteorology 30,
1354–1365.
Monteith, J.L., Unsworth, M.H., 1990. Principles of Environmental Physics, second
ed. Edward Arnold, London.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity ofunsaturated porous media. Water Resources Research 12, 513–521.
Nassar, I.N., Horton, R., 1989. Water transport in unsaturated nonisothermal salty
soil: II. Theoretical development. Soil Science Society of America Journal 53,
1330–1337.
Nimmo, J.R., Miller, E.E., 1986. The temperature dependence of isothermal moisturevs. potential characteristics of soils. Soil Science Society of America Journal 50,
1105–1113.
Noborio, K., McInnes, K.J., Heilman, J.L., 1996a. Two-dimensional model for water,
heat, and solute transport in furrow-irrigated soil: I. Theory. Soil Science Societyof America Journal 60, 1001–1009.
Noborio, K., McInnes, K.J., Heilman, J.L., 1996b. Two-dimensional model for water,
heat, and solute transport in furrow-irrigated soil: II. Field evaluation. Soil
Science Society of America Journal 60, 1010–1021.
Oke, T.R., 1978. Boundary Layer Climates. Methuen & Co Ltd., London.
Ortega-Farias, S., Antonioletti, R., Olioso, A., 2000. Net radiation model evaluation at
an hourly time step for Mediterranean conditions. Agronomie 20, 157–164.
Paulson, C.A., 1970. The mathematical representation of wind speed andtemperature proﬁles in the unstable atmospheric surface layer. Journal ofApplied Meteorology 9, 857–861.
Philip, J.R., de Vries, D.A., 1957. Moisture movement in porous materials undertemperature gradient. Transactions, American Geophysical Union 38, 222–
232.
Saito, H., Šimůnek, J., Mohanty, B.P., 2006. Numerical analysis of coupled water,
vapor, and heat transport in the vadose zone. Vadose Zone Journal 5, 784–800.

561

Satterlund, D.R., 1979. An improved equation for estimating long-wave radiationfrom the atmosphere. Water Resources Research 15, 1649–1650.
Scanlon, B.R., Christman, M., Reedy, R.C., Porro, I., Šimůnek, J., Flerchinger, G.N.,
2002. Intercode comparisons for simulating water balance of surﬁcial
sediments in semiarid regions. Water Resources Research 38 (12), 1323.
doi:10.1029/2001WR001233.
Scanlon, B.R., Reedy, R.C., Keese, K.E., Dwyer, S.F., 2005. Evaluation of
evapotranspiration covers for waste containment in arid and semiarid regionsin the southwestern USA. Vadose Zone Journal 4, 55–71.
Šimůnek, J., Suarez, D.L., 1993. The UNSATCHEM-2D code for simulating two
dimensional variably saturated water ﬂow, heat transport, carbon dioxidetransport, and solute transport with major ion equilibrium and kinetic
chemistry. Version 1.1. Research Report No. 128, US Salinity Laboratory,
USDA, ARS, Riverside, California, 218pp.
Šimůnek, J., Šejna, M., van Genuchten, M.T., 1998. The HYDRUS-1D software
package for simulating the one dimensional movement of water, heat, andmultiple solutes in variably-saturated media. Version 2.0. IGWMC-TPS-70.
International Groundwater Modeling Center, Colorado School of Mines, Golden,
CO.
Sun, S.F., 1982. Moisture and heat transport in a soil layer forced by atmosphericconditions. MS Thesis, Dept. of Civil Engineering, University of Connecticut, 77
pp.
Swinbank, W.C., 1963. Long-wave radiation from clear skies. Quarterly Journal of
the Royal Meteorological Society 89, 339–348.
van Bavel, C.H.M., Hillel, D.I., 1976. Calculating potential and actual evaporationfrom a bare soil surface by simulation of concurrent ﬂow of water and heat.
Agricultural Meteorology 17, 453–476.
van de Griend, A.A., Owe, M., 1994. Bare soil surface resistance to evaporation byvapor diffusion under semiarid conditions. Water Resources Research 30, 181–
188.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Science Society of America Journal 44,
892–898.