Interactive comment on “ Estimate soil moisture using trapezoidal relationship between remotely sensed land surface temperature and vegetation index ”

General Comment: The paper describes the estimation of soil moisture using the trapezoidal method and a numerical, energy balance-based assessment of the vertices of the trapezoid. The paper is sufficiently original to warrant publication in HESS. The paper is fairly well written, but it would improve the readability of the paper even more if greater attention was provided to improving language fluency. Before the paper can be considered finalised, I would like the authors to address the following issues. 1. The method’s dependence on ground-based data is one significant limitation of the method, restricting its use in areas where ground-based meteorological data (i.e., air temperature, wind speed, relative humidity) are available. This method could not be


Introduction
In the 1980's, it was found that, land surface temperature (T s ) and the fraction of vegetation cover, which is represented by Vegetation Indices (VI), e.g., Normalized Difference Vegetation Index (NDVI), typically show a strong negative relationship (e.g., Goward et al., 1985;Nemani and Running, 1989).Such a relationship has been widely used to investigate the moisture condition of land surfaces.Several studies focused on the slope of the T s /NDVI curve for providing information on vegetation and moisture conditions at the surface (e.g., Smith and Choudhury, 1991;Nemani et al., 1993).Their approach was later extended to use the information in the Correspondence to: W. Wang (w.wang@126.com)T s /VI scatter-plot space, whose envelope is considered to be in either a triangular shape (e.g., Price, 1990;Carlson et al., 1994), or a trapezoid shape (e.g., Moran et al., 1994).
The idea of a triangle in the T s /VI space has been used to develop the so called "triangle method", and has been applied by many researchers for estimating soil moisture and evapotranspiration (see Carlson, 2007).The central assumption of the triangle method is that, given a large number of pixels reflecting a full range of soil surface wetness and fractional vegetation cover, sharp boundaries (edges) in the data scatter plot reflect real physical limits: i.e., bare soil, 100 % vegetation cover, and lower and upper limits of the surface soil water content, e.g., completely dry or wet (field capacity), respectively.The dry and wet edges ultimately intersect at a (truncated) point at full vegetation cover.Then, based on the triangle, the relative value of surface soil water content and the surface energy fluxes at each pixel can be defined in terms of its position within the triangle.The advantage of the triangle method is its independence of ancillary data.The approach, however, has difficulty in defining the dry and wet edge, especially the dry edge.Even with a large number of remotely sensed observations, the boundaries of the triangle space are still hard to establish, because on one hand, there are situations when T s ∼ VI points scatter in a narrow range such as during rainy season or in areas with a narrow VI range; on the other hand, the T s ∼ VI relationship is much more complicated at large scale than at local scale and may vary at different parts due to heterogeneity in land surface properties and atmospheric forcing.Furthermore, as the triangle space is established empirically, the soil moisture estimates according to such an empirical triangle based on an image at one time are hard to be compared with those at another time.
Published by Copernicus Publications on behalf of the European Geosciences Union.Moran et al. (1994) proposed the idea of vegetation index/temperature (VIT) trapezoid for describing the relationship between the surface temperature and air temperature difference (T s − T a ) vs. the fractional vegetation cover (V C ), and developed the Water Deficit Index (WDI) for evaluating evapotranspiration rates of both full-cover and partially vegetated sites.However, although the idea of VIT trapezoid is well accepted, very few applications were found in the literature based on the idea of VIT trapezoid for estimating soil moisture, partly because of the difficulty in calculating T s − T a , partly because of the limitation of meteorological data requirements.In the present paper, we aim to simplify the idea of VIT trapezoid to T s ∼ VI trapezoid, and propose an iterative algorithm for quantifying the shape of the T s ∼ VI trapezoid, then estimate soil moisture based on the T s ∼ VI trapezoid.The method of establishing the T s ∼ VI trapezoid will be described in detail in Sect. 2. Then the method will be applied to the Walnut Gulch Experimental Watershed in Arizona, USA, for which, the data used and data pre-processing will be described in Sects.3 and 4, and the results will be presented in Sect. 5. Finally, some conclusions will be drawn in Sect.6.
2 Trapezoid method 2.1 The concept of (T s − T a ) ∼ V C trapezoid Idso et al. (1981) and Jackson et al. (1981) proposed the CWSI (Crop Water Stress Index) for detecting plant water stress based on the difference between canopy and air temperature.It is designed for full-cover vegetated areas and bare soils at local and regional scales.To overcome the difficulty of measuring foliage temperature in partially vegetated fields, Moran et al. (1994) proposed to use the shape of trapezoid to depict the relationship between the surface temperature and air temperature difference (T s −T a ) vs. the fractional vegetation cover (V C , ranging from 0 for bare soil to 1 for full-cover vegetation) (Fig. 1), so as to combine spectral vegetation indices with composite surface temperature measurements to allow application of the CWSI theory to partially vegetated fields without a priori knowledge of the percent vegetation cover.
Based on the trapezoid assumption and the CWSI theory, Moran et al. (1994) introduced the Water Deficit Index (WDI) for evaluating field evapotranspiration rates and relative field water deficit for both full-cover and partially vegetated sites.For a given pixel with measured surface temperature and air temperature difference, i.e., (T s − T a ) r , WDI is defined as: where T a is air temperature; T s is surface temperature; the subscripts min, max, and r refer to minimum, maximum, and measuring foliage temperature in partially vegetated fields, Moran et al. (1 shape of trapezoid to depict the relationship between the surface tempera difference (T s -T a ) vs. the fractional vegetation cover (V C , ranging from full-cover vegetation) (Fig. 1), so as to combine spectral vegetation indice temperature measurements to allow application of the CWSI theory to without a priori knowledge of the percent vegetation cover.Based on the trapezoid assumption and the CWSI theory, Moran et al. (19 Deficit Index (WDI) for evaluating field evapotranspiration rates and relati both full-cover and partially vegetated sites.For a given pixel with measured air temperature difference, i.e., (T s -T a ) r , WDI is defined as: where T a is air temperature; T s is surface temperature; the subscripts min, ma maximum, and measured values, respectively; and the minimum and maxim interpolated linearly on the cold edge and warm edge of the (T s -T a )~V C trap value of the pixel.Graphically, WDI is equal to the ratio of distances AC/AB The theoretical basis of (T s −T a ) ∼ V C trapezoid is the energy balance equation, i.e., where, R n is the net radiant heat flux density (W m −2 ), G is the soil heat flux density (W m −2 ), H is the sensible heat flux density (W m −2 ), and λE is the latent heat flux to the air (W m −2 ) and λ the heat of vaporization of water (kJ kg −1 ).
In their simplest forms, H can be estimated with a bulk transfer equation written in the form (Monteith and Unsworth, 2008): and λE can be expressed as It should be noticed that, the usage of surface temperature T s , instead of aerodynamic temperature T aero , in Eq. ( 3) may lead to some errors for calculating H because the differences between T aero and T s could range from 2 • C over uniform vegetation cover to 10 • C for partially vegetated areas (Kustas and Norman, 1996).To solve this problem, the traditional approach consists in modifying the aerodynamic resistance r a by adding an extra resistance (e.g., Kustas et al., 1989), commonly expressed as a function of the friction velocity u* and the dimensionless bulk parameter B −1 , which is reflected in Eq. ( 10).
Then, combining Eqs. ( 2), (3), and (4), we obtain the equation for temperature difference between air and land surface: As suggested by Moran et al. (1994), for the (T s − T a ) ∼ V C trapezoid, its four vertices correspond to (1) well-watered full-cover vegetation, (2) water-stressed full-cover vegetation, (3) saturated bare soil, and (4) dry bare soil.Using the energy balance equations, Moran et al. (1994) computed the values of the four vertices of the trapezoid as the following: 1.For full-covered and well-watered vegetation (Point 1) where r cm is the minimum canopy resistance, i.e., canopy resistance at potential evapotranspiration.
2. For full-covered vegetation with no available water (Point 2) where r cx , is the canopy resistance associated with nearly complete stomatal closure.
3. For saturated bare soil (Point 3), where canopy resistance r c = 0, we have 4. For dry bare soil (Point 4), where r c = ∞ (analogous to complete stomatal closure), and λE = 0, we have The (T s − T a ) ∼ V C trapezoid considers that relationship between (T s − T a ) and V C .Now we think about the issue in another way that, with a given value of T a , how T s is related with V C .To analyze this T s ∼ V C relationship, we use Eqs.(6) ∼ (9) to calculate the T s for the four extreme cases (or trapezoid vertices) by moving T a in Eqs.(6) ∼ (9) to the right side of the equations.At the same time, V C is replaced by a vegetation index (VI).So that, we modify the structure of the trapezoid, obtaining a simplified T s ∼ VI trapezoid with the horizontal axis as the VI, and the vertical axis as T s .We therefore refer the algorithm proposed here to as T s -VI trapezoid method.
To obtain the values of T s with Eq. ( 6) ∼ (9), we need to know r a , r c (including r cm and r cx ), R n , G for the four vertices separately, as shown in the following section.

Calculation of the components in the formula for
four vertices of T s ∼ VI trapezoid

Aerodynamic resistance (r a )
The water vapor aerodynamic resistance r a (s m −1 ) can be estimated with the following equation (Thom, 1975): where d is displacement height (m), given by d = 0.667 h, and h is the height of vegetation (Garratt, 1992), which should be given as an input.z 0m is the roughness lengths for momentum (m), given by z 0m = h/8 (Garratt, 1992).For bare soil surface, z 0m is commonly taken to be 0.01 m (Shuttleworth and Wallace, 1985).
z 0h is the roughness lengths for heat (m), given by Here, kB −1 is a dimensionless parameter.Kustas et al. (1989) showed that kB −1 is a linear function of the product of u z and T s − T a , given by where S kB is an empirical coefficient, which varies somewhere between 0.05 and 0.25.
k is the von Karman constant (k = 0.41); ψ h and ψ m are the stability corrections for heat and momentum transfer (unitless).ψ h and ψ m are calculated differently depending on the atmospheric stability, which could be indicated by the Monin-Obukhov length L (Monin and Obukhov, 1954), given by where g = 9.8 m s −2 , ρ is the air density (kg m −3 ), C p the air specific heat at constant pressure (1004 J kg −1 K −1 ), u * is the friction velocity defined by u * = ku z ln[z−d/z 0m ]−ψ m .For neutral conditions (L = 0), ψ h = ψ m = 0.For stable situations (L > 0), expressions of ψ h and ψ m are (Webb, 1970): For unstable conditions (L < 0), expressions of ψ h and ψ m are (Paulson, 1970):

Net radiant heat flux density (R n )
Net radiation is defined as the difference between the incoming and outgoing radiation fluxes including both long-and shortwave radiation at the surface of Earth.Net radiant heat flux density (R n ) (W m −2 ) can be expressed as (Brutsaert, 2005): where α is surface shortwave albedo, which is derived from MODIS product MCD43A3; -R s is solar radiation, estimated jointly by solar constant, solar inclination angle, geographical location and time of year, atmospheric transmissivity, ground elevation, etc.The basic formula for estimating R s is (Zillman, 1972): where S 0 is the solar constant at the atmospheric top (1367 W m −2 ), θ the solar zenith angle, e 0 is the vapor pressure.In consideration of the effects of topography on the incident short-wave radiation (R s ), the solar zenith angle (θ) is corrected using Digital Elevation Model (DEM) data (Duffie and Beckman, 1991) with the following formula: cosθ = sin(δ) sin(ϕ) cos(s) − sin(δ) cos(ϕ) sin(s) cos(r) where φ is the latitude (positive in the Northern Hemisphere); s is the slope, and r is the slope orientation, both derived from DEM; δ is solar declination, and ω solar hour angle, given by where DOY is the day of year, and t is the time when the satellite Terra which carries MODIS pass over the region (∼10:30 a.m.).
There are MODIS land surface emissivity products (MOD11) for MODIS bands 31 and 32, i.e. ε 31 and ε 32 .But for extreme cases, i.e., bare soil and full vegetation coverage, the direct use of there products are not appropriate.Instead, according to Synder et al. (1998), we set ε s = 0.93 for bare soil, and ε s = 0.993 for full vegetation coverage.
In our algorithm, R n is not directly solved with the Eq. ( 16), because T s is considered as an unknown variable.Instead, we replace the term R n in Eqs. ( 6) ∼ (9) with the Eq. ( 16) respectively, so that we get four quartic equations for T s at four vertices separately.Then the quartic equations are solved with the iterative algorithm which is shown later in Sect.2.4 and Fig. 2, by doing so, all the values of T s for the four vertices are obtained.

Soil heat flux density G
G is normally considered to be linearly related to R n .Several studies have shown that the value of G/R n typically ranges between 0.4 for bare soil and 0.05 for full vegetation cover (Choudhury et al., 1987).Idso et al. (1975) conducted some experiments investigating the impacts of water content on the net radiation ∼ soil heat flux relationship over bare soil surface, and showed that G/R n ranges from 0.2 for wet bare soil to 0.5 for dry bare soil.

Canopy resistance (r c )
Canopy resistance (r c ), including r cm and r cx that refer to the minimum and maximum canopy resistances respectively, should be calculated for Point 1 and Point 2. According to Moran et al. (1994), r cm in Eq. ( 6) is calculated with r sm /LAI (LAI is the maximum possible leaf area index, r sm is minimum stomatal resistance).r cx in Eq. ( 7) is calculated with r sx /LAI (r sx is maximum stomatal resistance).
Values of minimum and maximum stomatal resistance (r sm and r sx , respectively) are published for many agricultural crops under a variety of atmospheric conditions.Moran et al. (1994) suggested that, if values are not available, reasonable values of r sm = 25 ∼ 100 s m −1 and r sx = 1000 ∼ 1500 s m −1 will not result in appreciable error, we set r sm = 100 s m −1 and r sx = 1500 s m −1 .As values of LAI of various vegetation types are mostly less than 8 (Scurlock et al., 2001), we set LAI = 8.Therefore, we have r cm = 12.5 s m −1 and r cx = 187.5 s m −1 .

Iterative procedure for calculating T s
Values of T s for the four vertices are obtained by an iterative procedure for each pixel.An initial value of r a is estimated

The Walnut Gulch Experimental Watershed
Data of the Walnut Gulch Experimental Watershed (WGEW) w WGEW is defined as the upper 148 km 2 of the Walnut Gulch drainag of the San Pedro catchment in southeastern Arizona (Fig. 3).It was d the United States Department of Agriculture (USDA) in the mi receives 250-500 mm of precipitation annually, with about two-third during a summer monsoon season.The potential evapotranspiratio annual rainfall.The topography can be described as gently rolling channels which are more pronounced at the eastern end of the catchm Soil types range from clays and silts to well-cemented boulder cong cm) soil textures being gravelly and sandy loams containing, on ave matter (Renard et al., 1993).The mixed grass-brush rangeland vege coverage.Grasses primarily cover the eastern half of the catchm bush-dominated.
Solve the quartic equations for T s by replacing r a in Eq. ( 6)~( 9) for each vertex Calculate T s -T a Calculate H with Eq. ( 5) Calculate kB -1 with Eq. ( 12) and L with Eq. ( 13) Correct the value of r a by updating ψ h and ψ m as in Eq. ( 14) or ( 15 by assuming neutral conditions, i.e., ψ h = ψ m = 0.With the initial r a , initial values of T s are obtained with Eq. ( 6) ∼ (9) for the four vertices.Then the iterative procedure is proceeded by iteratively changing H , kB −1 , r a , and in consequence, T s , until the value of T s is stable (i.e., the change of T s is less than 0.01 K, and the change of r a is less than 0.1 s m −1 ).Normally, it takes 5 to 10 iterations for T s to get stable.However, there are cases in which T s cannot reach a stable solution.In those cases, we use the T s value of the first iteration.While T s is derived, R n , G, H , and r a for each vertex are obtained as well.When T s is less than T a , we set H = 0.
The iterative procedure is conducted pixel by pixel, that is, the trapezoid is constructed separately for each pixel, and each trapezoid has its own values of T s .

MODIS data and ground observational data used
The moderate resolution imaging spectroradiometer (MODIS) instrument has been widely used for monitoring soil moisture because of its high spectral (36 bands) resolution, moderate spatial (250-1000 m) resolution, and various products for land surface properties.All standard MODIS data products are freely available at NASA Land Processes Distributed Active Archive Center (URL: https://lpdaac.usgs.gov/lpdaac/).MODIS products used in the present study are listed in Table 1.They are derived from the images of the MODIS sensor onboard Terra (~10:30 a.m.overpass).We selected MODIS data of ten cloud-free days approximately evenly distributed in the period from January to December in 2004.All the 1km resolution MODIS data are resampled to 500 m resolution with the nearest neighbor method.The reason of downscaling the 1km resolution data to the 500m resolution data is to keep more details with smaller pixel size.Meteorological data required here include air temperature T a , relative humidity μ, and wind velocity u, observed approximately at the time (11 a.m.) when the satellite Terra passes over the WGEW region.The T a , relative humidity μ, and wind velocity u, are observed at three sites.We take the average of the observations at 3 meteorological observation sites for μ and u.Observations of T a are pre-processed, which will be discussed in section 4.3.To evaluate the soil moisture estimation results, soil moisture observations at 16 sites and precipitation data at 87 sites are used.The locations of these sites are plotted in Fig. 4. As some gauging sites are located on the edge of the watershed, to include the observations at these sites for evaluation, our study area is slightly larger than WGEW.We used soil moisture observation data in 10 dates in 2004 when cloud-free MODIS data are available.All the soil moisture data are observed at 5 cm below the surface.The statistics of the soil moisture observation we used are listed in  3 Case study area and data used

The Walnut Gulch Experimental Watershed (WGEW)
Data of the Walnut Gulch Experimental Watershed (WGEW) was used in the present study.The WGEW is defined as the upper 148 km 2 of the Walnut Gulch drainage basin in an alluvial fan portion of the San Pedro catchment in southeastern Arizona (Fig. 3).It was developed as a research facility by the United States Department of Agriculture (USDA) in the mid-1950s.This rangeland region receives 250-500 mm of precipitation annually, with about two-thirds of it as convective precipitation during a summer monsoon season.The potential evapotranspiration is approximately ten times the annual rainfall.The topography can be described as gently rolling hills incised by steep drainage channels which are more pronounced at the eastern end of the catchment near the Dragoon Mountains.Soil types range from clays and silts to well-cemented boulder conglomerates, with the surface (0-5 cm) soil textures being gravelly and sandy loams containing, on average, 30 % rock and little organic matter (Renard et al., 1993).The mixed grass-brush rangeland vegetation ranges from 20 to 60 % in coverage.Grasses primarily cover the eastern half of the catchment, while the western half is bush-dominated.

MODIS data and ground observational data used
The moderate resolution imaging spectroradiometer (MODIS) instrument has been widely used for monitoring soil moisture because of its high spectral (36 bands) from the images of the MODIS sensor onboard Terra (∼10:30 a.m.overpass).We selected MODIS data of ten cloud-free days approximately evenly distributed in the period from January to December in 2004.All the 1 km resolution MODIS data are resampled to 500 m resolution with the nearest neighbor method.The reason of downscaling the 1 km resolution data to the 500 m resolution data is to keep more details with smaller pixel size.Meteorological data required here include air temperature T a , relative humidity µ, and wind velocity u, observed approximately at the time (11:00 a.m.) when the satellite Terra passes over the WGEW region.The T a , relative humidity µ, and wind velocity u, are observed at three sites.We take the average of the observations at 3 meteorological observation sites for µ and u.Observations of T a are preprocessed, which will be discussed in Sect.4.3.To evaluate the soil moisture estimation results, soil moisture observations at 16 sites and precipitation data at 87 sites are used.The locations of these sites are plotted in Fig. 4. As some gauging sites are located on the edge of the watershed, to include the observations at these sites for evaluation, our study area is slightly larger than WGEW.We used soil moisture observation data in 10 dates in 2004 when cloud-free MODIS data are available.All the soil moisture data are observed at Hydrol.Earth Syst.Sci., 15,[1699][1700][1701][1702][1703][1704][1705][1706][1707][1708][1709][1710][1711][1712]2011 www.hydrol-earth-syst-sci.net/15/1699/2011/ www.hydrol-earth-syst-sci.net/15/1699/2011/ Hydrol.Earth Syst.Sci., 15, 1699-1712, 2011 and EVI = 2.5 (ρ NIR − ρ red ) (ρ NIR + 6.0 ρ red − 7.5 ρ blue + 1) , where ρ blue , ρ red , and ρ NIR represent reflectance at the blue (0.45-0.52 µm), red (0.6-0.7 µm), and Near-Infrared (NIR) (0.7-1.1 µm) wavelengths, respectively.In the present study, EVI is used due to its less sensitivity to background reflectance variations (Liu and Huete, 1995).The MODIS EVI values varied from 0.07 to ∼0.7 for major land cover types from hyperarid deserts to dense forests at 1-km resolution (Huete et al., 2011).
The MODIS VI product attempts to retrieve cloud-free, near-nadir, top-of-canopy greenness at 16-day interval.However, due to the global nature of the algorithm, some problems and uncertainties persist at local scales, mostly associated with residual clouds, shadows, aerosols, atmospheric correction performance, and view sun angle geometries, resulting in nonbiological artifacts and noise in the VI values (Huete et al., 2011).Therefore many researchers have tried to denoise the MODIS VI data.Jennifer and McDermid (2009) compared six NDVI time series noise-reduction techniques, and found that the asymmetric Gaussian, Double logistic, and 4253H twice filter perform very well in general.As the EVI tends to have less negatively-biased noise and more erroneous spikes than the NDVI, in which case noise-reduction techniques maintaining the upper envelope of values such as the double logistic and asymmetric Gaussian function fitting techniques may not be the most effective choice (Jennifer and McDermid, 2009) whereas 4253H twice filter has the ability of eliminating spurious drops and spikes (Velleman, 1980), we tried to apply the 4253H twice filter to reduce the noise in EVI time series.Figure 5 shows the denoising effect with 4253H twice filter, which applies a series of running medians of varying temporal window size and a weighted average filter (e.g., Hanning filter), with re-roughing, to the EVI time series at a randomly selected pixel in WGEW over about one year.From Fig. 5 we see that both low values and high values are smoothed with 4253H twice filter.
While it works with the noise-reduction techniques, it is possible that the application of these techniques might be causing a lost of valid information because the fluctuation of EVI may be related to soil moisture changes.For instance, Wang et al. (2007) found that, in growing seasons, the lag time for NDVI to respond to soil moisture change is about 5 days or less at the semi-arid sites, and 10 days at the humid site.Because MODIS vegetation index product is 16-day composite data and our study area in an arid zone, the effects of soil moisture change should be reflected in the EVI products.Therefore, we decide to skip the denoising procedure, using the original EVI product instead.
series of running medians of varying temporal window size and a weighted ave Hanning filter), with re-roughing, to the EVI time series at a randomly selected pixel about one year.From Fig. 6 we see that both low values and high values are smoot twice filter.While it works with the noise-reduction techniques, it is possible that the app techniques might be causing a lost of valid information because the fluctuation of EV to soil moisture changes.For instance, Wang et al. (2007) found that, in growing seas for NDVI to respond to soil moisture change is about 5 days or less at the semi-arid s at the humid site.Because MODIS vegetation index product is 16-day composite da area in an arid zone, the effects of soil moisture change should be reflected in th Therefore, we decide to skip the denoising procedure, using the original EVI product i

Topographic correction of air temperature
With methods of estimating soil moisture using thermal satellite images, often b temperature and ground-based air temperature observations are needed.When applyi to mountainous regions, terrain effects have to be taken into account becaus significantly affect both land surface temperature and air temperature.To avoid the pr sloping terrain, some authors just eliminated those pixels in mountainous part (e.g 1994), while in some other cases, land surface temperature was corrected (e.g., Hassan the present study, we go the opposite way, i.e., instead of correcting land surface correct the air temperature. To make a successful air temperature interpolation, many factors should be tak such as the elevation difference between a pixel and monitoring stations, temperature geometric characteristics (slope, aspect) of each pixel cell, and vegetation coverag

Topographic correction of air temperature
With methods of estimating soil moisture using thermal satellite images, often both land surface temperature and ground-based air temperature observations are needed.When applying such methods to mountainous regions, terrain effects have to be taken into account because terrain would significantly affect both land surface temperature and air temperature.To avoid the problem of steeply sloping terrain, some authors just eliminated those pixels in mountainous part (e.g., Carlson et al., 1994), while in some other cases, land surface temperature was corrected (e.g., Hassan et al., 2007).In the present study, we go the opposite way, i.e., instead of correcting land surface temperature, we correct the air temperature.
To make a successful air temperature interpolation, many factors should be taken into account, such as the elevation difference between a pixel and monitoring stations, temperature vertical gradient, geometric characteristics (slope, aspect) of each pixel cell, and vegetation coverage.Moore et al. (1993) proposed a specific algorithm to calculate daytime temperature at different altitudes within a valley.Based on that, Bellasio et al. (2005) proposed a simplified equation in the form of where T i is the unknown atmospheric temperature (K) at a z i altitude (m), T b is the measured atmospheric temperature (K) at a z b altitude (m), β is the vertical temperature gradient (K m −1 ), C is a constant, LAI max and LAI i are, respectively, maximum leaf area index (LAI) and its value at z i , and S i is the ratio between direct shortwave radiation on the actual surface (with its slope and aspect) and direct shortwave radiation on a horizontal free surface.
The above equation did not consider the impacts of wind.But according to the research of McCutchan and Fox (1986), for their study area (an isolated, conical mountain with elevation ranging from 2743 to 3324 m), wind speeds greater Hydrol.Earth Syst.Sci., 15, 1699Sci., 15, -1712Sci., 15, , 2011 www.hydrol-earth-syst-sci.net/15/1699/2011/ 12 four vertices of the trapezoids constructed for all the pixels of the WGEW region in four days in four seasons in Fig. 7.All the estimated T s at each point are plotted in the form of box-and-whisker plot.The data points (solid dots) of T s vs. EVI are also plotted in the map.From Fig. 7, we see that the constructed trapezoids well characterize the T s ~EVI space, and basically all the T s ~EVI data points are set in the envelope of the trapezoids.

Calculation of WDI
Based on the constructed T s ~VI trapezoid for each pixel, using the MODIS land surface temperature and EVI data, we calculate the WDI for each pixel p, than 5 m s −1 negate any slope, elevation or aspect effects that are present at low wind speed.We approximate this wind effect with a coefficient e −u/3 (u is the wind speed, m s −1 ), in consequence, we obtain a modified equation of Eq. ( 17) as Therefore, when there are air temperature observations at several sites, we can conduct air temperature correction in the following three steps: 1. Correct the observations to a flat plane at a base level All the temperature data are corrected to a flat plane at a base level (the lowest elevation z 0 of the observation sites), considering the effects of not only the elevation difference, but also the effects of wind, slope, and aspect.This is basically a reverse correction of Eq. ( 18), i.e., where a,b is the temperature observation corrected to the base level at site i, β is the temperature lapse rate (K m −1 ), z i is the elevation of site i, and z 0 is the elevation of the base level (m).
2. Interpolate temperature for each pixel p using observations on the flat plane at the base level Use the corrected air temperature observations a,b to interpolate the air temperature for all pixels with a spatial interpolation method (e.g., the inverse distance weighting interpolation method) to get interpolated air temperature T p a,I for each pixel p on the flat plane at the base level.
3. Topographic correction for each pixel p to its real position using Eq. ( 18), where T b is replaced by T p a,I .Bellasio et al. (2005) suggested to set LAI max = 10.As values of LAI of various vegetation types are mostly less than 8 (Scurlock et al., 2001), here for the arid zone WGEW we set LAI max = 8, C = 2 and β = 0.0065 K m −1 .β = 0.0065 K m −1 represents a long term average.It would be better to investigate the air temperature difference at sites of difference elevation to derive an accurate dynamic short-term estimate of β for a specific area.But the value of 0.0065 could be a reasonable estimate when not enough observations with different elevations are available to make a more accurate estimate.From the scatter plot of WDI vs. soil moisture observation in Fig. 9, we see that from the perspective of a whole year, WDI estimates derived with the T s ~VI trapezoid method has a negative correlation (correlation coefficient r=-0.673) with surface soil moisture, which indicates that WDI estimates can be used to detect the temporal variation in soil moisture.Especially on the scale of the watershed, the average WDI is strongly negatively related (correlation coefficient r=-0.924) to the average soil moisture observation, as shown in Fig. 10.Similar phenomena have been observed by some other researchers as well.For instance, Pellenq et al. (2003) noticed that the point-to-point comparison between observations and simulations shows a poor correlation, but a good correlation is obtained when averaging the simulated and observed soil moisture over a length of 100 m.
The comparison between the WDI estimates with ground observations of each date (Table 3) shows that there is basically no correlation between WDI estimates and surface soil moisture observations.This is partly because of the scale effect, i.e., point soil moisture observations are essentially different from pixel averaged soil moisture estimates due to sub pixel variability, partly due to the lower spatial variability than the temporal variability in soil moisture that makes it more difficult in using WDI to detect the spatial variation than to detect the temporal variation.Comparing the statistics of soil moisture observations in Table 2, we see that the average Coefficient of Variation (CV) for soil moisture observations at 16 sites in the 10 DOYs is 0.732 (ranging from 0.429 to 1.187), much larger than the average CV, 0.510, for observed soil moisture in any given date (ranging from 0.325 to 0.677).In consequence, we can use WDI to detect the temporal variation in soil moisture, but it seems that it would be hard to detect spatial variation in a day, especially for a small watershed with where T s is surface temperature obtained from MODIS; the subscripts min and max refer to minimum and maximum values; and the minimum and maximum values of T s are interpolated linearly on the dry edge and wet edge of the T s ∼ VI trapezoid based on the T s values calculated at four vertices for the specific VI value of the pixel p.
The WDI maps in ten DOYs for WGEW are illustrated in Fig. 7.

Comparison with soil moisture observation and precipitation
Using the surface soil moisture observations at 16 sites in 10 dates, we evaluate WDI estimates in several ways: (1) compared WDI estimates with ground observations of each site in 10 dates (Fig. 8); (2) compare the average of WDI estimates with the average ground observations of all sites in 10 dates (Fig. 9); (3) compare the WDI estimates with ground observations of all sites in each date (Table 3).
From the scatter plot of WDI vs. soil moisture observation in Fig. 8, we see that from the perspective of a whole year, WDI estimates derived with the T s ∼ VI trapezoid method has a negative correlation (correlation coefficient r = −0.673)with surface soil moisture, which indicates that WDI estimates can be used to detect the temporal variation in soil moisture.Especially on the scale of the watershed, the average WDI is strongly negatively related (correlation coefficient r = −0.924) to the average soil moisture observation, as shown in Fig. 9. Similar phenomena have been observed by some other researchers as well.For instance, Pellenq et al. (2003) noticed that the point-to-point comparison between observations and simulations shows a poor correlation, coefficient.The p-value is the t-test result for testing significance null hypothesis is: correlation coefficient = 0.A p-value of les correlation at a 0.05 significance level.)From the scatter plot of WDI vs. soil moisture observation perspective of a whole year, WDI estimates derived with the T s ~V correlation (correlation coefficient r=-0.673) with surface soil mo estimates can be used to detect the temporal variation in soil moist watershed, the average WDI is strongly negatively related (corre average soil moisture observation, as shown in Fig. 10.Similar ph some other researchers as well.For instance, Pellenq et al. (200  comparison between observations and simulations shows a poor co obtained when averaging the simulated and observed soil moisture o The comparison between the WDI estimates with ground ob shows that there is basically no correlation between WDI esti observations.This is partly because of the scale effect, i.e., poi essentially different from pixel averaged soil moisture estimates due to the lower spatial variability than the temporal variability in difficult in using WDI to detect the spatial variation than to detect the statistics of soil moisture observations in Table 2, we see that th (CV) for soil moisture observations at 16 sites in the 10 DOYs is 0. much larger than the average CV, 0.510, for observed soil moistur 0.325 to 0.677).In consequence, we can use WDI to detect the temp it seems that it would be hard to detect spatial variation in a day, esp but a good correlation is obtained when averaging the simulated and observed soil moisture over a length of 100 m.
The comparison between the WDI estimates with ground observations of each date (Table 3) shows that there is basically no correlation between WDI estimates and surface soil moisture observations.This is partly because of the scale effect, i.e., point soil moisture observations are essentially different from pixel averaged soil moisture estimates due to sub pixel variability, partly due to the lower spatial variability than the temporal variability in soil moisture that makes it more difficult in using WDI to detect the spatial variation than to detect the temporal variation.Comparing the statistics of soil moisture observations in Table 2, we see that the average Coefficient of Variation (CV) for soil moisture observations at 16 sites in the 10 DOYs is 0.732 (ranging from 0.429 to 1.187), much larger than the average CV, 0.510, for observed soil moisture in any given date (ranging from 0.325 to 0.677).In consequence, we can use WDI to detect the temporal variation in soil moisture, but it seems that it would be hard to detect spatial variation in a day, especially for a small watershed with low spatial soil moisture variability at the 500 m pixel scale.
Despite of the poor performance for characterizing the spatial variability of soil moisture with WDI, by a visual inspection of the WDI maps of the WGEW region of the 10 dates in Fig. 7, we can still see a clear spatial pattern of soil moisture distribution, which indicates that, to some extent, soil moisture variability could be depicted by WDI maps.
We analyzed the impacts of precipitation on soil moisture by calculating the correlation between WDI and Antecedent Precipitation (AP) of different number of days, and between soil moisture observation and AP of different number of days.The results are illustrated in Fig. 10, which show that WDI and soil moisture observation have similar levels of correlation with AP (one is positive, another is negative), and the maximum correlation occurs when approximately 10-day AP is taken into account.The scatter plot is shown in Fig. 11.The result indicates that, as expected, the temporal variation www.hydrol-earth-syst-sci.net/15/1699/2011/ Hydrol.Earth Syst.Sci., 15, 1699-1712, 2011

15
Note: Meanings of r and p-value are the same as in Fig. 9.
Despite of the poor performance for characterizing the spatial variability of soil moisture with WDI, by a visual inspection of the WDI maps of the WGEW region of the 10 dates in Fig. 8, we can still see a clear spatial pattern of soil moisture distribution, which indicates that, to some extent, soil moisture variability could be depicted by WDI maps.
We analyzed the impacts of precipitation on soil moisture by calculating the correlation between WDI and Antecedent Precipitation (AP) of different number of days, and between soil moisture observation and AP of different number of days.The results are illustrated in Fig. 11, which show that WDI and soil moisture observation have similar levels of correlation with AP (one is positive, another is negative), and the maximum correlation occurs when approximately 10-day AP is taken into account.The scatter plot is shown in Fig. 12.The result indicates that, as expected, the temporal variation of soil moisture (either reflected by ground observations, or by WDI estimates) is significantly dominated by precipitation process.Despite of the poor performance for characterizing the spatial variability of soil moisture with WDI, by a visual inspection of the WDI maps of the WGEW region of the 10 dates in Fig. 8, we can still see a clear spatial pattern of soil moisture distribution, which indicates that, to some extent, soil moisture variability could be depicted by WDI maps.
We analyzed the impacts of precipitation on soil moisture by calculating the correlation between WDI and Antecedent Precipitation (AP) of different number of days, and between soil moisture observation and AP of different number of days.The results are illustrated in Fig. 11, which show that WDI and soil moisture observation have similar levels of correlation with AP (one is positive, another is negative), and the maximum correlation occurs when approximately 10-day AP is taken into account.The scatter plot is shown in Fig. 12.The result indicates that, as expected, the temporal variation of soil moisture (either reflected by ground observations, or by WDI estimates) is significantly dominated by precipitation process. of soil moisture (either reflected by ground observations, or by WDI estimates) is significantly dominated by precipitation process.

Conclusions
Considerable efforts have been put on using the relationship between surface temperature and vegetation index (T s ∼ VI) to estimate surface soil moisture in the last two decades.
Since the publication of the paper by Moran et al. (1994), where they defined the trapezoidal relationship between the surface temperature and air temperature difference (T s − T a ) vs. the fractional vegetation cover (V C ), the shape of trapezoid has been commonly accepted as one of the ways of characterizing the T s ∼ VI relationship.However, in the algorithm proposed by Moran et al. (1994), when they calculated the value of T s − T a , no consideration was taken about the feedback effect of changes in T s − T a on some variables such as r a and R n which are used in the calculation of T s − T a .In addition, there is a problem of applying the method for areas with complex terrains and limited availability of ground meteorological observations.In the present study, we simplified the T s − T a versus V C trapezoid to the T s ∼ VI trapezoid for each pixel, and proposed an algorithm to iteratively update the values of quantities such as R n and r a so to keep the T s changing until it reaches a stable value.Then the Water Deficit Index (WDI) is calculated for each pixel based on the location of the data point of MODIS remotely sensed surface temperature versus Enhanced Vegetation Index (EVI) in the T s ∼ VI trapezoid.Using ground-based observations at the Walnut Gulch Experimental Watershed (WGEW) in Arizona, USA, the capability of using WDI to estimate soil moisture is evaluated by comparing it with soil moisture observations and antecedent precipitation.The result shows that, T s ∼ VI trapezoid based Hydrol.Earth Syst.Sci., 15, 1699Sci., 15, -1712Sci., 15, , 2011 www.hydrol-earth-syst-sci.net/15/1699/2011/ WDI can well capture temporal variation in surface soil moisture, but the capability of detecting spatial variation is poor for such a semi-arid region as WGEW.
The advantages of the T s ∼ VI trapezoid method include: (1) in comparison with the triangle method, which sets the dry and wet edges empirically, the trapezoid method determines the dry and wet edges based on a solid physical background, and hence, the WDI values for different periods are comparable, which makes WDI an appropriate index for monitoring the temporal change in soil moisture; (2) the method is applicable for even a small region with a narrow range of soil wetness and fractional vegetation cover, which is not possible with commonly used methods such as the triangle method or the SEBAL model (Bastiaanssen et al., 1998).
There are some limitations of the T s ∼ VI trapezoid method as well, including: (1) the method requires groundbased data, which restricts its use in areas where groundbased meteorological data (i.e., air temperature, wind speed, relative humidity) are poor or not available.For areas with poor to no data coverage, we have to either use limited observations to do interpolation, or take advantage of regional climate model such as the Weather Research and Forecasting (WRF) Model (http://www.wrf-model.org/).( 2) Some parameters, including those seem to be sensitive, such as S kB for calculating kB −1 , were determined empirically or by a trial and error procedure.In addition, the albedo values for four vertices should be different due to the difference of land surface, but here they are considered identical as that provided by the MODIS MCD43 product because of the difficulty in quantitatively describing the difference.These are topics for future research.

Fig. 1
Fig.1The hypothetical trapezoidal shape based on the relation between (T s -T a ) an vegetation cover (V C ).

Fig. 1 .
Fig. 1.The hypothetical trapezoidal shape based on the relation between (T s − T a ) and the fractional vegetation cover (V C ).

-
T s and T a are the surface temperature and air temperature (K), respectively; -C v is the volumetric heat capacity of air (1295.16J K −1 m −3 ); -VPD (vapor pressure deficit of the air) (hPa) is calculated as a difference between saturation vapour pressure Hydrol.Earth Syst.Sci., 15, 1699-1712, 2011 www.hydrol-earth-syst-sci.net/15/1699/2011/ e s and actual vapour pressure e a (hPa), given by (WMO, 2008) e S = 6.112 exp 17.62 T a T a + 243.12 (T a is the air temperature in • C, T a = T a − 273.15), e a = µe s , (µ is observed relative humidity) -is the slope of the curve of saturation water vapour pressure versus air temperature, calculated with (WMO, 2008) = 4098 • e s 237.3 + T a 2 γ the psychrometric constant (hPa k −1 ), given by (WMO, 2008) γ = 0.646 + 0.0006 T a r a the aerodynamic resistance to heat transfer (s m −1 ); r c the canopy resistance to vapor transport (s m −1 );

Fig. 2
Fig.2 Iterative procedure for calculating T s of the four verti

Fig. 2 .
Fig. 2. Iterative procedure for calculating T s of the four vertices of T s ∼ VI trapezoid.

Fig. 3
Fig. 3 Digital elevation model (DEM) of Walnut Gulch Experimental Watershed

_
Meteorological observation station ( Soil moisture observation site !Rain gage

Fig. 4
Fig.4 Locations of ground-based observation sites in WGEW

Fig. 6
Fig.6 Effects of EVI denoising preprocessing for a randomly selected pixel (Note: n horizontal axis indicate the 25 consecutive 16-day composite data sets starting from th 2003 to the first one in 2005)

Fig. 5 .
Fig. 5. Effects of EVI denoising preprocessing for a randomly selected pixel (Note: numbers of the horizontal axis indicate the 25 consecutive 16-day composite data sets starting from the last dataset in 2003 to the first one in 2005).
Fig.9WDI estimates vs. ground observations at 16 sites in 10 dates (Note: r is the correlation coefficient.The p-value is the t-test result for testing significance of the correlation coefficient.The null hypothesis is: correlation coefficient = 0.A p-value of less than 0.05 indicates significant correlation at a 0.05 significance level.)

Fig. 8 .
Fig. 8. WDI estimates vs. ground observations at 16 sites in 10 dates.(Note: r is the correlation coefficient.The p-value is the t-test result for testing significance of the correlation coefficient.The null hypothesis is: correlation coefficient = 0.A p-value of less than 0.05 indicates significant correlation at a 0.05 significance level.) Fig.10The average WDI estimates vs. the average ground obs r and p-value are the same as in Fig.9.)

Fig. 9 .
Fig. 9.The average WDI estimates vs. the average ground observations in 10 dates.(Meanings of r and p-value are the same as in Fig. 8.)

Fig. 12 Fig. 10 .
Fig.11 Correlation coefficient (r) between (a) soil moisture observation and AP of different number of days, and (b) WDI and AP of different number of days

Fig. 12 Fig. 11 .
Fig.11 Correlation coefficient (r) between (a) soil moisture observation and AP of different number of days, and (b) WDI and AP of different number of days

Table 2
. All the ground-based observational data, including meteorological observations, soil moisture observations and land cover data are obtained from the website of United States Department of Agriculture (USDA) Southwest Watershed Research Center (URL: http://www.tucson.ars.ag.gov/dap/).In addition, SRTM Digital Elevation Model (DEM) data (URL: http://srtm.csi.cgiar.org/)are used.
resolution, moderate spatial (250-1000 m) resolution, and various products for land surface properties.All standard MODIS data products are freely available at NASA Land Processes Distributed Active Archive Center (URL: https://lpdaac.usgs.gov/lpdaac/).MODIS products used in the present study are listed in Table 1.They are derived *MCD43A3 product is produced every 8 days, using data from the last 16 days.

Table 2 .
Statistics of surface soil moisture observations in 10 DOYs in 2004.SiteSoil moisture in 10 DOYs at each site DOY Soil moisture at 16 sites in each DOY Note: DOY is the Day Of Year.SD is the standard deviation; CV = SD/Mean.

Table 2
Statistics of surface soil moisture observations in 10 DOYs in 2004 Soil moisture in 10 DOYs at each site Soil moisture at 16 sites in each DOY

Table 3 .
Correlation between WDI estimates with surface soil moisture observations.
Note: Meanings of r and p-value are the same as in Fig.8.