Agricultural Systems 72 (2002) 215–239
www.elsevier.com/locate/agsy

Error analysis of soil temperature
simulations using measured and estimatedhourly weather data with 2DSOIL
D.J. Timlina,*, Ya. Pachepskyb, B.A. Acocka,1,
J. Šimunekc, G. Flerchingerd, F. WhislereaUSDA-ARS Alternate Crops and Systems Laboratory, Bldg 007, Rm 116,
10300 Baltimore Ave, Beltsville, MD 20705, USAb
USDA-ARS Hydrology Laboratory, Rm, 10300 Baltimore Ave, Beltsville, MD 20705, USAc
USDA/ARS US Salinity Laboratory, 450 W Big Springs Rd, Riverside, CA 92507-4617, USAd
Northwest Watershed Research Center, 800 Park Blvd, Suite 105, Boise, ID 83712-7716, USAe
Department of Plant & Soil Sciences, Mississippi State University, PO Box 9555,
Mississippi State, MS 39762, USA
Received 6 September 2000; received in revised form 30 April 2001; accepted 29 May 2001

Abstract
Many crop simulation models use 1-h time steps for atmospheric, soil and plant processes but
often meteorological data are only available as daily summaries. The objective of this study was toinvestigate how errors in estimation of hourly values of solar radiation and air temperature aﬀect
errors in simulation of soil temperature using the model 2DSOIL. 2DSOIL is a two-dimensional
ﬁnite element model that simulates water ﬂow, chemical, water and solute uptake by plant roots,chemical equilibria processes, and gas and heat transport in soil. The standard deviations of thediﬀerences between hourly estimated air temperatures were about 2 C and 85 W m2 for solar
radiation. The mean diﬀerence in simulated and measured soil temperatures using measured hourly
weather data for all depths at both sites was less than 1 C and standard deviations were about
1–3 C indicating low bias. The range of errors was highest in the surface soils when the soil
was wetted after rainfall. Relative to simulated soil temperatures using measured hourly data,
simulated soil temperatures using estimated data were, on average over all depths, 2 C lower
and standard deviations ranged from 2 to 3 C. The errors were similar over all depths. Use of
estimated hourly air temperature and radiation generally resulted in underpredictions of soil
temperature by 2–3 C and increased error. Also maximum daily soil temperatures were
underestimated. Published by Elsevier Science Ltd.
Keywords: 2DSOIL; Soil temperature; Simulation; Down scaling; Weather data

* Corresponding author.
1
Present address: Adjunct Professor, Duke University Phytotron, Duke University, Durham,
NC27708, USA.
E-mail addresses: dtimlin@asrr.arsusda.gov (D.J. Timlin).
0308-521X/02/$ - see front matter Published by Elsevier Science Ltd.
PII: S0308-521X(01)00075-0

216

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

Nomenclature


"


s
"
’;
"s
"w
a
b"
Es
J
L
R
R00,n
R0,n
Rn
RNs
Rt
Ru
Ta
td
tdk
Tdk
tdn
Tdry
Tmax
tmaxhr
Tmin
Tmint
Twet
Ty
V
H
G
RNLS
Hrain

psychrometric constant, kPa ( C)1
solar declination
vapor pressure deﬁcit, kPa
soil moisture content for the grid nodes at the soil surface
humidity ratio
albedo of soil
actual water vapor pressure for day (assumed constant), kPa
local latitude saturation vapor pressure, kPa
saturation vapor pressure at wet bulb temperature, kPa
atmospheric transmission coeﬃcientslope of saturation water pressure curve, kPa ( C)1
potential evaporation rate from soil surface, cm day1
day of the year
latent heat of evaporation [2500.8–2.37 (Temperature)], J g1
actual radiation incident at earth’s surface, W m2
radiation incident at the top of the atmosphere at noon, W m2
potential radiation incident at earth’s surface at noon, W m2
actual radiation incident at earth’s surface at noon, W m2
net radiation on a bare soil surface, W m2
daily solar radiation integral W m2
net upward long wave radiation, W m2
air temperature,  C
daylength, hours
time of dusk, hours
air temperature at dusk,  Ctime of dawn, hours
dry bulb temperature1  C
maximum air temperature during the day,  Ctime of maximum air temperature measured from dawn, hours
minimum air temperature during the day,  C
minimum air temperature for the following day,  C
wet bulb temperature  C
air temperature at sunset on the previous day,  Cwind speed at 2 m height, km h1sensible heat ﬂux between the soil and air, J cm2 d1
net heat ﬂux into or out of the soil, J cm2 d1
re-radiation of heat from the soil to the air (J cm2 d1)
heat ﬂux into the soil with rain (J cm2 d1)

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

217

Soil temperature is an important environmental variable for plant growth anddevelopment. Soil temperatures aﬀect root growth (Koski et al., 1988), and biological processes such as nitrogen dynamics and pesticide transformation (Choi et al.,1988). Heat radiating from the surface soil impacts the crop canopy. Crop simulators must therefore be able to calculate soil temperatures accurately to be useful ascrop management tools.
Many crop simulation models use 1-h time steps for atmospheric, soil and plant
processes but often meteorological data are only available as daily summaries. Meteorological data that provide boundary conditions for soil temperature simulations typically include air temperature at 2 m, radiation, rainfall, and wind speed. Several
procedures were suggested to downscale daily weather data to hourly values (i.e.
Floyd and Braddock, 1984; Geng et al., 1985; Ephrath et al., 1996). These procedures were judged by the accuracy of estimating the weather variables per se.
Because of high nonlinearity and complexity of responses of biological systems toenvironmental variables, errors in input data can be magniﬁed through the simulation
progress. Such error propagation was observed in modeling studies with estimatedsoil hydraulic properties (Wösten and van Genuchten, 1988; Leenhardt, 1995).
Because of error propagation, the accuracy of the input estimation per se does nottell what errors can be expected in simulation results of interest. Wösten and van
Genuchten (1988) suggested that the prediction accuracy of the estimated soil parameters can be quantiﬁed in terms of the accuracy of hydrologic model output. Leenhardt (1995) noted that error propagation depends on the model in use, and therefore
the estimation procedure has to be evaluated in the context of a speciﬁc model. The
situation with using downscaled weather data in crop simulations bears a similarity
to using estimated soil parameters. The procedures to downscale weather data
should be judged by the accuracy of simulations carried out with a speciﬁc model.
A modular simulator of soil processes, 2DSOIL, has been developed from othermodels to be easily interfaced with crop models (Timlin et al., 1996). 2DSOIL is a
two-dimensional ﬁnite element model that simulates water ﬂow, chemical, water andsolute uptake by plant roots, chemical equilibria processes, and gas and heat transport in soil. 2DSOIL has also a meteorological module to simulate hourly energy
exchange at the soil–atmosphere interface that downscales daily summaries tohourly values. This meteorological module was adapted from the GLYCIM model(Acock and Trent, 1991) but never fully evaluated. The objective of this study was toevaluate this module in 2DSOIL using two criteria: a correspondence between
measured and simulated hourly meteorological variables, and a correspondence
between simulated and measured soil temperatures.

1. Materials and methods
1.1. Weather and soil temperature data
The weather and soil temperature data are from two experiments carried out onbare soil. The ﬁrst data set was from in Pullman, WA (Flerchinger and Saxton,

218

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

1989). The soil is a Palouse silt loam (ﬁne-silty mixed Mesic Pachic Ultic Haploxeroll). The soil temperature data were collected every 3 h at a conventionally tilled
plot. The temperatures were not averaged over that period but ‘‘point’’ measurements were output. The measurement depths were surface (to 0.01 m), 0.076, 0.152,
0.254, 0.381, 0.533, 0.686, 0.838, 1.067, 1.372, and 1.676 m. The data were collected
over the period from 1 March to 22 September 1987. The second data set was
obtained from the Plant Sciences Research Center at Mississippi State University(Khorsandi et al., 1997). The soil was a Marietta loam (ﬁne-loamy, siliceous, thermic
Fluvaquentic Eutrochrept). Only daily average soil temperatures were available atthis site although hourly soil temperatures were originally measured. Soil temperature was measured in ridge and furrow positions of a ridge tilled soil. The distance
between furrows was 90 cm and the furrow height was 15 cm. Soil temperatures
were measured to 60 cm depth and measurements were available at 0.05, 0.1, 0.25,
and 0.6 m. The data are from the period 20 May to 27 September 1991. Soil
hydraulic properties and hourly values of air temperature, solar radiation, and windvelocity were available at both sites.

2. The meteorological model
The model of soil–atmosphere interaction uses relationships of celestial geometry,crop canopy geometry, and radiation balance theory (Acock and Trent, 1991; Acocket al., 2000). The model here is described as a sequence of submodels. A descriptionof the symbols and their units are listed in the Nomenclature. The time unit ofvariables that have a time component, for example evaporation rate, is day (d1).
These are instantaneous rates of the variables, and these rates are updated on anhourly basis which is the time step of the meteorological module of 2DSOIL.
2.1. Incident radiation submodel.
Solar declination, , is calculated using the algorithm of Robertson and Russelo (1968)
 ¼ þ
1i¼5
X½i sinð0:01721ði  1ÞJÞ þ iþ4 cosð0:01721ði  1ÞJÞ

ð1Þ

i¼2

Here and below index i indicates that the value is valid over the hour centered at
ti=(i0.5) h and J is day of year. The coeﬃcients i are equal to 1= 0.3964;
2=3.631; 3=0.03838; 4=0.07659; 5=0.0; 6=22.97; 7=0.3885;
8=0.1587; 9=0.01021.
Daylength, td is deﬁned as the interval between sunrise and sunset; at both of thesetimes the center of the sun’s disk is 50% below the horizon. The equation relating
solar altitude and declination to latitude and hour (Smithsonian Meteorological
Tables, 1966, p. 495) is solved for the hour angle and then results are converted tohours from solar noon. Doubling the answer gives day length:

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239


0:014544 þ sin’sin
td ¼ arccos 24
cos’cos

219



ð2Þ

with ’, latitude, in degrees.
The equations for solar radiation incident on top of the earth’s atmosphere, R00,n,
are given in the Smithsonian Meteorological Tables (1966, p. 417). The conditional
solar constant, i.e. the radiation ﬂux on the top of the troposphere is assumed to be1325.4 W m2 (Budyko, 1974):
R00;n ¼ 1325:4

cosð’  ÞZ2

ð3Þ

where Z=1+0.01674 sin[0.01721(J93.5)], is the radius-vector of the earth.
The atmospheric transmission coeﬃcient, a, is estimated for diﬀerent latitudes andtimes of year from the data of Budyko (1974):
a ¼ 0:68 þ ð1:57  0:1Þ

145  J;1000

J  237;
1000
1454J < 237

a ¼ 0:68 þ ½ ð2:04  0:37 Þ  0:19
a ¼ 0:68

J4145
J5237
ð4Þ

where =’/30.
Estimated potential radiation incident on the earth’s surface at noon R0,n, (Wm2) is calculated as


1
0:02
1=cosð’ÞþaR0;n ¼ R00;n 0:93 
ð5Þ
2
cosð’  ÞThis equation was derived from an equation given in the Smithsonian Meteorological Tables (1966, p. 420). The adsorption of the extraterrestrial radiation isallowed to vary with atmospheric path length, i.e. solar altitude, according to thedata given by Miller (1981).
Estimated actual radiation incident at the earth’s surface at noon, Rn, (W m2)
under a cloudless sky is calculated from the total measured daily radiation, Rt, (Wm2) and the assumption that the radiation ﬂux density varies as a half sine wave
over the photoperiod:
Rn ¼

Rt
td 2

ð6Þ

Both time of dawn, tdn, and time of dusk, tdk, are derived from daylength, td:
tdn ¼ 12  td =2;

tds ¼ 12 þ td =2

ð7Þ

220

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

Actual radiation incident at the earth’s surface at time ti, Ri, (W m2) is calculated
according to the half sine wave pattern:ði  0:5Þ
Ri ¼ Rn sin
ð8Þ
td
Corrections are made for incomplete hours after dawn and before dusk to obtainhourly values.
2.2. Temperature–vapor pressure submodel
The time of maximum air temperature, measured from dawn (tdn), is:



td
Tmax
 arcsintmaxhr ¼
Rn ð0:0945  8:06 105 Rn þ 6:77 1 4 Tmax Þð9Þ

where the expression after ‘arcsin’ is less than one. The coeﬃcients were ﬁt to datafrom Mississippi (Acock and Trent, 1991). The air temperature at sunset is:Tmax  Tmintd
3
1 þ sinþ
ð10Þþ TminTdk ¼
2
2
tmaxhr
Diurnal temperature variation is approximated by a half sine wave during the day
and by either a logarithmic or a linear dependency at night. At time t=ti we have:Tmax  Tminðti  tdn Þ 3
1 þ sinþ
; tdn 4ti 4tdk ;Ta;i ¼ Tmin þ
tmaxhr22"#tdn ti
Ty  Tmin 2tdn
1þ
1 ; ti < tdn ; Ty > Tmin ;Ta;i ¼ Tmin þ TT
tdn  ti
Ta;i ¼ Tmin þ Ty  Tmin
; ti < tdn ; Ty 4Tmin ;
2tdn
#

"
 ti tdk
Tdk  Tmint
Tdk  Tmint 2ð24tdk ÞTa;i ¼ Tdk þ T 1 þ
1 ;
1þ
T
T
t > tdk ; Tmint < Tdk ;Ta;i ¼ Tdk þ ðTmint  Tdk Þ

ti  tdk; ti tdk ; Tmint 5Tdk
2ð24  tdk Þ

ð11Þ

These equations have given a good ﬁt to data sets from Arizona and Mississippi,
when the parameter T* was set at 5 C (Acock and Trent, 1991).
Calculations of water vapor pressure for the day depend on the availability of wet
and dry bulb temperature values. If those are not available, dew point temperature isassumed to be the minimum temperature. When wet and dry bulb temperatures areD.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

221

not known, the algorithm of Weiss (1977) is used to calculate water vapor pressure,", where Tmin is used for the dew point temperature:
17:27min" ¼ 0:61exp
ð12Þ
Tmin þ 237:3
The psychometric ‘‘constant’’,  is 0.0645 kPa  C1.
If wet and dry bulb temperatures are known, then calculations for the water vapor
pressure, ", include the value of the saturation vapor pressure at wet bulb temperature "w, humidity ratio , and latent heat of evaporation L (J g1):

17:27Twet
0:622"w;"w ¼ 0:61exp; ¼
Twet þ 237:3
101:3  "w
1:006 þ 1:846
L ¼ 2500:8  2:37Twet :  ¼ 62:81;
Lð0:622 þ Þ2" ¼ "w   Tdry  Twetð13ÞThe saturation vapor pressure, "s is calculated using air temperature (Ta). The
slope, b", of the saturation water vapor pressure curve is estimated using the value ofthe saturated vapor pressure, and ﬁnally, the vapor pressure deﬁcit, ", is then calculated from:
17:27Ta"s ¼ 0:61exp
Ta þ 237:3

17:27ðTa þ 1Þb" ¼ 0:61exp
 "sðTa þ 1Þ þ 237:3
" ¼ "s  "ð14ÞWhere " is from Eq. (11) or Eq. (12) depending on the availability of wet and dry
bulb temperatures.
2.3. Potential evapotranspiration submodel
The evapotranspiration submodel is based on the Penman (1963) equation. Net
upward long-wave radiation, Ru, is calculated using an approximation derived by
Linacre (1968). Linacre’s original approximation used the ratio of the mean numberof hours of bright sunshine to day length. This ratio has been replaced by the ratioof actual radiation at the earth’s surface at solar noon to the potential radiation.
Longwave radiation is calculated for each period of the day, ti, where radiation isgreater than zero using air temperature (Ta) at that time:Rn
Ru ¼ 11:2 ð100  Ta Þ
ð15Þ
R0;n

222

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

The albedo () of exposed soil is estimated from the data of Bowers and Hanks(1965) as a function of soil water content of the surface soil ():
 ¼ 0:3  0:5ð16Þ

Net radiation on the soil surface (RNs) for the time period, ti, is calculated fromalbedo (), net upward longwave radiation (Ru) and estimated actual radiation
incident at the earth’s surface [Eq. (8)]:RNs ¼ ð1  ÞRi  Ru

ð17Þ

Penman’s equation gives the potential evaporation rate from the soil, Es (g cm2
day1), for time period, ti:
24
Es ¼

b" 3600
RNs þ 109:375ð1: þ 0:419VÞ"
 L
b"
þ1


ð18Þ

The actual evapotranspiration will be determined by 2DSOIL using an iterative
process that depends on the current water content and transport properties of thesoil.
2.4. Boundary condition for the heat transport in soil
The heat balance at the soil surface is:G ¼ L Es þ RNs þ H

ð19ÞHere, L is the latent heat of vaporization of water (J g1), Es is evaporation rate (gcm2 day1), RNs is the radiant energy ﬂux (J cm2 day1 ), H is sensible heat ﬂux (Jcm2 day1 ) between the soil and air and G is the net heat ﬂux into or out of the soil
(J cm2 day1 ) . Units for heat are given as Joules in this section since these are theunits used in the soil module.
Sensible heat ﬂux (H) is calculated as a function of the temperature diﬀerence
between the soil surface and the air. For bare soil, the sensible heat transfer between
the soil and atmosphere occurs when a temperature diﬀerence exists between the soil
surface and air and the soil surface is warmer than the air:
H ¼ ð24:1 þ 8:38VÞ ðTa  Tsurf Þð20Þ

Here the expression for the coeﬃcient of heat conduction between the soil and theatmosphere is from Linacre, (1968, Appendix 2), V is windspeed at 2 m (km day1),
Ta is air temperature and Tsurf is soil surface temperature. When the temperature of
the soil is greater than the temperature of the air, some heat is re-radiated as long

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

223

wave radiation into the atmosphere, RNLS (Linacre, 1968). This re-radiation
becomes a component of H in Eq. (19). Re-radiation of heat from the soil to the air(Linacre, 1968, uses cal cm2 min1) is calculated as (Linacre, 1968):
RNLS ¼ ð0:00683 þ 0:00009 Ta Þ ðTa  Tsurf Þ Tsurf > Ta

ð21Þ

The units of RNLS are then converted to (J cm2 day1 ).
Heat inﬂux, Hrain (J cm2 day1 ) due to rainfall is calculated based on the intensity of rainfall and temperature of the rainwater (the latter assumed to be equal tothe air temperature):
Hrain ¼ IR CW ðTa  Tsurf Þð22Þ

where IR is rain intensity, g cm2 day1, and CW is the speciﬁc heat capacity ofwater, J g1 ( C)1.

3. Boundary conditions, and water ﬂow and heat transport parameters in 2DSOIL
2DSOIL is a ﬁnite element representation of matter and energy transport in soil.
The model uses adaptive time steps (Timlin et al., 1998). The minimum time step
depends on either convergence criteria for the iterative equations or the amount oftime that needs to pass before a boundary condition is updated. The default units ofthe time step can be any time value but all the parameters and coeﬃcients must
use the same unit. For the simulations in this paper, a time unit of day is used.
Radiation or potential evapotranspiration are updated each hour or 1/24 of aday. The shorter the time period between updates, the more accurate will be theinformation passed to the model, however, more detailed information would be
required.
The hydrology and heat transport components of 2DSOIL require, as input, soil
hydraulic parameters, the saturated soil hydraulic conductivity, Ksat (cm day1),moisture retention parameters, saturated water content, soil texture, and initialwater contents and temperatures. The soil hydraulic properties are described by thevan Genuchten equation (van Genuchten, 1980), calculations of thermal conductivity are from DeVries (1963), water ﬂow and heat transport are calculatedusing a two-dimensional ﬁnite element representation of the Richards equation(Šimůnek et al., 1994). All parameters for the model were chosen from the literatureor calculated from measured data.
The soil hydraulic properties for the Washington site included the saturated
hydraulic conductivity and Brooks-Corey parameters for the matric potential vs.
water content relationship. At the Mississippi site, saturated hydraulic conductivity
and measured values of matric potential and water content were available. van
Genuchten’s equation (van Genuchten, 1980) for moisture release and unsaturated
hydraulic conductivity was ﬁtted to the data. Parameters of the equation, along withsoil texture and saturated hydraulic conductivity are given in Table 1.

224

Table 1
Soil properties used for the Washington and Mississippi sitesa
Ksatsb (cm3 cm3) rcd (cm3 cm3) d (cm1) nd’;b (Mg m3) Sand and silt (kg kg1) Clay (kg kg1) Organic matter (kg kg1)

Washington
0.000–0.076
0.077–0.152
0.153–0.254
0.255–0.381
0.382–0.533
0.534–0.686
0.687–0.838
0.839–1.067
1.068–1.372
1.373–1.676

2.88
2.88
2.88
3.84
6.72
5.04
3.36
2.30
0.27
0.07

0.588
0.481
0.517
0.465
0.508
0.478
0.440
0.404
0.367
0.338

0.095
0.078
0.067
0.074
0.072
0.070
0.065
0.061
0.069
0.074

0.01290.01290.02570.03020.04270.03390.0271
0.02190.01540.0109

1.440
1.440
1.399
1.336
1.330
1.353
1.362
1.388
1.399
1.389

0.95
0.95
1.26
1.17
1.32
1.21
1.29
1.39
1.49
1.59

68.1
68.2
68.5
68.8
69.0
69.0
72.0
72.0
76.0
76.0

0.301
0.302
0.305
0.308
0.310
0.310
0.280
0.280
0.240
0.240

0.018
0.015
0.010
0.005
0.001
0.001
0.001
0.001
0.001
0.001

Mississippi
0.000–0.200 12.64 0.540
0.201–0.900 15.11 0.442
0.901–2.050 1.59 0.435

0.167
0.201
0.226

0.28940.37140.4442

1.326 1.40
1.237 1.60
1.226 1.40

77.0
69.0
68.0

0.210
0.300
0.315

0.020
0.010
0.005

a

The values s, r, , and n are parameters for the van Genuchten equation (van Genuchten, 1980), and ’b is bulk density.
Saturated water content.
cResidual water content.
dParameters for van Genuchten’s equation.
b

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

Depth (m)D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

225

The bottom boundary for water ﬂow was given as a free drainage condition, i.e.
ﬂux out the bottom boundary was equal to the unsaturated hydraulic conductivity
at the current matric potential (unit gradient). In absence of other information, aconstant boundary value of 20 C was used for temperature at the Mississippi site. A
time varying boundary condition was used at the Washington site. This boundary
condition was used at the Washington site since the temperatures at the 1.67 m
depth showed more variation over time than did the temperatures at the Mississippi
site. The temperature at the bottom boundary was set to the temperatures at the 1.67
m depth. Water and heat ﬂow are not coupled in 2DSOIL, and there is no provision
for vapor ﬂux.

4. Design and sequence of the simulations
Hourly values of temperature, radiation and wind velocity were calculated fromdaily weather data (total radiation, maximum and minimum temperatures, and daily
windrun) to evaluate the model used to estimate hourly meteorological data fromdaily summary data. In the second step, simulations of soil temperature were runusing measured hourly weather data as input. The measured soil temperature data
are point (instantaneous) measurements at Washington and averages of temperature
over 24 h at Mississippi. Instantaneous temperatures were output for the Washington data and daily averages for the Mississippi data. These averages were comparedto the measured temperatures. In order to evaluate errors of simulated soil temperatures, soil temperatures were also simulated using estimated hourly meteorological data. These simulations were carried out the same as the simulations using
measured hourly data. The comparison of simulation and measurements was carriedout for a long time sequence of measured values, 195 days for Washington and 100
days in Mississippi.
A bare, ﬂat soil surface was simulated for the Washington site. A bare, ridge andfurrow surface conﬁguration was used for the Mississippi site. No calibrations of
2DSOIL were made prior to the temperature simulations.
Relative error was used to compare errors in simulations of soil temperaturesusing measured and estimated hourly weather data. Relative error was deﬁned as thediﬀerence in soil temperatures simulated using measured hourly data, and soil temperatures simulated using estimated hourly values divided by soil temperatures
simulated using measured hourly data.

5. Results and discussion
5.1. Comparison of measured and calculated hourly air temperatures and solar
radiation
Measured hourly radiation and air temperature, and hourly values calculated
from daily weather data for both Washington (Fig. 1) and Mississippi (Fig. 2) have

226

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

similar diurnal distributions throughout the 24-h period. Measured solar radiation
for the ﬁve day period shown in Figs. 1 and 2 appears higher than predicted at
midday for both sites. The diﬀerences in radiation as a function of time of day weresimilar for both locations (Table 2). The diﬀerences seen at the peak radiation atnoon are noticeable because the measured radiations are somewhat larger than thecalculated radiations during most of the daylit hours. The mean diﬀerences in solar
radiation (measured–predicted) for all the data given in Table 2 were consistently
positive from about 11:00 until 18:00 h and increased to about +90 W m2 at
midday (Table 2). The diﬀerences were consistently negative before 11:00 h. The
positive diﬀerences, however are largely within the standard deviations of the mean
diﬀerences. The consistent positive diﬀerences do suggest a slight bias. The magnitudes of the predicted and measured temperatures at a particular time could be

Fig. 1. Measured and calculated hourly air temperature, and radiation for the Washington site for aselected period of 7 days.

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

227

diﬀerent by as much as 5–10 C and radiation by as much as 500 W m2 (Fig. 3)
over all the data.
The calculated air temperatures in the morning hours at both sites were, on average, lower than measured temperatures (Table 2). The mean bias in air temperatures
was greater for the Washington site than for the Mississippi site. The standard
deviations were comparable for the two sites, however and the diﬀerences for mosttimes of day were largely within the range of the standard deviations. The diﬀerences
between the measured and calculated hourly temperatures for the Mississippi site
between the hours of 13:00 to 14:00 were larger than the standard deviation. The
data in Fig. 2 for a 5-day period suggest that the actual temperatures were rising at a
slower rate than the calculated temperatures so that the calculated late-morning,
early afternoon air temperatures became higher than measured. The reason for thisFig. 2. Measured and calculated hourly air temperature, and radiation for the Mississippi site for aselected period of 7 days.

228

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

leveling oﬀ of temperature rise can be attributed to more eﬃcient mixing of heated
air from the ground level into the atmospheric boundary layer (Driedonks, 1981;
Ephrath et al., 1996). This can aﬀect the time of the maximum temperatures. Eprath
et al. (1996) reported that a sine-exponential model did not accurately predict thediurnal variation in air temperature. This was because maximum air temperatureoften occurred later than solar noon. They used a site speciﬁc parameter to lag
temperature with respect to time. Our approach was similar here [Eq. (9)] withrespect to calculating the time of maximum temperature.
Eq. (9) did give reasonable average estimates of the time of maximum temperatures for both sites. The observed times of maximum temperatures were highlyvariable, however and were probably aﬀected by meteorogical events such as
weather fronts (Eprath et al., 1996). This equation was developed for application of
the GLYCIM model to conditions in Mississippi, USA (Acock and Trent, 1991) but
Table 2Means and standard deviations (S.D.) of diﬀerences between measured and calculated hourly temperatures and radiation for the Washington (n=196 per hour) and Mississippi (n=128 per) sites as a functionof hour of daya
Site

Hour of day

Temperature ( C)

Radiation (W m2 )

Mean

S.E.

Mean

S.D.

7
8
9
10
11
12
13
14
15
16
17
18

2.3
2.6
2.5
2.0
1.3
0.6
0.0
0.4
0.6
0.7
0.5
0.7

1.88
2.10
2.03
1.90
1.68
1.34
1.20
1.22
1.39
1.71
2.07
2.34

170.1
148.5
101.5
46.6
7.0
51.1
89.2
95.8
90.6
96.2
94.8
92.9

69.84
81.78
87.31
99.51
93.27
97.37
103.8095.51
107.54
99.57
81.79
69.51

7
8
9
10
11
12
13
14
15
16
17
18

0.5
0.7
0.7
0.3
0.3
0.9
1.4
1.5
1.1
0.6
0.1
0.6

0.80
1.11
1.37
1.36
1.23
0.91
0.85
1.23
1.20
1.65
2.25
2.16

187.5
153.4
83.9
10.8
24.4
77.0
83.3
88.6
78.9
85.1
78.3
76.7

60.47
85.88
90.52
111.83128.51
126.47144.71133.51135.14125.21
97.38
55.34

Washington

MississippiaCalculated values are obtained from Eqs. (1)–(14).

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

229

was never evaluated for other conditions. A replacement for this equation could beincorporated into 2DSOIL using the method of Eprath et al. (1996). The inability tosimulate the percentage and timing of cloud cover will also add to the error. Cloud cover
will reduce the daily radiation integral however which allows some compensation.
Fig. 3 shows the distributions of the bias errors (measured- calculated) for the twosites. The ranges in error were similar for the two sites for either air temperature or

Fig. 3. Cumulative probabilities distributions of errors (bias) in calculated air temperatures and radiation
for the two sites. Absolute error (bias error) is deﬁned as measured-calculated temperature or radiation.

230

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

radiation. The error bias is larger for radiation than for air temperatures, however.
For all times of the day, 35% of the errors are positive for the Washington site and
53% for the Mississippi site. The errors in radiation are more evenly distributed
about zero than for temperature.

6. Soil temperature simulations using measured hourly air temperatures andradiation
Simulated soil temperatures from the Washington site capture the measured relative changes in soil temperature with time and depth (Fig. 4). The diﬀerences
between simulated and measured soil temperatures were generally larger during thedaylight hours than during the evenings. During several time periods, however, themeasured temperatures in the surface soil were as much as 15–20 C larger than thesimulated temperatures (Fig. 4b). The largest diﬀerences occurred during or immediately after days where there was rainfall. Similar diﬀerences in temperature areseen for deeper layers of the proﬁle.
The mean diﬀerences between measured and simulated soil temperatures for all
depths at the Washington site were all less than the standard deviation (Table 3).
The standard deviation for the mean diﬀerence was largest at the surface but similar
for the other depths. Overall there is no signiﬁcant bias in the simulated temperatures
and the standard deviation is in the range of 1–2 C. The larger range of temperaturesin the soil surface contributes to the increased error there.
There is a diurnal pattern to the errors for soil temperature simulations at the
Washington site as a function of DOY for the surface layer and 0.15 m depth(Fig. 5). The diurnal error indicates that the simulated soil temperatures were lower
than measured temperatures during the day and slightly higher during the night.

Table 3
Summary statistics for diﬀerences between measured and simulated soil temperatures ( C) at 10 depths at
the Washington site (n=1523 per depth)a
Depth

Mean diﬀerence

Maximum

Minimum

S.D.

R2

0
0.076
0.152
0.254
0.381
0.533
0.686
0.838
1.067
1.372
1.676

0.22
0.29
0.51
0.53
0.28
0.26
0.04
0.21
0.48
0.54
0.54

19.8
8.4
6.3
6.1
5.4
5.1
4.5
3.6
2.9
1.9
1

11.9
6.3
5.8
4.8
4.4
3.2
3.4
3.3
3.4
2.8
2.2

5.29
1.90
1.89
1.77
1.55
1.25
1.20
1.12
1.00
0.81
0.63

79
95
94
93
94
96
95
95
95
96
98

a

S.D. is standard deviation of the diﬀerences.

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

231

Fig. 4. Simulated and measured soil temperatures for the surface 0.01 and 0.15 m depths at Washington
for 4-time periods: (a) DOY 88–105, (b) DOY 118–130, (c) DOY 156–168, and (d) DOY 218–230. Soil
temperatures were simulated by 2DSOIL using measured hourly weather data.

232

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

These data represent a period with the largest errors for the entire period of simulations. One contributing factor to the errors may be the diﬀerence in the location of
the measurement instrument and the depth the simulation results were recorded.
There are large gradients in temperature at the soil surface and small diﬀerences insensor placement can result in large diﬀerences in measured temperature. The ﬁrst
layer in the simulation extended from the surface to 0.02 m. Compared to a sensor in
the ﬁrst cm of soil, the simulated values in the 0–0.02 m depth may have been more
damped during the hottest part of the day when the soil was very wet and lost heatfaster in the evening. Also, the surface soil 0.01 m was quite variable in terms of
meteorological conditions.

Fig. 5. Diﬀerences between measured and simulated temperatures at the Washington site as a function oftime for a selected period of 13 days for the soil surface and 0.15 m depth. Soil temperatures were simulated by 2DSOIL using measured hourly weather data.

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

233

The weather measurement interval and method may also contribute to error. The
soil temperature measurements at the Washington site represented instantaneous
values while the input weather data was averaged over 1 h. The averaging of weather
data results in the loss of information regarding the peak temperature and radiation
occurring over the time of measurement. The high measured soil temperatures may
have been in response to a very short period of high solar radiation.

Fig. 6. Simulated and measured soil temperatures at Mississippi site for the ridge and furrow zones. Soil
temperatures were simulated by 2DSOIL using measured hourly weather data.

234

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

Simulations by 2DSOIL gave reasonable results for simulated temperatures forthe ridges and furrows for Mississippi (Fig. 6). There is minimal bias in the simulated temperatures for the Mississippi site, the mean diﬀerences are all less than 1 C
and all less than the standard deviation (Table 4). There are no large diﬀerencesamong the errors for the four depths. The standard deviations, however were larger
for the Mississippi data than for the Washington data. This reﬂects the smaller
number of data for the Mississippi site. The ranges of diﬀerences for the two sites
are similar (Tables 3 and 4). There was less error for the 0.05 m depth at Mississippi
than for the 0.01 m depth at Washington.
Some errors in simulated temperature values may occur when the evaporation
process begins to slow as the surface soil dries. This is due to the method 2DSOIL
uses to determine the boundary condition when simulating surface soil drying.
Evaporation from a wet soil will initially occur as a constant ﬂux condition using theevaporation rate as a prescribed ﬂux. If, during an iteration, the surface node cannot
support the prescribed ﬂux, the boundary condition is changed to a constant head.
The value of the soil matric potential at this constant head is an input datum andapproximates the matric potential at air dry water content. When the model
switches to a constant head boundary condition, the actual rate of water loss fromthe surface cell is less than the potential rate given as the evaporation rate. Energy
not used to evaporate water then goes to heat the soil. The true matric potential atwhich the evaporation becomes ﬂux limited is diﬃcult to determine from the available data. The evaporation rate is a function of the radiative demand as well as
surface water content. It is more likely that errors will be larger between simulated
and measured values when the evaporation process is in this transition phase. Oncethe soil dries the simulated and measured temperatures became more similar.

Table 4
Summary statistics for diﬀerences between measured and simulated soil temperatures at the Mississippi
site (n=19 per depth)a
Row position

Depth (m)Mean diﬀerence ( C)

Maximum

Minimum

S.D.

0.05
0.10
0.25
0.60

0.8
0.8
0.7
0.7

5.2
5.1
3.9
2.7

4.0
3.7
2.3
0.7

2.68
2.48
1.87
0.93

0.05
0.10
0.25
0.60

0.7
0.4
0.5
0.2

5.1
4.5
3.5
2.6

3.7
3.2
2.2
1.4

2.67
2.49
1.93
1.10

Furrow

Ridge

aS.D. is standard deviation of the diﬀerences.

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

235

7. Simulated soil temperatures using calculated hourly air temperature and radiation
Use of hourly air temperatures and radiation that are calculated from max-min air
temperatures and daily radiation integrals can potentially introduce additional
errors into temperature simulations. Fig. 7 shows the cumulative probability ofrelative error for the Washington soil temperatures. The relative error ranges from200 to 64% although only 35 out of 1500 values have relative errors more than theabsolute value of 50%. More than 70% of the relative errors are greater than zero
which means the simulated soil temperatures using calculated data are largely
underestimated. Generally, soil temperatures in the higher range are underestimatedfor simulations with calculated hourly weather data. This is probably related to theunderestimation of hourly radiation.
Table 5 lists the average diﬀerences between measured and simulated soil temperatures and standard deviations for the Washington data by depth. The mean
diﬀerences as a function of depth were not large and varied from 1.7 to 2.5 C .The
diﬀerences also increased slightly with depth. The data indicate a slight bias in thesimulated temperatures when calculated hourly weather data are used because soil
temperatures were underestimated by 1–2 C. The standard deviations for the mean
diﬀerences between temperatures simulated using measured and calculated hourly

Fig. 7. Cumulative probability distribution of relative errors in simulated soil temperatures at Washington using measured and calculated hourly weather data for the surface 0.01, 0.15 and 0.53 m depths.
Relative error is calculated as (simulated with measured hourly–simulated with calculated hourly)/simulated with measured hourly.

236

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

weather data (Table 5) are higher than the errors for the mean diﬀerences of themeasured and simulated temperatures where measured hourly weather data wereused (Table 3).
The cumulative probabilities of relative error for the Mississippi data are shown in
Fig. 8. Only the probabilities for the ridge position are shown; there were not large
diﬀerences between the cumulative distributions for the ridge and furrow positions.
The relative errors range from 7.5 to 8% and more than 80% of the errors are
positive. This is similar to the Washington simulated results. The relative errors are
smaller for Mississippi than for the Washington site possibly because the soil temperatures from Mississippi were daily averages and not hourly values.
Soil temperatures simulated using calculated hourly weather data were similar tosoil temperatures simulated using measured hourly data except at the higher rangesof soil temperatures. This is partially due to the fact that the early morning andmaximum hourly value of radiation near noon are often underestimated. There wasmore error for surface temperatures than subsurface soil temperatures when calculated hourly data were used as input. This is probably because of the large impact ofradiation on surface temperature and evaporation. The impact of radiation isdamped for subsurface layers.
Heat energy is also added to the soil with rainfall. When daily weather data areused in models, rainfall data is used as a daily summation, and in 2DSOIL, rainfall
and associated heat is added to the soil at night. When hourly measured data areused, rainfall is added at the measurement time. Depending on air temperature,
diﬀerent amounts of heat could be added to the soil. A better estimate heat energy
added to the soil with rain for rain during daylight periods may be possible if thetemperature of the precipitation was assumed to be the same as the wetbulb temperature rather than as air temperature.
This study has concentrated on use of calculated air temperatures and solar
radiation. Hourly wind velocity is also calculated from daily wind runs. High wind
Table 5
Mean diﬀerences and standard deviations for simulated temperatures (%) by soil depth at the Washington site where measured and calculated hourly meteorological data were used (n=1523 per depth)a
Depth (m)Mean diﬀerence ( C)

S.D.

0.00
0.08
0.15
0.25
0.38
0.53
0.69
0.84
1.07
1.37
1.68

1.7
1.8
1.8
1.9
1.9
2.0
2.1
2.1
2.2
2.3
2.5

2.90
2.30
2.09
2.03
2.06
2.15
2.27
2.40
2.59
2.84
3.06

a
The diﬀerence was calculated as soil temperatures simulated using measured hourly soil temperatures
simulated using calculated hourly meteorological data. S.D. is standard deviation of the diﬀerences.

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

237

Fig. 8. Cumulative probability distribution of relative errors in simulated soil temperatures at Mississippi
for ridge positions using measured and calculated hourly weather data. Relative error is calculated as(simulated with measured hourly–simulated with calculated hourly)/simulated with measured hourly.

velocities can increase evapotranspiration and hence aﬀect soil temperature. The
model, however, does account for wind velocity to calculate heat transfer between
the soil and air. Inability to model this eﬀectively can lead to additional errors during
periods of high wind. This should possibly be addressed in the future.

8. Summary and conclusion
The simulation results and comparison with measured soil temperature data using
measured hourly weather data show that 2DSOIL calculates soil temperatures wellon bare soil at both sites with low average bias and error. The comparison of simulation and measurements was carried out for a long time sequence of measured
values (150 days). The errors were largest during and after rainfall periods where
hourly soil temperature values were available. The diﬀerences were also largest at
the soil surface.
Use of calculated hourly weather data can lead to errors in simulation of soil
temperatures. These errors are propagated in the simulation and aﬀect subsurface aswell as surface layers temperature simulations. The largest errors were at the surface
where there is a boundary between the soil and air. The additional error due to theuse of calculated hourly weather data can range from 2 to 3 C This is in addition to238

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

the error of 1–2 C due to the simulation model. Soil temperatures were also underpredicted by 2–3 C when calculated weather data were used. Depending on the
purpose of the simulations, this error may be acceptable. The range of errors due to
use of calculated weather data is not much larger than the error due to the simulation method and so should be acceptable for many purposes. Spatial variability insoil temperature in a ﬁeld, however, may be 3–4 C due to variation in irrigation andwind (Yates et al., 1988). This is a similar range in error. If it is important to predict
noon or afternoon soil or air temperatures for the purpose of estimating plant stressthen measured air temperatures and radiation would be more desirable.

Acknowledgements
The authors gratefully acknowledge the assistance of Dr. Steven Evett and Dr.
Jonathan Ephrath who have reviewed earlier versions of this manuscript.
References
Acock, B.A., Pachepsky, Ya., Timlin, D.J., 2000. Soil–atmosphere Boundary Setting: SetSurf modules.
(Chapter 5). In: Timlin, D.J., Pachepsky, Ya., Van Genuchten, M.Th. (Eds.), 2DSOIL, a Modular
Simulator of Soil and Root Processes. Remote Sensing and Modeling Laboratory Publication No. 2,
USDA Agriculture Research Service, Beltsville, MD.
Acock, B., Trent, A., 1991. The Soybean Simulator, GLYCIM: Documentation for the Modular Version
91. Agric. Exp. Stn., University of Idaho, Moscow, ID.
Bowers, S.A., Hanks, R.J., 1965. Reﬂection of radiant energy from soils. Soil Sci. 100, 130–138.
Budyko, M.I., 1974. Climate and Life. Miller, D.H. (Ed., translated). Academic Press, New York.
Choi, J.S., Fermanian, T.W., Wehner, D.J., Spomer, L.A., 1988. Eﬀect of temperature, moisture, and soil
texture on DCPA Degradation. Agron. J. 80, 108–113.
Driedonks, A.G.M. 1981. Dynamics of the Well Mixed Boundary Layer. Scientiﬁc Report 81–2, Royal
Netherlands Meteorological Institute, De Bilt.
De Vries, D.A., 1963. Thermal properties of soils. In: van Wijk, W.R. (Ed.), Physics of the Plant Environment. North-Holland, Amsterdam, pp. 210–235.
Ephrath, J., Goudriaan, E.J., Marani, A., 1996. Modelling diurnal patterns of air temperature, radiation,wind speed and relative humidity by equations from daily characteristics. Agric. Sys. 51, 377–393.
Flerchinger, G.N., Saxton, K.E., 1989. Simultaneous heat and water model of a freezing snow-residualsoil system II. Field veriﬁcation. Transactions of the ASAE 32, 573–578.
Floyd, R.B., Braddock, R.D., 1984. A simple method for ﬁtting average diurnal temperature curves.
Agric. For. Meteorology 38, 217–229.
Geng, S., Penning de Vries, F.W.T., Supit, I., 1985. Analysis and simulation of weather variables. II.
Temperature and solar radiation. Simulation report CABO-TT No. 5 Wageningen, The Netherlands.
Khorsandi, F.M., Boone, Y.L., Weerakkody, G., Whisler, F.D., 1996. Validation of the soil temperature
subroutine HEAT in the cotton simulation model GOSSYM. Agron J. 89, 415–420.
Koski, A.J., Street, J.R., Danneberger, T.K., 1988. Prediction of Kentucky bluegrass root growth using
degree-day accumulation. Crop Sci. 28, 848–850.
Leenhardt, D., 1995. Errors in the estimation of soil water properties and their propagation through ahydrological model. Soil Use and Managment 11, 15–21.
Linacre, E.T., 1968. Estimating the net radiation ﬂux. Agr. Meterol. 5, 49–63.
Miller, D.H., 1981. Energy at the Surface of the Earth. Academic Press, New York.
Penman, H.L., 1963. Vegetation and Hydrology. C.A.B. Tech. Comm. 53. Commonwealth Agricultural
Bureaux, Harpenden, Hertfordshire, UK.

D.J. Timlin et al. / Agricultural Systems 72 (2002) 215–239

239

Robertson, G.W., Russelo, G.A., 1968. Astrometeorological Estimator. Agric. Meteorol. Tech. Bull. 14.
Plant Res. Inst., C.E.F. Ottowa 3, Ontario.
Šimunek, J., Vogel, T., van Genuchten, M.Th., 1994. The SWMS_2D code for simulating water ﬂow andsolute transport in two-dimensional variably saturated media. Version 1.2. Res. Rep. 132. United States
Salinity Laboratory, Riverside, CA.
Smithsonian Meteorological Tables. 1963. 6th Edition. Lust, R.J. (Ed.). Smithsonian Institution,
Washington, DC.
Timlin, D.J., 1996. Pachepsky, Ya., Acock, B.A. A design for a modular, generic soil simulator to interface with plant models. Agronomy J. 88, 162–169.
van Genuchten, M.Th, 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
Wösten, J.H.M., van Genuchten, M.Th, 1988. Using texture and other soil properties to predict theunsaturated soil hydraulic functions. Soil Sci. Soc. Am. J. 52, 1762–1770.
Weiss, A., 1977. Algorithms for the calculation of moist air properties on a hand calculator. Trans.
A.S.A.E 20, 1133–1136.
Yates, S.R., Warrick, A.W., Matthias, A.D., Musil, S., 1988. Spatial variability of remotely sensed surface
temperatures at ﬁeld scale. Soil Soc. Am. J. 52, 40–45.