Freely available online through the author-supported open access option.

Abstract

The spatiotemporal development of soil surface temperatures (SST) depends on water availability in the near-surface soil layer. Because the soil loses latent heat during evaporation and water available for evaporation depends on soil hydraulic properties (SHP), the temporal variability of SST should contain information about the near-surface SHP. The objective of this study was to investigate the uncertainties of SHP derived from SST. The HYDRUS-1D code coupled with a global optimizer (DREAM) was used to inversely estimate van Genuchten–Mualem parameters from infrared-measured SST and time domain reflectometry (TDR)-measured water contents. This approach was tested using synthetic and real data, collected during September 2008 from a harrowed silty loam field plot in Selhausen, Germany. The synthetic data illustrated that SHP can be derived from SST and that additional soil water content measurements reduce the uncertainty of the estimated SHP. Unlike for the synthetic experiment with a vertically homogeneous soil profile, a layered soil profile had to be assumed to derive SHP from the real data. Therefore, the uncertainty of SHP derived from real data was considerably larger. Water retention curves of undisturbed soil cores were similar to those estimated from SST and TDR data for the deeper undisturbed soil. The retention curves derived from SST and TDR data for the harrowed topsoil layer were typical for a coarse-textured soil and deviated considerably from the retention curves of soil cores, which were typical for a fine-textured soil and similar to those from the subsoil.

An 11-d field experiment was performed to estimate van Genuchten–Mualem hydraulic parameters from infrared (IR)-measured soil surface temperatures and water contents measured by time domain reflectometry (TDR). Either IR data or IR and TDR data were used for optimization. Uncertainties in the parameters estimated from both approaches were characterized with a global optimizer.

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 atmosphere and are required to properly simulate these processes.

Classical methods usually identify soil hydraulic properties representative for small soil volumes (Dane and Topp, 2002). Such measurements are costly and time consuming and need 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 et al., 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 of water 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 find 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 hydraulic properties (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 de Vries (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 et al., 1973; Cahill and Parlange, 1998; Scanlon et al., 2003; Goss and Madlinger, 2007; Shokri et al., 2008). The general agreement in 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, remotely sensed, infrared-based measurements of soil surface temperatures. For that purpose, the optimization algorithm was linked to a hydrologic model that considers coupled movement of water and 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 hydraulic properties. 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 estimation of soil hydraulic properties with short measurement periods, (ii) to characterize uncertainties in soil hydraulic properties estimated from SST and soil water content data, and (iii) to determine the effect of soil tillage on near-surface soil hydraulic properties.

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 for variably saturated water flow and convection–dispersion-type equations for heat transport. The governing equations for water flow and heat transport are solved numerically using Galerkin-type linear 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, which dries 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 the upper 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 Description

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 volumes of liquid water) in soils can be described as 
\[q_{1}{=}q_{1h}{+}q_{1T}{=}{-}K_{1h}\left(\frac{{\partial}h}{{\partial}z}{+}1\right){-}K_{1T}\frac{{\partial}T}{{\partial}z}\]
[1]
and 
\[q_{\mathrm{v}}{=}q_{\mathrm{v}h}{+}q_{\mathrm{v}T}{=}{-}K_{\mathrm{v}h}\frac{{\partial}h}{{\partial}z}{-}K_{\mathrm{v}T}\frac{{\partial}T}{{\partial}z}\]
[2]
where Klh and Kvh are the isothermal hydraulic conductivities for liquid water and water vapor fluxes (m s−1), respectively, KlT and KvT are the thermal hydraulic conductivities for liquid water and water vapor fluxes (m2 s−1 K−1), respectively, qlh and qlT are the isothermal and thermal liquid water fluxes (m s−1), respectively, qvh and qvT are the isothermal and thermal water vapor fluxes (m s−1), respectively, h is the pressure head (m), z (m) is the vertical coordinate, and T (K) is the temperature.
Inserting the liquid water and water vapor fluxes into a mass conservation equation yields the governing equation for one-dimensional water flow in variably saturated porous media, which describes the change in the total water content with time (Saito et al., 2006): 
\begin{eqnarray*}&&\frac{{\partial}\mathrm{{\theta}}}{{\partial}t}{=}\frac{{\partial}}{{\partial}z}\left[K_{1h}\frac{{\partial}h}{{\partial}z}{+}K_{1h}(h){+}K_{1T}\frac{{\partial}T}{{\partial}z}\right.\ \\&&\left.\ {+}K_{\mathrm{v}h}\frac{{\partial}h}{{\partial}z}{+}K_{\mathrm{v}T}\frac{{\partial}T}{{\partial}z}\right]\end{eqnarray*}
[3]
where θ is the total volumetric water content, consisting of the sum 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.

Thermal and Isothermal Water Vapor Conductivities

The thermal (KvT) water vapor hydraulic conductivity is defined as (Saito et al., 2006) 
\[K_{\mathrm{v}T}{=}\frac{D}{\mathrm{{\rho}}_{\mathrm{w}}}\mathrm{{\eta}}H_{\mathrm{r}}\frac{\mathrm{d{\rho}}_{\mathrm{sv}}}{\mathrm{d}T}\]
[4]
where D is the vapor diffusivity in the soil (m2 s−1), ρw is the density of water (kg m−3), η 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 the relative humidity (dimensionless), and ρsv is the saturated vapor density (kg m−3).
The isothermal vapor hydraulic conductivity (Kvh) (m s−1) is defined as 
\[K_{\mathrm{v}T}{=}\frac{D}{\mathrm{{\rho}}_{\mathrm{w}}}\mathrm{{\eta}}H_{\mathrm{r}}\frac{Mg}{RT}\]
[5]
where M is the molecular weight of water (kg mol−1), g is the gravitational acceleration (m s−2), and R is the universal gas constant (J mol−1 K−1).

Thermal and Isothermal Liquid Water Conductivities

Following Noborio et al. (1996b), the thermal hydraulic conductivity, KlT (m2 K−1 s−1), for liquid water flux is 
\[K_{1T}{=}K_{1h}\left(hG_{\mathrm{w}T}\frac{1}{\mathrm{{\gamma}}_{0}}\frac{\mathrm{d{\gamma}}}{\mathrm{d}T}\right)\]
[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 water at 25°C (kg s−2), and GwT (dimensionless) is a gain factor (Nimmo and Miller, 1986) that depends on the soil water content and the type of soil and accounts for temperature-induced changes in the soil water retention.
The isothermal hydraulic conductivity (Klh [m s−1]) itself is derived using the pore-size distribution model of Mualem (1976) from the analytical retention function of van Genuchten (1980): 
\[K_{1h}(h){=}K_{\mathrm{s}}S_{\mathrm{e}}^{l}\left[1{-}(1{-}S_{\mathrm{e}}^{1/m})^{m}\right]^{2}\]
[7]
where Ks is the saturated hydraulic conductivity (m s−1), Se is the effective saturation (dimensionless), and l and m are empirical parameters.

Heat Transport

Based on the mathematical model of Philip and de Vries (1957) and Philip (1957), the soil heat flux is described using 
\[q_{\mathrm{h}}{=}{-}\mathrm{{\lambda}}(\mathrm{{\theta}}_{1})\frac{{\partial}T}{{\partial}z}{+}C_{\mathrm{w}}Tq_{1}{+}C_{\mathrm{v}}Tq_{\mathrm{v}}{+}L_{0}q_{\mathrm{v}}\]
[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 the volumetric 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 the water vapor flux, which is expressed as an equivalent volume flux of liquid water (m s−1), and L0 (J m−3) is the volumetric latent heat of 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: 
\begin{eqnarray*}&&\frac{{\partial}C_{\mathrm{p}}T}{{\partial}t}{+}L_{0}\frac{{\partial}\mathrm{{\theta}}_{\mathrm{v}}}{{\partial}t}{=}\frac{{\partial}}{{\partial}z}\left[\mathrm{{\lambda}}(\mathrm{{\theta}}_{1})\frac{{\partial}T}{{\partial}z}\right]{-}C_{\mathrm{w}}\frac{{\partial}q_{1}T}{{\partial}z}\\&&{-}L_{0}\frac{{\partial}q_{\mathrm{v}}}{{\partial}z}{-}C_{\mathrm{v}}\frac{{\partial}q_{\mathrm{v}}T}{{\partial}z}\end{eqnarray*}
[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, and the solid phases (Saito et al., 2006).

Apparent Thermal Conductivity

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 water flow velocity: 
\[\mathrm{{\lambda}}(\mathrm{{\theta}}_{1}){=}\mathrm{{\lambda}}_{0}(\mathrm{{\theta}}_{1}){+}\mathrm{{\beta}}C_{\mathrm{w}}{\vert}q_{1}{\vert}\]
[10]
where β is the thermal dispersivity (m). Chung and Horton (1987) described the thermal conductivity, λ01), according to 
\[\mathrm{{\lambda}}_{0}(\mathrm{{\theta}}_{1}){=}b_{1}{+}b_{2}\mathrm{{\theta}}_{1}{+}b_{3}\mathrm{{\theta}}_{1}^{0.5}\]
[11]
where b1, b2, and b3 are empirical parameters (W m−1 K−1).

Energy Balance Boundary Condition

The solution of the partial differential equations describing water and heat flux requires knowledge of the initial and boundary conditions. For the boundary condition at the soil–air interface, the surface heat and water fluxes are computed using a surface energy balance from measured meteorologic data (i.e., rainfall, air temperature, wind speed, air humidity, and radiation) and the state variables 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 water and heat flow equations themselves, the surface fluxes and state variables 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 surface energy balance is given as 
\[G{=}R_{\mathrm{net}}{-}H{-}L_{\mathrm{w}}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−2 s−1). The variables G and Rnet are defined to be positive downward (into the soil profile) and LwE and H are defined to be positive upward (out of the soil profile).
The net radiation, Rnet, is the sum of the net shortwave and net longwave radiation (sum of incoming and outgoing shortwave and longwave radiation) and was measured half-hourly for the duration of the field experiment. Following van Bavel and Hillel (1976), the latent heat flux, LwE, and the sensible heat flux, H, from a bare soil can be defined as 
\[L_{\mathrm{w}}E{=}{-}L_{\mathrm{w}}\left(\frac{\mathrm{{\rho}}_{\mathrm{s}}{-}\mathrm{{\rho}}_{\mathrm{a}}}{r_{\mathrm{a}}}\right)\]
[13]
 
\[H{=}C_{\mathrm{a}}\frac{T_{\mathrm{s}}{-}T_{\mathrm{a}}}{r_{H}}\]
[14]
where ρ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 (s m−1) is the atmospheric resistance to water vapor flow, Ca is the 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 function of vapor density gradients and vapor diffusivity in the soil, the reduction 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).

The energy balance equation was used to provide surface fluxes for heat as well as for water (Camillo et al., 1983). Solving the energy balance for each time step during the numerical simulation yields the heat flux density G and the evaporation rate E, which are used as heat and water flux boundary conditions for the top of the soil profile. In the case of precipitation, the sum of precipitation and evaporation 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: 
\[\frac{E}{\mathrm{{\rho}}_{\mathrm{w}}}{=}q_{\mathrm{v}}{+}q_{1}{\vert}_{z{=}0}\]
[15]
and the calculated heat flux, G, defines the heat flux into the soil: 
\[G{=}{-}q_{\mathrm{h}}{\vert}_{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 the water vapor flux at the soil surface. The total heat flux into the soil, G, comprises the conductive and convective heat fluxes and the latent heat flux due to vapor flux, L0qv, in the soil profile; G is written as 
\[G{=}{-}L_{0}q_{\mathrm{v}}{-}q_{h,\mathrm{nonlatent}}\]
[17]
where qh,nonlatent is the heat flux in the soil that is not caused by latent heat flow. Thus, qh,nonlatent increases (in absolute terms) with qv to provide energy for evaporation within the soil profile. In a dry soil where most of the water is evaporated within the soil profile and the evaporation rate approximates E ∼ ρwqv|z = 0, it follows from the surface energy balance that −qh,nonlatent = Rnet− H.

For a wet soil profile, qv|z = 0 ∼ 0, E ∼ ρwql|z = 0 and −qh,nonlatent = RnetHL0ql|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.

Field Measurements

Field data were collected from 15 to 26 Sept. 2008 on a bare soil at the Selhausen field test site, near the Research Center Jülich, in Jülich, Germany. The soil at the experimental test site had been kept bare since September 2006 and was regularly treated with herbicides 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 an organic 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 determined in 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 part of 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 the laboratory (Table 1). To install the TDR probes, we excavated three 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 to a TDR 100 cable tester and the recorded waveforms were stored on a CR1000 data logger (both Campbell Scientific, Logan, UT). We used the PCTDR software (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 h with 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 the brightness 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 to 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 with clouds, 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 standard data for the duration of the experiment. Half-hourly values of the meteorologic data were used as input data for the model.

Parameterization and Inverse Optimization

A 90-cm soil profile was discretized for numerical simulations into 100 finite elements and 101 nodes. Zero pressure head and temperature 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 profile, 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, all W 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 numerical simulations 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 beginning of 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 were used only to define the initial conditions as well as possible but were not further used in the optimization process. Water content measurements at the 7- and 15-cm depths were used for inversion.

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. 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 (a parent) 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 parameter set given the measured data. By updating the parameter sets in the chains based on their log-likelihoods, the distribution of the parameter sets converges to the posterior distribution given the measured data set.

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 field plot. For this single-objective (SO) approach, the log-likelihood function lSST was defined as 
\begin{eqnarray*}&&l_{\mathrm{SST}}({\Theta}{\vert}\mathbf{\mathrm{{\hat{Y}}}}_{\mathrm{SST}},\mathrm{{\zeta},{\phi}}){=}\\&&{-}\frac{N}{2}[\mathrm{ln}(2\mathrm{{\pi}}){+}\mathrm{ln}(\mathrm{{\sigma}})]\\&&{-}0.5{{\sum}_{i{=}1}^{N}}\frac{[y_{i,\mathrm{SST}}({\Theta}{\vert}\mathrm{{\zeta},{\phi}}){-}{\hat{y}}_{i,\mathrm{SST}}]^{2}}{\mathrm{{\sigma}}_{\mathrm{SST}}^{2}}\end{eqnarray*}
[18]
where Θ is a vector of model parameters (van Genuchten–Mualem parameters, in this case), ζ and ϕ are the measured initial and boundary 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 the measured 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 the greatest 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 modeled and measured surface temperatures (σSST = 1.95°C). This value 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 water content measurements at 7- and 15-cm depths, expressed as the mean of the 7- and 15-cm-depth TDR measurements. A multiobjective optimization (MO) was used in this case, with one objective function defined 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 as 
\begin{eqnarray*}&&l_{\mathrm{WC}_{j}}({\Theta}{\vert}\mathbf{\mathrm{{\hat{Y}}}}_{\mathrm{WC}_{j}},\mathrm{{\zeta},{\phi}}){=}{-}\frac{N}{2}[\mathrm{ln}(2\mathrm{{\pi}}){+}\mathrm{ln}(\mathrm{{\sigma}})]\\&&{-}0.5{{\sum}_{i{=}1}^{N}}\frac{[y_{i,\mathrm{WC}_{j}}({\Theta}{\vert}\mathrm{{\zeta},{\phi}}){-}{\hat{y}}_{i,\mathrm{WC}_{j}}]^{2}}{\mathrm{{\sigma}}_{\mathrm{WC}_{j}}^{2}}\end{eqnarray*}
[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 at a depth of 7 cm and 0.006 cm3 cm−3 at a depth of 15 cm) comprise both measurement and model errors.
The total log likelihood was obtained as 
\[l_{\mathrm{totol}}{=}l_{\mathrm{SST}}{+}l_{\mathrm{WC}7}{+}l_{\mathrm{WC}15}\]
[20]
Based on the calculated likelihood of every parameter set, the Markov 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: 
\[\mathrm{{\alpha}}{=}\mathrm{min}[\mathrm{exp}(l_{\mathrm{new}}{-}l_{\mathrm{old}}),1]\]
[21]
where α is the Metropolis acceptance probability, and lnew and lold represent the likelihoods of the proposal and the parent position of the Markov chain, respectively. Equation [21] shows that a higher likelihood in the proposal position (lnew) than the parent position (lold) results in α = 1, forcing the Markov chain to move to 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 efficient and 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 of computation 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 the reference 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 and boundary conditions, the same meteorologic conditions, and the same 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.

Results and Discussion

Synthetic Data Experiment

In Fig. 2 , posterior probability density functions (pdfs) of the parameters obtained using the SO and MO approaches are shown. The horizontal axes of the plots correspond with the initial 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 distributions or the CI95 compared with the width of the prior distributions (Fig. 2; Table 2) show a relatively good ability to identify the parameters θ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 measurements of only soil surface temperatures (SO) or on measurements of soil 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, the uncertainty intervals can also be plotted directly around the hydraulic 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 for the 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 the MO 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, and is somewhat narrower, i.e., from 0.072 to 0.106 cm3 cm−3 for the MO approach.

The optimization results of the van Genuchten–Mualem parameters were then used in the HYDRUS-1D forward simulation to compare 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 more significant. Figure 4 shows that the MO approach using measurements of both soil surface temperatures and water contents yielded a 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).

The synthetic experiment was designed to illustrate the usefulness of soil surface temperatures to derive the parameters of 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 water content data are used in addition to soil surface temperatures, the reference parameter set used to generate the synthetic measurement data was reproduced reasonably well using both the SO and the MO approaches.

Field Data Experiment

The measured field data were used to parameterize a model accounting for a soil profile with vertically uniform hydraulic properties. This assumption implies that an effective parameter set can be found representing a vertically heterogeneous soil profile and adequately predicting soil surface temperatures and soil 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 content dynamics properly. This is similar to the findings of Yang et 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 profile can be approximated assuming a homogeneous soil profile. Because the estimation of effective van Genuchten–Mualem hydraulic properties for a homogeneous soil profile 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), indicating that different soil hydraulic properties would have to be assigned to the two layers. We therefore subdivided the soil profile into a layer representing the harrowed upper part of the soil (0–10 cm) and a layer 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 profile 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 profile was divided into two horizons, with the parameters estimated in the first step assigned to 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 water contents 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 parameters of the other layer fixed, until the estimated hydraulic properties of both layers stopped changing between subsequent steps. The iterative process converged after three iterations. This 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 parameters of the soil layer of which the parameters were estimated was set to the 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 and water content measurements at the 7-cm depth (the MO approach), while keeping the soil hydraulic parameters for the deeper part of the soil profile constant. Except for the last step of the iteration process, we used the DREAM algorithm only until convergence, 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 profile because our additional modeling showed that soil surface temperatures did not contain sufficient information to infer soil water dynamics of the second, deeper layer of the soil profile.

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 parameter n 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 n cannot 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 use of 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, the parameter uncertainty was also considerably larger than that in the synthetic experiment. The CI95 in the synthetic experiment calculated from the SO approach for θr ranged from 0.037 to 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 concentration of 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.

Water retention curves resulting from the best van Genuchten–Mualem parameter sets and their 95% confidence intervals are shown in Fig. 8 for the SO and MO approaches. As expected, examining the posterior parameter distributions (Fig. 7), the uncertainty in the optimized retention curve was considerably reduced when soil water content data were considered together with soil surface temperature measurements in the optimization process. Figure 8 also compares the estimated water retention curves against those determined in the laboratory (Table 1). Both the SO- and MO-derived water retention curves for the surface soil layer deviated considerably from the retention curve determined in the laboratory. Both estimated retention curves are characteristic of a soil with a coarser texture (i.e., larger α parameter). The optimized retention curve for the topsoil layer also deviates considerably from the retention curve derived for the deeper soil layer, which shows more resemblance to that determined in the laboratory. This indicates larger pore sizes in the topsoil layer than what would normally be observed for a soil with this texture. The bulk density of the topsoil was decreased by harrowing at the beginning of the experiment, which may explain the effect on α. Despite the larger uncertainty of the SO-estimated hydraulic parameters, it is still possible to infer that the hydraulic properties of the disturbed topsoil layer are considerably different from the hydraulic properties determined on undisturbed soil samples.

Due to the lower bulk density of the harrowed upper soil layer than the untilled deeper layer, the saturated hydraulic conductivity and the saturated water content were expected to be larger for the upper soil layer; however, fitted Ks and θs in the lower soil layer were higher than in the upper soil layer (Table 4). This may be due to a lack of information in the wet part of the water retention curve, especially for the upper layer of the soil (measurements at 7 cm) where measured water contents did not exceed values around 0.24 cm3 cm−3. As a consequence, the estimated saturated hydraulic conductivity and water content are an extrapolation of hydraulic properties that were derived for drier soil conditions. For aggregated or tilled soils with a considerable interaggregate porosity, this extrapolation may lead to a considerably underestimation of the real saturated hydraulic conductivity and water content. A measurement period covering a wider range of dry and wet soil conditions may give better results concerning the saturated hydraulic conductivity and the saturated water content.

Figure 9 shows the water contents and surface temperatures simulated using parameters determined by both SO and MO approaches. Both approaches show a good agreement between measured and modeled water contents and soil surface temperatures. The RMSE for the water content at the 7-cm depth is 0.012 cm3 cm−3 for the SO approach and 0.011 cm3 cm−3 for the MO approach. At the 15-cm depth, the RMSE calculated from the SO approach is 0.0046 cm3 cm−3 and 0.0037 cm3 cm−3 for the MO approach. The soil water dynamics at the 7-cm depth following the rain 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.

Figure 9 also shows surface temperatures simulated using soil hydraulic parameters determined in the laboratory. Simulated surface temperatures significantly underestimated peak measured temperatures, which shows 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 that were obtained from the SO and MO approaches were similar, the reduction 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 hydraulic parameters, we determined whether soil temperatures measured at depths of 3 and 6 cm could be simulated properly with the parameterized model. Because both best parameter sets (from the SO and MO approaches) predicted nearly identical temperatures at both depths, only soil temperatures simulated using the MO parameters are plotted in Fig. 11 . It shows that modeled 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.

Summary and Conclusions

Estimating soil hydraulic properties using inverse modeling of the informational 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 estimated soil 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 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 profile that could be used to describe field-measured soil temperatures and water contents at different depths was not successful. This shows that models that try to describe 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 profile, it was found that the hydraulic 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 effect of 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 contain enough information to infer the hydraulic properties of the deeper part of the soil profile, 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 obtained water retention curves for the deeper soil layer corresponded better with those derived on laboratory columns. This indicates that soil surface temperature measurements should be combined with information 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 be obtained 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 the variability 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 and corresponding water retention curves that were obtained using only soil surface temperatures or a combination of temperature and 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 applicability of a particular method but can only provide a first approximation about the actual performance of the method under field conditions. Considering a layered soil profile, both parameter sets estimated from field-measured data (using MO and SO approaches) led to similar 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 be applied 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 of soil 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 subject of further research.

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, for providing the source code of the global optimizer DREAM.