Special Section: Patterns

Chris an Steenpass*
Jan Vanderborght
Michael Herbst
Jirka Šimůnek
Harry Vereecken

Es ma ng Soil Hydraulic
Proper es from Infrared
Measurements of Soil Surface
Temperatures and TDR Data
The spa otemporal development of soil surface temperatures (SST) depends on water
availability in the near-surface soil layer. Because the soil loses latent heat during evaporaon and water available for evapora on depends on soil hydraulic proper es (SHP), thetemporal variability of SST should contain informa on about the near-surface SHP. The
objec ve of this study was to inves gate the uncertain es of SHP derived from SST. The
HYDRUS-1D code coupled with a global op mizer (DREAM) was used to inversely es mate van Genuchten–Mualem parameters from infrared-measured SST and me domainreflectometry (TDR)-measured water contents. This approach was tested using synthe cand real data, collected during September 2008 from a harrowed silty loam field plot in
Selhausen, Germany. The synthe c data illustrated that SHP can be derived from SST andthat addi onal soil water content measurements reduce the uncertainty of the es mated
SHP. Unlike for the synthe c experiment with a ver cally homogeneous soil profile, alayered soil profile had to be assumed to derive SHP from the real data. Therefore, theuncertainty of SHP derived from real data was considerably larger. Water reten on curvesof undisturbed soil cores were similar to those es mated from SST and TDR data for the
deeper undisturbed soil. The reten on curves derived from SST and TDR data for the harrowed topsoil layer were typical for a coarse-textured soil and deviated considerably fromthe reten on curves of soil cores, which were typical for a fine-textured soil and similar to
those from the subsoil.
Abbrevia ons: MO, mul objec ve; SO, single objec ve; SST, soil surface temperatures; TDR, me
domain reflectometry.

An 11- d field experiment was
per formed to estimate van
Genuchten–Mualem hydraulic parameters from infrared (IR)-measured
soil surface temperatures and watercontents measured by me domainreflectometry (TDR). Either IR dataor IR and TDR data were used for
optimization. Uncertainties in theparameters estimated from both
approaches were characterized witha global op mizer.

C. Steenpass, J. Vanderborght, M. Herbst, and H.
Vereecken, Ins tute of Chemistry and Dynamics
of the Geosphere (ICG)–Agrosphere, ResearchCenter Jülich, Jülich, Germany; and J. Šimůnek,Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92507. *Corresponding
author (c.steenpass@fz-juelich.de).
Vadose Zone J. 9:910–924
doi:10.2136/vzj2009.0176
Received 30 Nov. 2009.
Published online 20 Oct. 2010.
© Soil Science Society of America
5585 Guilford Rd. Madison, WI 53711 USA.

All rights reserved. No part of this periodical may
be reproduced or transmi ed in any form or by anymeans, electronic or mechanical, including photocopying, recording, or any informa on storage andretrieval system, without permission in wri ng fromthe publisher.

2010, Vol. 9

The dynamics of soil water content at the soil surface provides important
insights into physical processes governing the soils’ interactions with the atmosphere. These
interactions are important in a variety of models that predict soil water and energy balances and plant growth (Vereecken et al., 2008). Soil hydraulic properties are important
controlling factors for water and energy exchange between the soil and the atmosphereand are required to properly simulate these processes.
Classical methods usually identify soil hydraulic properties representative for small soilvolumes (Dane and Topp, 2002). Such measurements are costly and time consuming andneed to be repeated for a large number of soil cores and locations, however, because soil
hydraulic properties vary considerably in space. This variation is caused not only by spatial
heterogeneities within the soil structure but also by soil compaction or soil tillage (Mapa etal., 1986; Murphy et al., 1993; Osunbitan et al., 2005). On the other hand, remote sensing
techniques provide an opportunity to observe the state of the soil surface for relatively large
areas and may thus be used to derive effective soil hydraulic properties at a much larger scale(Burke et al., 1997; Chanzy et al., 1995; Demarty et al., 2005; Das and Mohanty, 2006).
Camillo et al. (1986) and Burke et al. (1998) used ground-based passive microwave data to
inversely estimate soil hydraulic properties. They concluded that hydraulic properties may
be inferred from ground-based remotely sensed data if the data set covers a broad range ofwater contents from dry to wet. Additionally, Camillo et al. (1986) and Burke et al. (1997,1998) concluded that a short period of intensive measurements (Camillo et al. [1986] used
3 d, Burke et al. [1998] used 15 d) was sufficient to derive hydraulic properties.
Several methods allow the inversion of field-measured or remotely sensed data.
Optimization algorithms, such as the Shuffled Complex Evolution algorithm developed
at the University of Arizona (SCE-UA; Duan et al., 1992), are designed to reliably findwww.VadoseZoneJournal.org | 910

the global minimum of the objective function in a space of optimized parameters. Because measurements are generally error-prone,
however, the “best” parameter set corresponding to the global minimum is uncertain, which may lead to a considerable uncertainty
in the model outputs (Vrugt et al., 2003).
In previous studies, either observed soil surface temperatures
were used to infer evapotranspiration rates (Ben-Asher et al.,
1983; Olioso et al., 1996, 1999; Mauser and Schädlich, 1998) or
measured evaporation fluxes were used to estimate soil hydraulicproperties (Jhorar et al., 2002). Additionally, soil surface temperatures strongly depend on the water content in the top few
centimeters of the soil. This means that soil surface temperatures
may be used to directly infer soil hydraulic properties (Demarty
et al., 2005).
To estimate soil hydraulic properties from measured surface temperatures, it is important that water and energy fluxes, as well as
their mutual interactions, are adequately modeled. Philip and deVries (1957) pointed out that the movement of water vapor in
the soil depends on both thermal and hydraulic gradients. The
influence of water vapor in the soil on soil water dynamics has
been studied in numerous works (e.g., Rose, 1968a,b; Jackson etal., 1973; Cahill and Parlange, 1998; Scanlon et al., 2003; Gossand Madlinger, 2007; Shokri et al., 2008). The general agreementin these studies is that water vapor fluxes become more important
for the total soil water dynamics as the soil dries out. Milly (1982),
Kondo et al. (1990), Scanlon and Milly (1994), Saito et al. (2006),
and Bitelli et al. (2008), among others, therefore included water
vapor flow in their models to obtain improved predictions of temperatures, water dynamics, and evaporation rates.
We developed a modeling approach for the inverse estimation of
the van Genuchten–Mualem parameters of soil hydraulic properties (van Genuchten, 1980) from relatively short periods (following
Camillo et al., 1986; Burke et al., 1998) of ground-based, remotelysensed, infrared-based measurements of soil surface temperatures. For that purpose, the optimization algorithm was linked
to a hydrologic model that considers coupled movement of waterand energy (Saito et al., 2006). The use of a recently developed
optimization algorithm (Vrugt et al., 2008a) allowed determination of the information content of soil surface temperatures for
estimation of soil hydraulic parameters. In addition, we tested
whether combining local water content measurements with soil
surface temperatures would improve our estimates of soil hydraulicproperties. Because subsurface hydraulic properties, and the corresponding soil water dynamics, are important for water availability
and heat transport close to the soil surface, time domain reflectometry (TDR) measurements in the 7- and 15-cm depths were
made. The general applicability of the approach was first tested
using numerically generated data. The same approach was then also
applied to a real data set. This study was conducted (i) to analyze
the information content of soil surface temperatures for estimationof soil hydraulic properties with short measurement periods, (ii)to characterize uncertainties in soil hydraulic properties estimatedfrom SST and soil water content data, and (iii) to determine theeffect of soil tillage on near-surface soil hydraulic properties.

6 Materials and Methods
For forward modeling, Version 3.0 of the well-tested hydrologic
model HYDRUS-1D (Šimůnek et al., 2005, 2008) was used. The
HYDRUS-1D code numerically solves the Richards equation forvariably saturated water flow and convection–dispersion-type
equations for heat transport. The governing equations for waterflow and heat transport are solved numerically using Galerkin-typelinear finite element schemes.
For inverse estimation of soil hydraulic properties, soil surface
temperatures have to be modeled as accurately as possible. Various
studies (e.g., de Vries, 1975; Milly, 1984; Cahill and Parlange, 1998;
Saito et al., 2006) have pointed out that soil temperature calculations can be erroneous if thermal and isothermal water vapor fluxes
are neglected. Especially in the upper part of the soil profile, whichdries out very fast, water vapor plays a major role in the total water
movement (Cahill and Parlange, 1998). Following the approach of
Saito et al. (2006), we used a modified HYDRUS-1D Version 3.0
that accounts for water vapor flow and the effects of temperature
gradients on water flow and uses an energy balance equation as theupper boundary condition for the coupled heat and water flow equations. Only basic equations are given below. For a comprehensive
description, see Saito et al. (2006). Similar approaches were developed by Camillo et al. (1983), Camillo and Gurney (1986), Nassar
and Horton (1989, 1992), and Milly (1982, 1984).

Model Descrip on
Liquid Water and Water Vapor Movement
Following Saito et al. (2006), the liquid water, ql (m s−1), and water
vapor, qv (m s−1), fluxes (expressed in terms of equivalent volumesof liquid water) in soils can be described as⎛ ∂h ⎞∂Tq l = q l h + q l T = −K l h ⎜⎜⎜ + 1⎟⎟⎟ − K l T⎝ ∂z ⎠
∂z

[1]

and

q v = q v h + q v T = −K v h∂h∂T− K vT∂z∂z[2]where Klh and Kvh are the isothermal hydraulic conductivities forliquid water and water vapor fluxes (m s−1), respectively, KlT and
KvT are the thermal hydraulic conductivities for liquid water andwater vapor fluxes (m2 s−1 K−1), respectively, qlh and qlT are theisothermal and thermal liquid water fluxes (m s−1), respectively,
qvh and qvT are the isothermal and thermal water vapor fluxes (ms−1), respectively, h is the pressure head (m), z (m) is the verticalcoordinate, and T (K) is the temperature.

www.VadoseZoneJournal.org | 911

Inserting the liquid water and water vapor fluxes into a massconservation equation yields the governing equation for onedimensional water flow in variably saturated porous media, which
describes the change in the total water content with time (Saito
et al., 2006):
∂θ
∂ ⎡
∂h
∂T
=
⎢ K lh+ K lh ( h ) + K lT
∂ t ∂ z ⎣⎢
∂z
∂z
[3]
∂h
∂T ⎤
+K v h+ K vT⎥∂z∂ z ⎥⎦
where θ is the total volumetric water content, consisting of thesum of the volumetric liquid water content, θ l, and the water vapor
content, θ v, (θ = θ l + θ v, both expressed as a volume of liquid water
per volume of bulk soil [m3 m−3]), and t (s) is time.

type of soil and accounts for temperature-induced changes in thesoil water retention.
The isothermal hydraulic conductivity (Klh [m s−1]) itself is derived
using the pore-size distribution model of Mualem (1976) from theanalytical retention function of van Genuchten (1980):
m ⎤2⎡
K l h ( h ) = K s Sel ⎢1 − (1 − Se1/m ) ⎥
⎢⎣
⎥⎦

[7]

where Ks is the saturated hydraulic conductivity (m s−1), Se is theeffective saturation (dimensionless), and l and m are empiricalparameters.

Heat Transport
Thermal and Isothermal Water Vapor Conduc vi es
The thermal (KvT) water vapor hydraulic conductivity is definedas (Saito et al., 2006)

q h = −λ ( θ l )dρ
DK vT =
η H r sv
dT
ρw

[4]

where D is the vapor diff usivity in the soil (m2 s−1), ρ

w is the den−3
sity of water (kg m ), η is the enhancement factor (dimensionless)(Cass et al., 1984) accounting for increased thermal vapor fluxes
due to increased temperature gradients in the air phase, Hr is therelative humidity (dimensionless), and ρ sv is the saturated vapor
density (kg m−3).
The isothermal vapor hydraulic conductivity (Kvh) (m s−1) isdefined as

K vh =MgD
ρ sv H r
ρw
RT

[5]

tational acceleration (m s−2), and R is the universal gas constant(J mol−1 K−1).

Thermal and Isothermal Liquid Water Conduc vi es
Following Noborio et al. (1996b), the thermal hydraulic conductivity, KlT (m2 K−1 s−1), for liquid water flux is[6]

where Klh is the isothermal liquid hydraulic conductivity (m s−1), γ
is the surface tension (kg s−2), γ 0 is the surface tension of soil waterat 25°C (kg s−2), and G wT (dimensionless) is a gain factor (Nimmo
and Miller, 1986) that depends on the soil water content and the∂T+ C wTq l + C vTq v + L0 q v∂z[8]where qh (J m−2 s−1) is the total heat flux, λ(θ l) (J m−1 s−1 K−1)
is the apparent thermal conductivity, which is a function of thevolumetric water content θ l, T is the temperature (K), Cw and Cv(J m−3 K−1) are the volumetric heat capacities for water and vapor,respectively, ql (m s−1) is the liquid water flux, qv (m s−1) is thewater vapor flux, which is expressed as an equivalent volume fluxof liquid water (m s−1), and L0 (J m−3) is the volumetric latent heatof vaporization of liquid water (Saito et al., 2006).
Combining the heat flux and the energy conservation equations
yields the governing equation for heat flow, describing the change
in heat storage with time:∂ C pT

where M is the molecular weight of water (kg mol−1), g is the gravi-⎛
1 d γ ⎞⎟⎟⎟
K l T = K l h ⎜⎜⎜ hG w T⎜⎝
γ 0 d T ⎠⎟Based on the mathematical model of Philip and de Vries (1957)
and Philip (1957), the soil heat flux is described using

∂t

+ L0

∂ q lT∂θ v ∂ ⎡∂T ⎤=⎢λ ( θ l )⎥ −C w⎢⎥∂t ∂z ⎣∂z ⎦∂z∂qv
∂ q vT
− L0
−C v
∂z
∂z

[9]

where Cp (J m−3 K−1) is the heat capacity of the soil, consisting of
the volumetric heat capacities of the liquid water, water vapor, andthe solid phases (Saito et al., 2006).

Apparent Thermal Conduc vity
According to Šimůnek et al. (2005) and Hopmans et al. (2002),
the apparent thermal conductivity, λ(θ l), can be described as a
combination of the thermal conductivity of the porous medium
and the macrodispersion, defined as a linear function of the waterflow velocity:λ ( θ l ) = λ 0 ( θ l )+ β C w q l

www.VadoseZoneJournal.org | 912[10]

where β is the thermal dispersivity (m). Chung and Horton (1987)
described the thermal conductivity, λ0(θ1), according to
λ 0 ( θ l ) = b1 +b2 θ l +b3 θ 0.5
l[11]

where b1, b2 , and b3 are empirical parameters (W m−1 K−1).

volumetric heat capacity of the air (J m−3 K−1), rH (s m−1) is the
aerodynamic resistance to heat flow, and Ta (K) is the air temperature. Because vapor flow in the soil is modeled as a functionof vapor density gradients and vapor diff usivity in the soil, thereduction of the vapor flux due to resistance to vapor flow within
the soil is explicitly accounted for and not parameterized as a soil
resistance in a model that predicts vapor transfer from the soil
surface to the atmosphere (van de Griend and Owe, 1994).

Energy Balance Boundary Condi on
The solution of the partial differential equations describing waterand heat flux requires knowledge of the initial and boundary conditions. For the boundary condition at the soil–air interface, thesurface heat and water fluxes are computed using a surface energybalance from measured meteorologic data (i.e., rainfall, air temperature, wind speed, air humidity, and radiation) and the statevariables T, h, and θ v at the soil surface. Because the calculated
mass and energy fluxes, which are used as boundary conditions
for the solution of the water and heat flow equations, depend on
the soil surface state variables and thus the solution of the waterand heat flow equations themselves, the surface fluxes and statevariables had to be derived iteratively. Two iteration loops were
found to be a good tradeoff between accuracy and computational
effort (Saito et al., 2006).

The energy balance equation was used to provide surface fluxes forheat as well as for water (Camillo et al., 1983). Solving the energybalance for each time step during the numerical simulation yields
the heat flux density G and the evaporation rate E, which are usedas heat and water flux boundary conditions for the top of the soil
profi le. In the case of precipitation, the sum of precipitation andevaporation rates is used as a flux boundary condition for water flow.
The calculated evaporation flux, E, defines the water and vapor flux
at the soil surface:E = q +q
v
l
ρw
z =0

[15]

and the calculated heat flux, G, defines the heat flux into the soil:
The surface energy balance is given as

G = Rnet − H − Lw E

[12]where G (W m−2) is the heat flux into the soil, Rnet is the net radiation (W m−2), H is the sensible heat flux (W m−2), Lw is the latent
heat of vaporization (J kg−1), and E is the evaporation rate (kg m−2s−1). The variables G and Rnet are defined to be positive downward(into the soil profi le) and LwE and H are defi ned to be positive
upward (out of the soil profi le).
The net radiation, Rnet, is the sum of the net shortwave and net
longwave radiation (sum of incoming and outgoing shortwave andlongwave radiation) and was measured half-hourly for the durationof the field experiment. Following van Bavel and Hillel (1976), thelatent heat flux, LwE, and the sensible heat flux, H, from a bare soil
can be defined as⎛ ρ − ρ a ⎞⎟⎟Lw E = −Lw ⎜⎜⎜ s⎜⎝ ra ⎠⎟⎟

[13]

Ts −TarH[14]

H = Cawhere ρ s (kg m−3) and Ts (K) are the water vapor density and temperature at the soil surface, ρ a (kg m−3) is the water vapor density
measured in the air at a certain height above the soil surface, ra (sm−1) is the atmospheric resistance to water vapor flow, Ca is theG = − q h z =0

[16]

The negative sign follows from definitions of G and qh; fluxes are
positive upward and negative downward.
The evaporation rate, E, includes both the liquid water flux toward
the soil surface, which is evaporated at the soil surface, and thewater vapor flux at the soil surface. The total heat flux into thesoil, G, comprises the conductive and convective heat fluxes andthe latent heat flux due to vapor flux, L0 qv, in the soil profi le; Gis written as
G = −L0 q v − qh ,nonlatent[17]

where qh,nonlatent is the heat flux in the soil that is not caused bylatent heat flow. Thus, qh,nonlatent increases (in absolute terms)
with q v to provide energy for evaporation within the soil profile.
In a dry soil where most of the water is evaporated within thesoil profile and the evaporation rate approximates E ? ρ wq v| z
= 0 , it follows from the surface energy balance that −qh,nonlatent= R net − H.
For a wet soil profi le, qv|z = 0 ? 0, E ? ρ wql|z = 0 and −qh,nonlatent= Rnet − H − L0 ql|z = 0. In a wet soil, the latent heat is consumed
at the soil surface. Thus, including vapor flow in a soil water flow
model leads to the prediction of higher soil surface temperatures
when vapor flow in the soil contributes significantly to the evaporation rate.

www.VadoseZoneJournal.org | 913

Field Measurements
Field data were collected from 15 to 26 Sept. 2008 on a bare soilat the Selhausen field test site, near the Research Center Jülich, inJülich, Germany. The soil at the experimental test site had been
kept bare since September 2006 and was regularly treated withherbicides to prevent the accumulation of weeds. According to the
USDA, the soil was classified as a silt loam (13% sand, 70% silt,
and 17% clay; Herbst et al., 2009). Bauer et al. (2008) measured anorganic matter content of 2.2 to 2.5% (v/v) for the upper 20 cm of
the soil. From the site, small undisturbed soil cores (100 cm3) were
sampled. Water retention curves of the samples were determinedin the laboratory using sand suction tables and pressure cells. The
saturated hydraulic conductivity was derived using a constant-head
permeameter. The soil hydraulic parameters, i.e., the Mualem–van
Genuchten parameters and the saturated conductivity that were
derived from the laboratory measurements, are given in Table 1.
We installed 72 TDR probes (two-rod TDR probes, 25-cm rod
length, 2.3-cm spacing between the rods, 11.5-m coaxial cable
length) at 36 locations in a field plot of 36 m2 (6 by 6 m), which
was harrowed before the installation of the TDR probes. The bulk
density of the harrowed upper part of the soil (0–10 cm) was 1.35
g cm−3 and 1.69 g cm−3 for the lower (10–90 cm), untilled partof the soil. Note that the bulk density of the harrowed top layer
was considerably smaller than the bulk density of the soil samples
that were taken from the topsoil layer of an undisturbed field plot
and from which the soil hydraulic parameters were derived in thelaboratory (Table 1). To install the TDR probes, we excavatedthree trenches with six small pits at each side of each trench. One
week before the start of the experiment, we horizontally installed
two TDR probes at depths of 7 and 15 cm in each of the 36 pits.
Because the coaxial TDR cables might disturb the infrared surface
temperature measurements, we used the trenches to bury the cables(at 10-cm depths) to the TDR instruments outside the test plot.
Before installation, all TDR probes were calibrated with water–
air measurements. The installed TDR probes were connected toa TDR 100 cable tester and the recorded waveforms were stored
on a CR1000 data logger (both Campbell Scientific, Logan, UT).
We used the PCTDR soft ware (Campbell Scientific) to analyze
the waveforms for the dielectrical permittivity and the empirical
relationship of Topp et al. (1980) to convert the dielectrical permittivity to water content. The soil water content was monitored
on a half-hourly basis. Additionally, at each location, soil temperatures at the 3- and 6-cm depths were measured manually every 6 hwith a TypeE thermocouple (Li 8100–201, LI-COR Biosciences,
Lincoln, NE).
Additionally, an infrared (IR) camera (VarioCAM, Infra Tec
GmbH, Dresden, Germany) was installed to measure the soil
surface temperature every 5 min. The IR camera measures thebrightness temperature in the spectral range from 8 to 14 μm. It
has a resolution of 320 by 240 pixels and an absolute measurement accuracy of ±1.5 K in the temperature range from −10 toTable 1. Soil hydraulic parameters† for the upper (0–10 cm) and lower(10–90 cm) parts of the soil profi le as determined in the laboratory
using the sand suction table and the pressure cell method; 100-cm3 soil
samples were taken before the harrowing of the plot in an undisturbed
part of the soil profi le. The saturated hydraulic conductivity was determined using the constant-head method.

cm

θs
3
—— cm cm−3 ——

cm−1

0–10

0.00000.35300.0042

10–90

0.00000.37190.0690

Depth

θr

α

n

Ks

BD

cm d−1g cm−31.324291.6800

1.61401.252550.8800

1.6920† θ r and θ s, residual and saturated volumetric water contents, respectively; α and n,
shape parameters; Ks, saturated hydraulic conductivity; BD, bulk density.

50°C. The IR camera was installed on an auto hoist 11 m above
the ground to cover the 6- by 6-m measurement plot.
During the field experiment, the sky was partially covered withclouds, which can be seen in the net radiation data in Fig. 1. The
air temperature, relative humidity (combined temperature–air
humidity sensor CS215, Campbell Scientific), wind speed (CSAT3,
Campbell Scientific), and net radiation (1 x SP Lite, Kipp & Zonen,
Delft, the Netherlands) were monitored by an automatic weather
station installed close to the test site. All meteorologic data were
measured 2 m above the ground. Figure 1 shows the meteorologic

Fig. 1. Meteorological data measured at the Selhausen test site during
the field experiment from 15 to 26 Sept. 2008 (Day of Year [DOY]259–270). The upper plot shows the irrigation events, the middle plot
the wind speed (solid black line) and the relative humidity (dotted redline), and the lower plot the air temperature (dotted red line) and the
net radiation (solid black line).

www.VadoseZoneJournal.org | 914

standard data for the duration of the experiment. Half-hourly valuesof the meteorologic data were used as input data for the model.

The a priori distribution was assumed to be a uniform distribution
with ranges of the hydraulic soil parameters that cover most soil
types (see Table 2). The DREAM algorithm runs multiple Markov
chains in parallel for different starting points and every chain (aparent) generates new candidate points (proposals) for every single
model evaluation. Whether or not the new candidate parameter
sets are accepted is based on the log-likelihood of the parameterset given the measured data. By updating the parameter sets inthe chains based on their log-likelihoods, the distribution of theparameter sets converges to the posterior distribution given themeasured data set.

Parameteriza on and Inverse Op miza on
A 90-cm soil profi le was discretized for numerical simulations
into 100 finite elements and 101 nodes. Zero pressure head andtemperature gradients (i.e., free drainage) were used as bottom
boundary conditions for water flow and heat transport, respectively. The energy balance and the soil surface boundary condition
were described above. Because pressure head and temperature
gradients increase toward the top of the soil profi le, the spatial
discretization was gradually refined toward the soil surface, where
the grid size was 0.3 cm. We used the Chung and Horton (1987)
parameters of a loamy soil (b1 = 24.3, b2 = 39.3, b3 = 153.4, allW m−1 K−1) implemented in HYDRUS-1D to describe the thermal properties of the soil at the test site. The thermal dispersivity,β, has not been measured often for unsaturated soils because
it only becomes important when the water flux velocity is high(Hopmans et al., 2002). We took β = 5 cm, which is the default
value in HYDRUS-1D. The initial conditions for the numericalsimulations were derived from TDR water content measurements
at 7-, 15-, 45-, 60-, and 90-cm depths and temperature measurements at 3-, 6-, 15-, 45-, 60-, and 90-cm depths at the beginningof the experiment. Water contents at 45-, 60-, and 90-cm depths
and temperatures at 15-, 4-5, 60-, and 90-cm depths were measured
a few meters away from the experimental plot with horizontally
installed TDR and temperature probes. These measurements wereused only to defi ne the initial conditions as well as possible butwere not further used in the optimization process. Water content
measurements at the 7- and 15-cm depths were used for inversion.

The posterior distributions of soil hydraulic parameters were
derived for two different data sets. The first data set contained only
SST measurements, which were spatially averaged across the fieldplot. For this single-objective (SO) approach, the log-likelihood
function lSST was defined as()

l SST Θ Yˆ SST , ζ, φ =N ⎡ln ( 2 π ) + ln ( σ )⎦⎤
2 ⎣
⎤2
N ⎡y
⎢⎣ i,SST (Θ ζ, φ ) − yˆi,SST ⎥⎦− 0.5∑
2
σSST
i=1
−

[18]

where Θ is a vector of model parameters (van Genuchten–Mualem
parameters, in this case), ζ and φ are the measured initial andboundary conditions, respectively, ŶSST is a vector of the observed
system behavior (in this case, the SST measurements), N denotes
the number of SST measurements, yi,SST and ŷi,SST are the ith
SST measurement value and the SST model prediction, respectively, and σSST is the standard deviation of residuals between themeasured and simulated soil surface temperatures for the optimal
parameter set, i.e., for the parameter set with the greatest likelihood. The optimal parameter set or the parameter set with thegreatest log-likelihood does not depend on σSST; however, σSST
determines the spreading of the posterior parameter distribution.
Therefore, a two-step approach was used to determine the posterior
parameter distribution. In the first step, the optimal parameter
set was determined using a dummy value of σSST. For the optimal
parameter set, σSST was calculated from the residuals between the

The Differential Evolution Adaptive Metropolis (DREAM)
algorithm (Vrugt et al., 2008a,b) was used to estimate the van
Genuchten–Mualem soil hydraulic parameters from measured
soil surface temperatures and water contents and to determine
their posterior probability distribution given the observed data set.
The algorithm starts with an initial population of points within
the possible parameter space (the a priori parameter distribution),
which is sampled by the Latin hypercube method to effectively represent the parameter space with a low number of sampling points.

Table 2. The reference set of the van Genuchten–Mualem soil hydraulic parameters† used to generate measurements for the synthetic experiment andthe optimized values obtained using single-objective (SO) and multiobjective (MO) approaches. Numbers in parentheses are the lower and upper bandsof the 95% confidence interval. The parameter range gives the upper and lower limits allowed for each parameter during the optimization.
θs3——————————— cm cm−3 ——————————

cm−1

Reference parameter

0.09000.40000.09001.900010.0000

Soil surface
temperature (SO)0.0910 (0.0371–0.1051)

0.3838 (0.2905–0.4438)

0.1015 (0.0500–0.1908)

1.7069 (1.3283–1.9771)

8.688 (1.3489–42.4521)

Soil surface temperature+ water content (MO)0.0860 (0.0582–0.0976)

0.4014 (0.3258–0.4466)

0.1063 (0.0689–0.1764)

1.8375 (1.5345–1.9882)

11.6054 (2.4394–28.7872)

Parameter range

0.000–0.150

0.250–0.450

0.003–0.200

1.100–2.000

1.00–500.00

Measured data

θr

α

n

Ks
cm d−1† θ r and θ s, residual and saturated volumetric water contents, respectively; α and n, shape parameters; Ks, saturated hydraulic conductivity.

www.VadoseZoneJournal.org | 915

modeled and measured surface temperatures (σSST = 1.95°C). Thisvalue was subsequently used to derive the posterior parameter distribution. It must be noted that σSST comprises both measurement
and model errors.
The second data set contained soil surface temperature and watercontent measurements at 7- and 15-cm depths, expressed as themean of the 7- and 15-cm-depth TDR measurements. A multiobjective optimization (MO) was used in this case, with one objectivefunction defi ned using surface temperatures and the other two
using water contents at different depths. The log-likelihood was
calculated for water contents measured at different depths j = 7,
15 cm asNl WC j Θ Yˆ WC j , ζ, φ = − ⎡⎣ ln (2 π) + ln (σ)⎤⎦
2⎡⎤2N ⎢ yi,WC (Θ ζ, φ) − yˆi,WC ⎥jj⎦− 0.5∑ ⎣
2
σ
i=1
WC j

(

)

[19]

where σWCj are the standard deviations of the residuals between
the measured and simulated water contents at a depth of j = 7
and j = 15 cm, which were calculated with the same approach as
described for σSST. Again, the values of σWCj (0.005 cm3 cm−3 ata depth of 7 cm and 0.006 cm3 cm−3 at a depth of 15 cm) compriseboth measurement and model errors.
The total log likelihood was obtained asl total = l SST + l WC7 + l WC15

[20]

Based on the calculated likelihood of every parameter set, theMarkov chain either moves to a new location when the likelihood
is higher or remains at the old position when the likelihood is lower.
Additionally, a random change in the chain’s position can occur:

α = min ⎡⎣ exp ( l new − l old ) ,1⎤⎦[21]

where α is the Metropolis acceptance probability, and lnew andlold represent the likelihoods of the proposal and the parent position of the Markov chain, respectively. Equation [21] shows that ahigher likelihood in the proposal position (lnew) than the parent
position (lold) results in α = 1, forcing the Markov chain to moveto the proposal position. For α < 1, a random number r is drawn
from a uniform distribution between 0 and 1 and the chain moves
to the proposal position if r < α. This Random Walk Metropolis
method makes the DREAM optimization algorithm an efficientand robust sampler of the parameter space. For inverse modeling,a total of 50,000 model evaluations distributed on five central processing units was used. This setup resulted in approximately 2 d ofcomputation time for an optimization run.
To test the ability of the DREAM algorithm to infer the van
Genuchten–Mualem soil hydraulic parameters from soil surface

temperature data, a synthetic experiment was performed first. A
HYDRUS-1D forward simulation with known van Genuchten–
Mualem parameters (Table 2), subsequently referred to as thereference parameters, was first performed to generate synthetic soil
surface temperatures and water contents at two different depths.
The synthetic data set was generated using the same initial andboundary conditions, the same meteorologic conditions, and thesame flow domain as for the real Selhausen data set. A random
error with known standard deviation (σSST = 0.75°C, σWC7 =0.005 cm3 cm−3, σWC15 = 0.005 cm3 cm−3) was added to the simulated data, which were then used as a synthetic measurement data
set for the inverse optimization of the van Genuchten–Mualem
parameters of soil hydraulic properties.

6 Results and Discussion
Synthe c Data Experiment

In Fig. 2, posterior probability density functions (pdfs) of theparameters obtained using the SO and MO approaches areshown. The horizontal axes of the plots correspond with theinitial parameter distribution. The width of the posterior pdf
compared with the width of the a priori parameter distribution
is a measure for how strongly the posterior parameter distribution is constrained by the measurements. The true and optimized
parameter sets are given in Table 2 together with the 95% confidence intervals (CI95). The width of the posterior distributionsor the CI95 compared with the width of the prior distributions(Fig. 2; Table 2) show a relatively good ability to identify theparameters θ r and Ks. In relation to the initial a priori parameter
ranges, confidence intervals for parameters n, θ s , and α are larger.
Table 2 shows that both optimizations based on measurementsof only soil surface temperatures (SO) or on measurements ofsoil surface temperatures and subsurface water contents (MO)produced parameter sets close to the reference parameters. The
univariate posterior distributions of the optimized parameters
do not contain information about correlations between parameters or clustering of parameter combinations in the parameter
space. Instead of investigating the parameter uncertainty, theuncertainty intervals can also be plotted directly around thehydraulic functions. In Fig. 3, the optimized water retention
curves (obtained for the optimized parameters) are shown
together with the 95% confidence range of water contents forthe given pressure heads. A visual comparison of optimized soil
water retention curves obtained with the SO and MO approaches(Fig. 3) illustrates that the use of additional water content measurements yielded slightly better results. The main difference
between the two approaches is the uncertainty in the optimized
retention curves, which is clearly larger for the SO than for theMO approach, especially for lower pressure heads, e.g., at h =−1000 cm, the 95% uncertainty interval of the water content, θ,ranges from 0.063 to 0.131 cm3 cm−3 for the SO approach, andis somewhat narrower, i.e., from 0.072 to 0.106 cm 3 cm−3 forthe MO approach.

www.VadoseZoneJournal.org | 916

Fig. 2. Probability density functions (pdfs) of the optimized van Genuchten–Mualem parameters for thesynthetic experiment: (A) shape parameter n, (B) saturated hydraulic conductivity (Ks), (C) shape parameterα, (D) residual volumetric water contents (θr), and (E)saturated volumetric water contents (θs) for the synthetic experiment. Red (wide bars) histograms show the
pdfs resulting from single-objective (SO) optimization(soil surface temperature data only), blue (narrow bars)
histograms show the pdfs for a multiobjective (MO)optimization (soil surface temperature and water content
data). In the upper part of every plot, the 95% confidence
intervals are given for the SO (upper line) and MO (lowerline) approaches. Vertical black lines indicate the true
parameter values.

The optimization results of the van Genuchten–Mualem parameters were then used in the HYDRUS-1D forward simulation tocompare the modeled soil water contents and soil surface temperatures with those generated using the reference parameter set.
Figure 4 shows a comparison between the generated (using the reference parameter set) and modeled (using optimized parameters)
soil surface temperatures and water contents. Both MO and SO
approaches provided a similar fit of soil surface temperatures with
only marginal differences. The obtained root mean squared errors(RMSEs), 0.7425°C for SO and 0.7427°C for MO, are close to
the random error that was added to the forward simulated temperatures. Differences in the soil water dynamics are much moresignificant. Figure 4 shows that the MO approach using measurements of both soil surface temperatures and water contents yieldeda better fit to measured values at both depths (RMSE for 7 cm =0.0051 cm3 cm−3, RMSE for 15 cm = 0.0046 cm3 cm−3, which is
close to the random error that was added to the forward simulations). Parameters estimated from only soil surface temperature
data resulted in an overestimation of soil water contents at both
depths. This bias led to a larger RMSE (RMSE ? 0.01 cm3 cm−3)
for the SO approach (Fig. 5) than the random error (0.005 cm3
cm−3); however, the prediction of the soil water dynamics with the
optimized parameters from the soil surface temperature measurements was still reasonable (Fig. 4).

www.VadoseZoneJournal.org | 917

Fig. 3. Soil water retention curves for the synthetic experiment: (A)the red line is the optimized water retention curve based on only soil
surface temperature measurements (SST) and (B) the blue line is the
optimized water retention curve based on soil surface temperature
and water content measurements (SST + WC). The dashed black line
represents the reference soil water retention curve. The shaded areas
represent the 95% confidence intervals of the water contents.

The synthetic experiment was designed to illustrate the usefulness of soil surface temperatures to derive the parametersof soil hydraulic properties. The results show that soil surface
temperatures contain information that can be used to derive
the parameters of soil hydraulic properties. Even though uncertainties in the optimized parameters become smaller when watercontent data are used in addition to soil surface temperatures, thereference parameter set used to generate the synthetic measurement data was reproduced reasonably well using both the SO andthe MO approaches.

Field Data Experiment
The measured field data were used to parameterize a modelaccounting for a soil profi le with vertically uniform hydraulicproperties. Th is assumption implies that an effective parameterset can be found representing a vertically heterogeneous soil profi le and adequately predicting soil surface temperatures and soil

Fig. 4. (A) Soil water contents (WC) and (B) soil surface temperatures (SST) calculated for Days of the Year (DOY) 259 to 270 using
the van Genuchten–Mualem parameters estimated using the singleobjective (SO) approach (SST data, red lines) and the multiobjective(MO) approach (SST + WC data, blue lines). Solid and open circles
show synthetic water contents at 7- and 15-cm depths, respectively,
and squares denote synthetic soil surface temperatures.

water contents at different depths. The soil hydraulic parameters
were optimized using the DREAM algorithm and, analogous
to the synthetic experiment, both SO and MO approaches. The
final optimized values are shown in Table 3 and a priori parameter
ranges during the optimization process are given in Table 4.
Figure 6 reveals that the two estimated parameter sets (using
either the SO or MO approaches) cannot describe soil water
Fig. 5. Modeled temperatures and water contents(y axis) plotted against thesynthetic measurements (xaxis). The results for parameters estimated inversely
using the single-objective(SO) approach (from soil
surface temperatures, SST)
are shown in the top row.
The results for parameters
estimated inversely using
the multiobjective (MO)approach (from soil surface
temperatures combined with
subsurface soil water contents, SST + WC) are shownin the bottom row.

www.VadoseZoneJournal.org | 918

Table 3. Van Genuchten–Mualem parameter sets† of an assumed
homogeneous soil profi le that were estimated from measured datausing single-objective (SO) and multiobjective (MO) approaches.
Measured data

θr

θs

α

n

Ks

—— cm3 cm−3 ——cm−10.09040.30590.13101.306154.9872

Soil surface
0.1381temperature + watercontent (MO)0.33130.1965

1.19622.8787

Soil surface
temperature (SO)cm d−1† θ r and θ s, residual and saturated volumetric water contents, respectively; αand n, shape parameters; Ks, saturated hydraulic conductivity.

content dynamics properly. This is similar to the findings of Yanget al. (2005), who showed that effective soil hydraulic properties
cannot represent layered soils well. They pointed out that soil
hydraulic parameters cannot be adjusted in such a way that a layered, heterogeneous soil profi le can be approximated assuming ahomogeneous soil profi le. Because the estimation of effective van
Genuchten–Mualem hydraulic properties for a homogeneous soil
profi le did not yield sensible results, parameter uncertainties were
not analyzed.
Before the installation of the TDR probes, the plot was harrowed. As
a consequence, the bulk density of the topsoil layer (1.35 g cm−3) was
considerably lower than that of the subsoil (1.69 g cm−3), indicatingthat different soil hydraulic properties would have to be assigned tothe two layers. We therefore subdivided the soil profile into a layer
representing the harrowed upper part of the soil (0–10 cm) and alayer representing the non-harrowed subsoil (10–90-cm depth).
A sequential iterative approach was used to determine the parameters of the two soil layers. In the first step, a homogeneous soil
profi le was assumed. Its hydraulic properties were estimated using
soil surface temperatures and water contents at a depth of 7 cm.
During the second step, the soil profi le was divided into two horizons, with the parameters estimated in the first step assigned toFig. 6. Soil water content measurements at the 15- (circles) and7-cm depths (filled dots) and predictions (lines) for Days of the Year
(DOY) 259 to 270 assuming a homogeneous soil profile that wasparameterized using only soil surface temperature measurements(SST) and using a combination of soil surface temperatures and soil
water content measurements (SST + WC).

the surface soil layer (0–10 cm) and the parameters of the subsurface layer (10–90 cm) derived from water content measurements
at the 15-cm depth. Because soil surface temperatures and watercontents are influenced by the hydraulic properties of both layers,
the hydraulic properties of the individual soil layers were successively updated in subsequent steps, while keeping the parametersof the other layer fi xed, until the estimated hydraulic propertiesof both layers stopped changing between subsequent steps. The
iterative process converged after three iterations. Th is approach
represents an MO inversion because both water contents and soil
surface temperatures were used during this parameter optimization process. For each step, the prior distribution of the parametersof the soil layer of which the parameters were estimated was set tothe initial prior distributions. For the last step, the posterior distributions of the van Genuchten–Mualem soil hydraulic parameters
for the upper part of the soil were estimated using only soil surface
temperatures (the SO approach) and soil surface temperature andwater content measurements at the 7-cm depth (the MO approach),while keeping the soil hydraulic parameters for the deeper part ofthe soil profi le constant. Except for the last step of the iteration
process, we used the DREAM algorithm only until convergence,

Table 4. Van Genuchten–Mualem parameter sets† and their 95% confidence intervals of a layered soil profi le (0–10 and 10–90 cm) that were estimatedfrom measured data using single-objective (SO) and multiobjective (MO) approaches. The soil hydraulic properties of the deeper soil layer (10–90 cm)
were estimated iteratively using time domain reflectometry measurements at the 15-cm depth. Numbers in parentheses are the lower and upper bandsof the 95% confidence interval. The parameter range gives the upper and lower limits allowed for each parameter during the optimization.
Depth

Measured data

cm
0–10

soil surface
temperature (SO)

θr

θs
3
—————————— cm cm−3 —————————

α

n

cm−1

Ks
cm d−10.1083 (0.0322–0.1500) 0.2811 (0.2102–0.4427) 0.1232 (0.0128–0.1916) 1.4669 (1.1971–1.9671) 16.2378 (1.0454–108.8930)

soil surface temperature 0.1132 (0.0732–0.1181) 0.3097 (0.2720–0.3413) 0.1862 (0.1413–0.1998) 1.8664 (1.6146–2.0000) 5.4761 (4.3361–35.7108)+ water content(7 cm) (MO)10–90

water content (15 cm) 0.0538

0.39020.03391.570123.6872

0 –90

parameter range

0.250–0.450

0.003–0.200

1.100–2.000

1.00–500.00

0.000–0.150† θ r and θ s, residual and saturated volumetric water contents, respectively; α and n, shape parameters; Ks, saturated hydraulic conductivity.

www.VadoseZoneJournal.org | 919

which greatly sped up computational time. Therefore, no confidence intervals for the best estimate of the hydraulic properties
for the lower soil layer (0–90 cm) were obtained (Table 4). The
focus was mainly on the upper part of the soil profi le because ouradditional modeling showed that soil surface temperatures did notcontain sufficient information to infer soil water dynamics of thesecond, deeper layer of the soil profi le.
Soil hydraulic parameters estimated from field-measured surface
temperatures and from both temperatures and water contents(Table 4) differ more than those estimated from the synthetic
data. The best estimate of the van Genuchten–Mualem parametern shows the largest difference between the two approaches (SO:n = 1.42; MO: n = 1.89). The posterior parameter distribution of n remained nearly uniform, however, indicating that ncannot be identified with the SO approach. The uncertainty of
all optimized parameters was considerably larger when only surface temperatures were used, compared with the combined useof surface temperatures and soil water content measurements(e.g., CI95 for Ks with SO: 1.05–108.89 cm d−1; CI95 for Ks
with MO: 4.34–35.71 cm d−1; Fig. 7). For the SO approach, theparameter uncertainty was also considerably larger than that in
the synthetic experiment. The CI95 in the synthetic experimentcalculated from the SO approach for θ r ranged from 0.037 toFig. 7. Probability density functions (pdfs) of the optimized van Genuchten–Mualem parameters for thesynthetic experiment: (A) shape parameter n, (B) saturated hydraulic conductivity (Ks), (C) shape parameterα, (D) residual volumetric water contents (θr), and (E)saturated volumetric water contents (θs) for the 0- to
10-cm layer for the field experiment. Red (wide bars)
histograms show the pdfs resulting from single-objective(SO) optimization (soil surface temperature data only),
blue (narrow bars) histograms show the pdfs for a multiobjective (MO) optimization (soil surface temperature
and water content data). In the upper part of every plot,
the 95% confidence intervals are given for the SO (upper
red line) and MO (lower blue line) approaches.

www.VadoseZoneJournal.org | 920

0.105 cm3 cm−3, whereas it ranged from 0.03
to 0.15 cm3 cm−3 in the field experiment. For
the α and n parameters, the additional use
of water content data led to a concentrationof the probability mass at the upper boundary of the parameter ranges. The value ranges
for these parameters could not be extended
due to numerical instabilities and unrealistic
behavior in the model.
Fig. 8. Soil water retention curves for the optimized van Genuchten–Mualem parameter sets: (A)Water retention curves resulting from theoptimized for the soil layer from 10 to 10 cm depth using only soil surface temperature (SST)
best van Genuchten–Mualem parameter setsmeasurements; (B) optimized for the soil layer from 0 to 10 cm depth using both soil surface temperature and water content (SST + WC) measurements at the 7-cm depth; and (C) determined
and their 95% confidence intervals are shownfor the soil layer from 10 to 90 cm depth using the sand bed and the pressure plate method (dashedin Fig. 8 for the SO and MO approaches. Asline) and optimized using time domain reflectometry measurements (solid line) at the 15-cm
expected, examining the posterior parameter
depth. Shaded areas represent the 95% confidence intervals in terms of water contents.
distributions (Fig. 7), the uncertainty in the
optimized retention curve was considerably
reduced when soil water content data were considered together
Due to the lower bulk density of the harrowed upper soil layer
with soil surface temperature measurements in the optimization
than the untilled deeper layer, the saturated hydraulic conprocess. Figure 8 also compares the estimated water retention
ductivity and the saturated water content were expected to be
curves against those determined in the laboratory (Table 1).
larger for the upper soil layer; however, fitted Ks and θ s in theBoth the SO- and MO-derived water retention curves for the
lower soil layer were higher than in the upper soil layer (Table
surface soil layer deviated considerably from the retention curve4). Th is may be due to a lack of information in the wet part ofdetermined in the laboratory. Both estimated retention curves
the water retention curve, especially for the upper layer of theare characteristic of a soil with a coarser texture (i.e., larger αsoil (measurements at 7 cm) where measured water contents did
parameter). The optimized retention curve for the topsoil layer
not exceed values around 0.24 cm3 cm−3. As a consequence, thealso deviates considerably from the retention curve derived for
estimated saturated hydraulic conductivity and water content
the deeper soil layer, which shows more resemblance to thatare an extrapolation of hydraulic properties that were derived
determined in the laboratory. Th is indicates larger pore sizes infor drier soil conditions. For aggregated or tilled soils with a conthe topsoil layer than what would normally be observed for a soil
siderable interaggregate porosity, this extrapolation may lead towith this texture. The bulk density of the topsoil was decreased
a considerably underestimation of the real saturated hydraulicby harrowing at the beginning of the experiment, which may
conductivity and water content. A measurement period coverexplain the effect on α . Despite the larger uncertainty of theing a wider range of dry and wet soil conditions may give better
SO-estimated hydraulic parameters, it is still possible to inferresults concerning the saturated hydraulic conductivity and thethat the hydraulic properties of the disturbed topsoil layer are
saturated water content.
considerably different from the hydraulic properties determinedon undisturbed soil samples.
Figure 9 shows the water contents and surface temperatures
simulated using parameters determined by both SO and MO

Fig. 9. (A) Soil water contents (WC) and (B) soil surface temperatures (SST) calculated using the van Genuchten–Mualem parameters estimated from
SST (red lines) and SST + WC (blue lines) data. (B) also shows soil surface temperatures (dotted line) calculated using the soil hydraulic properties
determined in the laboratory on undisturbed soil samples. Solid and open circles represent measured water contents at depths of 7 and 15 cm, respectively, and squares denote measured soil temperatures. Measured surface temperatures and water contents are expressed as an average of all plots.

www.VadoseZoneJournal.org | 921

approaches. Both approaches
show a good agreement between
measured and modeled watercontents and soil surface temperatures. The RMSE for thewater content at the 7-cm depth
is 0.012 cm 3 cm−3 for the SO
approach and 0.011 cm 3 cm−3
for the MO approach. At the
15-cm depth, the RMSE calculated from the SO approach is
0.0046 cm 3 cm−3 and 0.0037
cm3 cm−3 for the MO approach.
The soil water dynamics atthe 7-cm depth following therain event around Day of Year
(DOY) 267 of the field experiment were less well described,
which probably indicates inadequacies of the unimodal van
Genuchten–Mualem model.

Fig. 10. Modeled temperatures and water contents (y axis) plotted against the field measurements (x axis). The
results for parameters estimated inversely from soil surface temperatures (SST) are shown in the top row. The
results for parameters estimated using both soil surface temperatures and soil water contents (SST + WC) are
shown in the bottom row.

Figure 9 also shows surface
temperatures simulated using
soi l hydrau l ic para meters
determined in the laboratory. Simulated surface temperatures
significantly underestimated peak measured temperatures, whichshows that not only soil water dynamics but also soil surface
temperatures strongly depend on the soil tillage.
As already noted above for the synthetic experiment, including
soil water content measurements in the parameter estimation
improved the fit of the soil water dynamics (see Fig. 10); however,
the improvement for the real data set was not as pronounced as it
was for the synthetic experiment. This can be seen in the RMSE
for the water content at, e.g., the 7-cm depth, which is reduced
from 0.012 cm3 cm−3 for the SO approach to 0.011 cm3 cm−3
for the MO approach. Because the optimal parameter sets thatwere obtained from the SO and MO approaches were similar, thereduction in the RMSE of the soil water content fit by the MO
approach is small (Fig. 10).
As a validation of our simulations and estimated soil hydraulicparameters, we determined whether soil temperatures measured
at depths of 3 and 6 cm could be simulated properly with theparameterized model. Because both best parameter sets (fromthe SO and MO approaches) predicted nearly identical temperatures at both depths, only soil temperatures simulatedusing the MO parameters are plotted in Fig. 11. It shows thatmodeled and measured soil temperatures are in good agreement.
Dampening of the amplitude of the diurnal temperature variation from the soil surface to depths of 3 and 6 cm was predicted
well by the model.

6 Summary and Conclusions
Estimating soil hydraulic properties using inverse modeling of theinformational content of the soil surface temperature was analyzed
in this study. To quantify the constraints that soil surface temperature has on the estimation of soil hydraulic parameters, we estimatedsoil hydraulic parameters using a SO approach, which considered
only soil surface temperatures, and a MO approach, which considered both soil surface temperatures and TDR-measured soil water

Fig. 11. Measured and simulated soil temperatures at depths of 3
and 6 cm for Days of Year (DOY) 259 to 270. Simulated values wereobtained with the van Genuchten–Mualem parameter set determinedusing the multiobjective approach.

www.VadoseZoneJournal.org | 922

contents. A synthetic experiment showed that, although the soil
surface temperature contains enough information to estimate soil
hydraulic parameters, their uncertainties can be reduced by additionally considering also water content measurements.
A field experiment was then performed to test whether hydraulic properties can be estimated using this approach under field
conditions. An attempt to find effective soil hydraulic properties
characterizing the entire soil profi le that could be used to describe
field-measured soil temperatures and water contents at different
depths was not successful. Th is shows that models that try todescribe water dynamics in layered soils using a single effective
hydraulic parameter set (as most land-surface models do) may
lead to wrong predictions of soil surface states (water contents
and temperatures) and fluxes across the soil–atmosphere boundary,which are important for characterizing soil–atmosphere interactions. Considering a two-layered soil profi le, it was found that thehydraulic properties of the soil surface layer were considerably different (e.g., α = 0.187 cm−1 for the MO approach) from those of
the deeper soil (α = 0.0339 cm−1) and from the hydraulic properties derived in the laboratory on undisturbed soil columns (α =0.069 cm−1). The hydraulic properties that were optimized for the
top layer of the silty loam soil were more typical of a coarse sandy
soil and characterized by a high value of the α parameter of the van
Genuchten–Mualem function. This indicates the important effectof soil tillage on the hydraulic properties of the soil surface layer
and its consequences for soil–atmosphere interactions.
Because the soil surface temperature measurements did not containenough information to infer the hydraulic properties of the deeper
part of the soil profi le, we parameterized the deeper part of the soil
using soil water contents measured using TDR probes installed at a
depth of 15 cm. Compared with the shallow soil layer, the obtainedwater retention curves for the deeper soil layer corresponded betterwith those derived on laboratory columns. This indicates that soil
surface temperature measurements should be combined withinformation about subsurface soil water contents to derive the soil
hydraulic properties of the topsoil layer. Soil surface temperatures
can be observed at a relatively large scale and could be used to infer
information about the effects not only of spatial, but also temporal,
variations in surface soil properties, such as those due to soil tillage, on soil hydraulic properties and soil–atmosphere interactions.
This information, however, should be combined with spatial information about the subsurface soil water contents, which could beobtained using networks of soil water content sensors (Rosenbaum
et al., 2010). How to combine these two sources of information
at larger scales and how to address the problems involved in thevariability of parameters at various scales requires further research.
The soil hydraulic parameters of the upper soil layer, inversely
estimated only from soil surface temperature measurements,
were associated with considerable uncertainty (e.g., CI95 for n= 1.33–1.98). Nevertheless, the soil hydraulic parameters andcorresponding water retention curves that were obtained using
only soil surface temperatures or a combination of temperatureand soil water content measurements were similar and could
both be distinguished from the soil hydraulic properties of the
deeper undisturbed soil. Comparing the results obtained using
the synthetic and real field experiments indicates that synthetic
experiments are important tools to outline the general applicabilityof a particular method but can only provide a first approximation
about the actual performance of the method under field conditions.
Considering a layered soil profi le, both parameter sets estimatedfrom field-measured data (using MO and SO approaches) led tosimilar predictions of near-surface soil water dynamics (RMSE forθ at 7 cm = 0012 cm3 cm−3 [SO approach] and 0.011 cm3 cm−3
[MO approach]), indicating that the proposed method can beapplied under field conditions. Drying–wetting cycles are required,however, to derive soil hydraulic properties from soil surface temperature measurements because it is not possible to relate the soil
surface temperature to the soil hydraulic properties in dry soils
where heat transfer is controlled by heat conduction.
Finally, the range of soil hydrologic conditions, i.e., the range ofsoil water contents and the range of boundary conditions (radiation, evaporation rates, and temperatures), during our experiment
was limited. Observations of soil surface temperatures in different
plots with different soil surface structure and for a longer time
period, containing more than one drying–wetting cycle, may considerably reduce the parameter uncertainty and will be the subjectof further research.

Acknowledgments

This work was financed by the Transregio Collaborative Research Center 32, Patterns in Soil–Vegetation–Atmosphere Systems: Monitoring, Modelling and Data
Assimilation. We want to thank Jasper Vrugt, University of California, Irvine, forproviding the source code of the global optimizer DREAM.

References
Bauer, J., M. Herbst, J.A. Huisman, L. Weihermüller, and H. Vereecken. 2008.
Sensi vity of simulated soil heterotrophic respira on to temperatureand moisture reduc on func ons. Geoderma 145:17–27.
Ben-Asher, J., A.D. Ma hias, and A.W. Warrick. 1983. Assessment of evaporaon from bare soil by infrared thermometry. Soil Sci. Soc. Am. J. 47:185–191.
Bitelli, M., F. Ventura, G.S. Campbell, L.R. Snyder, F. Gallega , and P.R. Pisa.
2008. Coupling of heat, water vapor and liquid water fluxes to computeevapora on in bare soils. J. Hydrol. 362:191–205.
Burke, E.J., R.J. Gurney, L.P. Simmonds, and T.J. Jackson. 1997. Calibra ng asoil water and energy budget model with remotely sensed data to obtainquan ta ve informa on about the soil. Water Resour. Res. 33:1689–1697.
Burke, E.J., R.J. Gurney, L.P. Simmonds, and P.E. O’Neill. 1998. Using a modelling approach to predict soil hydraulic proper es from passive microwave measurements. IEEE Trans. Geosci. Remote Sens. 36:454–462.
Cahill, A.T., and M.B. Parlange. 1998. On water vapor transport in field soils.
Water Resour. Res. 34:731–739.
Camillo, P.J., and R.J. Gurney. 1986. A resistance parameter for bare soil
evapora on models. Soil Sci. 141:95–105.
Camillo, P.J., R.J. Gurney, and T.J. Schmugge. 1983. A soil and atmosphericboundary layer model for evapora on and soil moisture studies. Water
Resour. Res. 19:371–380.
Camillo, P.J., P.E. O’Neill, and R.J. Gurney. 1986. Es ma ng soil hydraulic parameters using passive microwave data. IEEE Trans. Geosci. Remote Sens.
GE-24:930–936.

www.VadoseZoneJournal.org | 923

Cass, A., G.S. Campbell, and T.L. Jones. 1984. Enhancement of thermal water
vapor diﬀusion in soils. Soil Sci. Soc. Am. J. 48:25–32.
Chanzy, A., L. Bruckler, and A. Perrier. 1995. Soil evapora on monitoring: A
possible synergism of microwave and infrared remote sensing. J. Hydrol.
165:235–259.
Chung, S.-O., and R. Horton. 1987. Soil heat and water flow with a par alsurface mulch. Water Resour. Res. 23:2175–2186.
Dane, J.H., and G.C. Topp (ed.). 2002. Methods of soil analysis. Part 4. Physical methods. SSSA Book Ser. 5. SSSA, Madison, WI.
Das, N.N., and B.P. Mohanty. 2006. Root zone soil moisture assessment using
remote sensing and vadose zone modeling. Vadose Zone J. 5:296–307.
Demarty, J., C. O lé, J. Braud, A. Olioso, J.P. Frangi, H.V. Gupta, and L.A. Bas das. 2005. Constraining a physically based Soil–Vegeta on–Atmosphere
Transfer model with surface water content and thermal infrared brightness temperature measurements using a mul objec ve approach. Water Resour. Res. 41:W01011, doi:10.1029/2004WR003695.
de Vries, D.A. 1975. Heat transfer in soils. p. 5–28. In D.A. de Vries and N.H.
Afgan (ed.) Heat and mass transfer in the biosphere: I. Transfer processes
in the plant environment. Scripta Book Co., Washington, DC.
Duan, Q., V.K. Gupta, and S. Sorooshian. 1992. A shuﬄed complex evolu on
approach for eﬀec ve and eﬃcient global minimiza on. J. Op m. Theory
Appl. 76:501–521.
Goss, K.-U., and M. Madlinger. 2007. Es ma on of water transport based onin situ measurements of rela ve humidity and temperature in a dry Tanzanian soil. Water Resour. Res. 43:W05433, doi:10.1029/2006WR005197.
Herbst, M., N. Prolingheuer, A. Graf, J.A. Huisman, L. Weihermüller, and J.
Vanderborght. 2009. Characteriza on and understanding of bare soil respira on spa al variability at plot scale. Vadose Zone J. 8:762–771.
Hopmans, J.W., J. Šimůnek, and K.L. Bristow. 2002. Indirect es ma on ofsoil thermal proper es and water flux using heat pulse probe measurement: Geometry and dispersion eﬀects. Water Resour. Res. 38(1):1006,
doi:10.1029/2000WR000071.
Jackson, R.D., B.A. Kimball, R.J. Reginato, and F.S. Nakayama. 1973. Diurnalsoil-water evapora on: Time–depth–flux pa erns. Soil Sci. Soc. Am. Proc.
37:505–509.
Jhorar, R.K., W.G.M. Bas aanssen, R.A. Feddes, and J.C. Van Dam. 2002. Inversely es ma ng soil hydraulic func ons using evapotranspira on fluxes. J. Hydrol. 258:198–213.
Kondo, J., N. Saigusa, and T. Sato. 1990. A parameteriza on of evapora on
from bare soil surfaces. J. Appl. Meteorol. 29:385–389.
Mapa, R.B., R.E. Green, and L. Santo. 1986. Temporal variability of soil hydraulic proper es with we ng and drying subsequent to llage. Soil Sci.
Soc. Am. J. 50:1133–1138.
Mauser, W., and S. Schädlich. 1998. Modelling the spa al distribu on of
evapotranspira on on diﬀerent scales using remote sensing data. J. Hydrol. 212/213:250–267.
Milly, P.C.D. 1982. Moisture and heat transport in hystere c, inhomogeneous, porous media: A matric head-based formula on and a numericalmodel. Water Resour. Res. 18:489–498.
Milly, P.C.D. 1984. A simula on analysis of thermal eﬀects on evapora on
from soil. Water Resour. Res. 20:1087–1098.
Mualem, Y. 1976. A new model for predic ng the hydraulic conduc vity ofunsaturated porous media. Water Resour. Res. 12:513–521.
Murphy, B.W., T.B. Koen, B.A. Jones, and L.M. Huxedurp. 1993. Temporal varia on of hydraulic proper es for some soils with fragile structure. Aust. J.
Soil Res. 31:179–197.
Nassar, I.N., and R. Horton. 1989. Water transport in unsaturated nonisothermal salty soils: II. Theore cal development. Soil Sci. Soc. Am. J.
53:1330–1337.
Nassar, I.N., and R. Horton. 1992. Simultaneous transfer of heat, water, andsolute in porous media: I. Theore cal development. Soil Sci. Soc. Am. J.
56:1350–1356.
Nimmo, J.R., and E.E. Miller. 1986. The temperature dependence of isothermal moisture vs. poten al characteris cs in soils. Soil Sci. Soc. Am.
J. 50:1105–1113.
Noborio, K., K.J. McInnes, and J.L. Heilmann. 1996b. Two-dimensional model
for water, heat, and solute transport in furrow-irrigated soil: II. Field evalua on. Soil Sci. Soc. Am. J. 60:1010–1021.
Olioso, A., H. Chauki, D. Courault, and J.-P. Wigneron. 1999. Es ma on of
evapotranspira on and photosynthesis by assimila on of remote sensing data into SVAT models. Remote Sens. Environ. 68:341–356.
Olioso, A., O. Taconet, and M.B. Mehrez. 1996. Es ma on of heat and mass
fluxes from IR brightness temperature. IEEE Trans. Geosci. Remote Sens.
34:1184–1190.
Osunbitan, J.A., D.J. Oyedele, and K.O. Adekalu. 2005. Tillage eﬀects on bulk
density, hydraulic conduc vity and strength of a loamy sand soil in southwestern Nigeria. Soil Tillage Res. 82:57–64.
Philip, J.R. 1957. Evapora on, and moisture and heat fields in the soil. J. Meteorol. 14:354–366.

Philip, J.R., and D.A. de Vries. 1957. Moisture movement in porous materials under temperature gradient. Trans. Am. Geophys. Union 38:222–232.
Rose, C.W. 1968a. Water transport in soil with a daily temperature wave: I.
Theory and experiment. Aust. J. Soil Res. 6:32–44.
Rose, C.W. 1968b. Water transport in soil with a daily temperature wave: II.
Analysis. Aust. J. Soil Res. 6:45–57.
Rosenbaum, U., J.A. Huisman, A. Weuthen, H. Vereecken, and H. Bogena.
2010. Quan fica on of sensor-to-sensor variability of the ECH2O EC-5,
TE, and 5TE sensors in dielectric liquids. Vadose Zone J. 9:181–186.
Saito, H., J. Šimůnek, and B.P. Mohanty. 2006. Numerical analysis of coupled water, vapor, and heat transport in the vadose zone. Vadose Zone
J. 5:784–800.
Scanlon, B., K. Keese, R.C. Reedy, J. Šimůnek, and B. Andraski. 2003. Variaons in flow and transport in thick desert vadose zones in response topaleoclima c forcing (0–90 kyr): Monitoring, modeling, and uncertaines. Water Resour. Res. 39(7):1179, doi:10.1029/2002WR001604.
Scanlon, B.R., and P.C.D. Milly. 1994. Water and heat fluxes in desert soils: 2.
Numerical simula ons. Water Resour. Res. 30:721–733.
Shokri, N., P. Lehmann, P. Vontobel, and D. Or. 2008. Drying front and watercontent dynamics during evapora on from sand delineated by neutron radiography. Water Resour. Res. 44:W06418, doi:10.1029/2007WR006385.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2005. The HYDRUS-1D so ware package for simula ng the one-dimensional movement of water,
heat, and mul ple solutes in variably-saturated media. Version 3.0. HYDRUS So w. Ser. 1. Dep. of Environ. Sci., Univ. of Calif., Riverside.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2008. Development and applica ons of the HYDRUS and STANMOD so ware packages and related
codes. Vadose Zone J. 7:587–600.
Topp, G.C., J.L. Davis, and A.P. Annan. 1980. Electromagne c determina onof soil water content: Measurement in coaxial transmission lines. Water
Resour. Res. 16:574–582.
van Bavel, C.H.M., and D.I. Hillel. 1976. Calcula ng poten al and actual evapora on from a bare soil surface by simula on of concurrent flow of waterand heat. Agric. For. Meteorol. 17:453–476.
van de Griend, A., and M. Owe. 1994. Bare soil surface resistance to evapora on by vapor diﬀusion under semiarid condi ons. Water Resour. Res.
30:181–188.
van Genuchten, M.Th. 1980. A closed-form equa on for predic ng the hydraulic conduc vity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
Vereecken, H., J.A. Huisman, H. Bogena, J. Vanderborght, J.A. Vrugt, andJ.W. Hopmans. 2008. On the value of soil moisture measurements
in vadose zone hydrology: A review. Water Resour. Res. 44:W00D06,
doi:10.1029/2008WR006829.
Vrugt, J., H.V. Gupta, W. Bouten, and S. Sorooshian. 2003. A Shuﬄed Complex
Evolu on Metropolis algorithm for op miza on and uncertainty assessment of hydrologic model parameters. Water Resour. Res. 39(8):1201,
doi:10.1029/2002WR001642.
Vrugt, J., C.J.F. ter Braak, M.P. Clark, J.M. Hyman, and B.A. Robinson. 2008a.
Treatment of input uncertainty in hydrologic modeling: Doing hydrology
backward with Markov chain Monte Carlo simula on. Water Resour. Res.
39(8):1201, doi:10.1029/2002WR001642.
Vrugt, J., C.J.F. ter Braak, H.V. Gupta, and B.A. Robinson. 2008b. Equifinality informal (DREAM) and informal (GLUE) Bayesian approaches in hydrologic
modeling? Stochas c Environ. Res. Risk Assess. 23:1011–1026.
Yang, K., T. Koike, B. Ye, and L. Bas das. 2005. Inverse analysis of the roleof soil ver cal heterogeneity in controlling surface soil state and energy
par on. J. Geophys. Res. 110:D08101, doi:10.1029/2004JD005500.

www.VadoseZoneJournal.org | 924