Estimation of Surface Soil Moisture from Thermal Infrared Remote Sensing Using an Improved Trapezoid Method

Surface soil moisture (SM) plays a fundamental role in energy and water partitioning in the soil–plant–atmosphere continuum. A reliable and operational algorithm is much needed to retrieve regional surface SM at high spatial and temporal resolutions. Here, we provide an operational framework of estimating surface SM at fine spatial resolutions (using visible/thermal infrared images and concurrent meteorological data) based on a trapezoidal space defined by remotely sensed vegetation cover (Fc) and land surface temperature (LST). Theoretical solutions of the wet and dry edges were derived to achieve a more accurate and effective determination of the Fc/LST space. Subjectivity and OPEN ACCESS Remote Sens. 2015, 7 8251 uncertainty arising from visual examination of extreme boundaries can consequently be largely reduced. In addition, theoretical derivation of the extreme boundaries allows a per-pixel determination of the VI/LST space such that the assumption of uniform atmospheric forcing over the entire domain is no longer required. The developed approach was tested at the Tibetan Plateau Soil Moisture/Temperature Monitoring Network (SMTMN) site in central Tibet, China, from August 2010 to August 2011 using Moderate Resolution Imaging Spectroradiometer (MODIS) Terra images. Results indicate that the developed trapezoid model reproduced the spatial and temporal patterns of observed surface SM reasonably well, with showing a root-mean-square error of 0.06 m3·m−3 at the site level and 0.03 m3·m−3 at the regional scale. In addition, a case study on 2 September 2010 highlighted the importance of the theoretically calculated wet and dry edges, as they can effectively obviate subjectivity and uncertainties in determining the Fc/LST space arising from visual interpretation of satellite images. Compared with Land Surface Models (LSMs) in Global Land Data Assimilation System-1, the remote sensing-based trapezoid approach gave generally better surface SM estimates, whereas the LSMs showed systematic underestimation. Sensitivity analyses suggested that the trapezoid method is most sensitive to field capacity and temperature but less sensitive to other meteorological observations and parameters.


Introduction
Surface soil moisture (SM) plays an essential role in determining various land surface processes and the feedback between the Earth and the climate system [1].However, long-term and large scale in situ SM data are limited due to logistic and economic reasons.In addition, local-scale variations in soil properties, terrain, and vegetation cover result in a high spatial variability in surface SM, which makes it difficult to fully assess regional surface soil water conditions based on limited point measurements [2].
Remote sensing techniques provide a unique opportunity to capture land surface information over large geographic extents.Particularly, microwave remote sensing is often considered most effective in characterizing near surface SM from the space (e.g., [3][4][5]).It has been shown that the relationship between the L-band (frequency f = 1-2 GHz, wavelength λ = 30-15 cm) brightness temperature and the SM profile (up to 5 cm) is generally stronger than that at higher frequency [5].In addition, relatively longer wavelengths of microwave allow it to penetrate non-perceptible clouds and thus can furthest reduce the risk of losing data due to unfavorable climate conditions.On the other hand, lower frequencies (or longer wavelengths) result in coarser resolutions.For example, the Soil Moisture and Ocean Salinity (SMOS) satellite launched in November 2009 carries the low frequency L-band sensor, which provides surface SM retrievals every three days at 40 km resolution [6].The Advanced Microwave Scanning Radiometer-Earth Observation System (AMSR-E) on board the Aqua satellite provides a nearly daily global coverage but with a mean spatial resolution of ~56 km [7].Those coarse resolution microwave products offer valuable mean SM information over large areas, but lose detailed SM spatial variability because the signals are smoothed out at a coarser level.Moreover, numerous studies have shown that the microwave-based method may fail over surfaces with moderate to dense vegetation (e.g., [8]).
Although it is imperative to develop disaggregation methods to downscale microwave remote sensing SM estimates into finer resolutions [9], it is worthwhile to explore approaches to estimate SM moisture directly based on other types of remotely sensed information, i.e., visible/thermal infrared remote sensing images [10].Because the surface radiant temperature depends on surface soil water content and the distribution of vegetation, there should be a way for diagnosing SM by exploiting observations of land surface temperature (LST) in the thermal infrared band (8 to 13 microns) in combination with observations of surface greenness from the visible (380-760 nm) and near infrared spectrum (760 nm-1 micron).By interpreting the relationship between LST and the remotely sensed vegetation index (VI), a multitude of studies have shown that the envelope of data cloud (VI/LST) forms a triangle when a sufficiently large number of pixels are present with contaminated and outlier pixels being removed (see review by Carlson [11] for detail).In such a space, a higher VI value generally corresponds to a lower LST and a higher surface SM and vice versa.Furthermore, isolines of surface soil wetness were found in the VI/LST space [11], representing the same soil moisture availability.Once the boundaries of the triangle are defined, slopes of each soil wetness isoline can be determined by interpolation between the two boundaries that correspond to extreme dry and wet conditions, respectively [12][13][14].This so-called "triangle method" has been used extensively in estimating evapotranspiration (ET) (e.g., [12,15]) but is rarely applied to estimate surface SM.Whilst a few studies reported some promise in estimating surface SM based on the triangle method (e.g., [16,17]), none of them has transitioned to an operational framework.
There are, however, issues associated with the triangle method to be resolved.First, Moran et al. [18] claimed that the triangular space, in which the driest and wettest fully vegetated surfaces have the same surface temperature, does not account for the impact of water stress on canopy transpiration, and replaced the triangle with a trapezoid.In the trapezoidal space, the temperature of the extreme dry fully vegetated surface is higher than that of the wettest surface with full vegetation cover but is lower than extreme dry bare soil surface.Yang et al. [15] demonstrated that the trapezoid is more theoretically sound than using other shapes.Second, traditional approaches for determining the extreme boundaries of the triangle (or trapezoid) space are based on visual examination of the VI/LST data points (e.g., [19]), which implicitly requires that both conditions of extreme dry (lowest SM) and extreme wet points (highest SM) be present in a study scene to ensure the validity of the triangle (or trapezoid) method [11].This is often not the case in reality.In addition, substantial subjectivity results from the observed boundaries as well, because they can differ among operators for the same study area [12].Furthermore, visual examination of extreme boundaries would largely increase "labor work" and therefore severely reduces the operability of the method.Third, one key assumption of the triangle (or trapezoid) method is that among-pixel variation in surface energy flux is primarily controlled by that in water availability but is not affected by differences in atmospheric conditions [20].As a result, the triangle (or trapezoid) method is applicable for areas with relatively homogeneous meteorology.
In order to address the above-mentioned issues, this study provides an operational framework of estimating surface SM based on visible and thermal infrared information from satellite images and the trapezoid method.Rather than determining the trapezoidal boundaries from visual inspection, theoretical solutions were derived to achieve a more accurate and effective determination of the VI/LST space.
The theoretical derivation of the extreme boundaries also allows a per-pixel determination of the VI/LST space.Consequently, the assumption of uniform atmospheric forcing over the entire scene is no longer required.In addition to surface SM retrieval, this theoretically defined trapezoidal space may also benefit other applications, such as estimation of surface ET (e.g., [21]) and surface temperature (e.g., [15]).Here, we focus on SM, as the progress of SM studies is indeed far behind that of ET and temperature retrieval from visible/thermal infrared remote sensing.The developed approach will be validated with a well-designed SM observation network on the Tibet Plateau [22], and compared with SM estimates from land surface models [23].Meanwhile, the strength of theoretical boundaries over observed boundaries will be discussed.

Trapezoidal VI/LST Space
Theoretically, four critical points related to four extreme conditions define a trapezoid (Figure 1).Points A and B represent the driest surfaces without vegetation and with full vegetation coverage (Fc), respectively.As a result, points A and B constitute the warm edge of the trapezoid space.On the warm edge, the surface has the largest water stress for the entire range of Fc (from 0 to 1), and surface SM for those surfaces is assumed to be minimum (residual soil water content, θR).Accordingly, evaporation from soil surfaces and vegetation transpiration due to nearly complete stomatal closure on the theoretical warm edge are both taken to be zero.Though in reality, this situation would not occur frequently, the theoretical warm edge does provide a prescribed boundary condition given a certain meteorological field to confine ET/SM estimates for surfaces with generic temperatures and fractional vegetation cover.On the contrary, points C and D represent fully vegetated and bare soil surface without water stress, respectively, and segment CD is referred to as the cold edge.On the cold edge, the surface SM is at its maximum value (field capacity, θF) and ET occurs at its potential rates.Soil wetness isolines representing constant soil moisture availability were found to exist in the Fc/LST space based on inverse modelling of a soil-vegetation-atmospheric transfer model [16,17] (Figure 1).The slope of each isoline is derived by interpolating the slope of the warm edge and that of the cold edge [12][13][14], in terms of the temperature difference between the pixel and cold edge (a) and that between the pixel and warm edge (b).It is noted that slopes of the soil moisture isolines in the framework used in this study range from the lowest (negative) to zero (see Section 2.2).This is because theoretically the temperature of sunlit soil should not be lower than that of vegetation along the same soil moisture isoline.As a result, surface SM (θ, in volumetric water content (VWC)) for each pixel can be calculated from: where Tmin is the surface temperature for the cold edge.Ts,max and Tc,max are surface temperatures for the driest bare soil surface and driest fully vegetated surface, respectively.Note that the term a/(a + b) in Equation ( 1) represents surface soil moisture availability [11].

Determination of Boundaries
Theoretical solutions were proposed to avoid subjectivity and uncertainties in determining the trapezoidal space [12,21].For bare soil surfaces, the radiation budget and energy budget can be described as where Rn,s is the net radiation for bare soil surfaces (W m −2 ); Sd is the downward shortwave radiation (W m −2 ), which is calculated following [24]; αs is the bare surface albedo; σ is the Stefan-Boltzmann constant (5.67×10 −8 W m −2 K −4 ); εs is the bare surface emissivity (e.g., 0.95) [25]; εa is the atmospheric emissivity, which is a function of air pressure (Pa) and vapor pressure (ea) [26].Ta is the air temperature and Ts is the bare surface temperature (K); G, Hand LE are ground heat flux, sensible heat flux and latent heat flux, respectively (W m −2 ); ρ is the air density (kg m −3 ); Cp is the specific heat of air at constant pressure (J kg −1 K −1 ); ra,a is the aerodynamic resistance to heat transfer between Zom + d (Zom is the canopy roughness length for momentum transfer, and d is the zero displacement height) and the reference height (m s −1 ); ra,s is the aerodynamic resistance to heat flow in the boundary immediately above the soil surface (m s −1 ).Estimation of those aerodynamic resistance terms follows [27].
The first two terms of Taylor's formula of the upward longwave radiation (εsσTs 4 ) at Ta is written as [12] Substituting Equation ( 6) into Equation ( 4) and combine it with Equation ( 5), we can get an explicit expression of Ts, For extreme dry surfaces where there is the largest water stress, ET is completely suppressed (LEs = 0).As a result, the surface temperature for the driest bare surface (point A in Figure 1), Ts,max, is, where the ratio of G/Rn,sis taken to be 0.35 for bare soil surfaces [28].
In a similar vein, the temperature for the driest fully vegetated surface (point B in Figure 1), Tc,max, can be derived as, where αc is the albedo for fully vegetated surfaces; εc is the vegetation emissivity (0.98) [25]; ra,c is the aerodynamic resistance to heat transfer between canopy and the reference height (m s −1 ).The two albedo terms (αs and αc) in Equations ( 8) and ( 9) can be estimated either by extending the upper envelope of the Fc/α space intersecting with Fc = 0 and Fc = 1, respectively, or from the look-up table for each specific type of surface provided by [29].For the cold edge, the highest ET corresponds to the lowest sensible heat, and therefore the temperature gradient between the land surface and reference height would approach zero.As such, air temperature is taken as the surface temperature for the lower theoretical limiting edge (i.e., Ts,min = Tc,min = Ta).This simplification is an operational way to determine the theoretical lower limiting edge and obviates the need for vapor pressure deficit in the work of [18].However, this simplification may not be suited for surfaces with either extremely high vapor pressure deficit or strong advection effect [21].A flowchart of the trapezoidal model is shown in Figure 2.

Site Description and SM Measurement
In situ measurements of surface SM to validate the algorithm were obtained from a well-designed soil moisture observation network located around Naqu town, Tibet (31°N-32°N, 91.5°E-92.5°E),covering an area of ~100 km ×100 km with a mean elevation of 4650 m a.s.l.(the Tibetan Plateau Soil Moisture/Temperature Monitoring Network, SMTMN) (Figure 3) [22,30].The area is fairly smooth with rolling hills despite remarkable relief in some parts.The region has a typical sub-humid monsoon climate.Annual precipitation is about 400-500 mm, with ~80% occurring during the monsoon season (June to August).The main components of soil texture are sand and silt, and exhibit a large spatial heterogeneity [31].Organic carbon content was observed to be high within the top soil layer and declines dramatically with depth [31].Alpine grasslands cover about 93.6% of the land surface, with a mean vegetation height of several centimeters [22].Soil thawing and freezing are typical phenomena in the region, which generally occur around May and November, respectively [22].Soil moisture at four depths (0-5 cm, 10 cm, 20 cm and 40 cm) were measured at 30 sites along four branches of the national/provincial roads (Figure 3) from 01/08/2010/ to 31/08/2011 with capacitance probes (http://www.decagon.com/)with an accuracy of ±3% VMC and resolution of 0.1% VMC.Another 20 sites were added into the measurement since July 2011.The measurement time interval was 30 min, and each record reflected the average state of SM over the past half-hour.Sensor calibration was performed using the gravimetric method in laboratory.Detailed description of the study region and SM measurement can be found in [22,32].

Remote Sensing Data
Moderate Resolution Imaging Spectroradiometer (MODIS) data were used in this study because of its high temporal resolution (daily) and accessible spatial resolution (250-1000 m).These data were provided by the NASA Data Center [33].Daily LST was obtained from the global daily land surface temperature and emissivity product (MOD11A1), and broadband surface albedo was calculated based on narrowband surface reflectance from dataset MOD09GA and the algorithm proposed by [34].The normalized difference vegetation index (NDVI) was derived from red and near-infrared bands following [35], and Fc was subsequently estimated from, where NDVImax and NDVImin are NDVI for fully vegetated and bare surfaces, which are taken to be 0.85 and 0.15, respectively [36] ( Donohue et al., 2013).The coefficient n is a function of leaf orientation distribution within the canopy and was set a value of 2 following [37].It should be noted that all three parameters in Equation (10) were set somehow objectively, as different studies reported different parameter values (e.g., [38]).However, the sensitivity analysis suggests that the estimated surface SM is not sensitive to variations in these parameters (see Section 4.4).All remote sensing data were re-projected into the Universal Transverse Mercator (UTM) projection and resampled to a spatial resolution of 1000 m.

Other Data Sources
In addition to MODIS data, other two datasets are needed to run the algorithm.The first one is meteorological data, which was provided by the China Meteorology Data Center [39] at the Naqu station (Figure 3) to represent the regional mean climate condition (i.e., wind speed and vapor pressure).The air temperature for each pixel was corrected by considering the altitude effect based on DEM data (SRTM [40]) and a lapse rate of −6.5°C/km.The daily meteorological measurements were downscaled to the hourly timescale based on the algorithm proposed by [41].Then, the hourly meteorological data at satellite overpass time were used in subsequent calculation.The second dataset is soil texture, which was measured at each site.Field capacity (θF) was estimated to be the water content retained in the soil at −0.02MPa of suction pressure, which is midway of most reported θF values (−0.01 MPa to −0.33 Mpa) [42,43].The VG-M model was adopted to describe the soil water retention curve whose parameters were estimated from measured soil texture and organic content using the method given by [44].The residual water content (θR) is given directly by [44].Both meteorological and soil texture data at the site level were spatially interpolated using the Inverse Distance Weighting (IDW) method to obtain spatially-continuous inputs over the study domain.
SM output from four land surface models (LSM), including Noah 2.7 [45], Mosaic [46], Variable Infiltration Capacity model (VIC, [47]), and Community Land Model 2.0 (CLM, [48]) with a 1° spatial resolution and 3-h temporal resolution in Global Land Data Assimilation System-1 (GLDAS-1) [23] of the upper centimeters of the soil profile (10 cm for Noah, 2 cm for Mosaic, 10 cm for VIC, and 4.5 cm for CLM) was used for comparison with the developed algorithm and ground-based measurements.Surface SM output from these LSMs of two grid cells corresponding to the study site were spatially averaged for this comparison.

Results and Discussion
The current study only considered non-frozen periods (1 August2010~31 October 2010 and 1 June 2011~31 August 2011) for model testing because: (1) only liquid water content was measured within frozen soils during winter; and (2) soil freezing and thawing processes lead to a more complex water and heat transfer dynamics within soil and between soil and the atmosphere, which may not be captured by the trapezoid framework.In addition, we estimated the theoretical VI/LST trapezoidal space for each pixel to minimize uncertainties rising from non-homogeneous atmospheric forcing over the entire study region.

Validation with Site Observations
As shown in Figure 4, the site-level estimates and observations agree reasonably well with each other during the entire study period, with the coefficient of determination (R 2 ) of 0.73 and the root-mean-square error (RMSE) of 0.06 m 3 •m −3 , which falls within the AMSR-E mission accuracy requirement (i.e., the RMSE ≤ 0.06 m 3 •m −3 ).However, it is noted that the trapezoid method underestimated surface SM at the high end, resulting in an overall underestimation of 0.02 m 3 •m −3 (mean bias, defined as estimated mean minus observed mean).This may be related to the underestimation of field capacity of the surface soil layer or overestimation/underestimation of the warm/cold edge.According to the trapezoidal principle, actual soil water content is bounded by two critical values, which are the field capacity and the residual water content, respectively (Figure 1).The estimated field capacities from [44] were all lower than 0.6 m 3 •m −3 , whereas the observed water content reached 0.7 m 3 •m −3 .On the other hand, it is possible that the observed high SM reflected water content between field capacity and saturated water content, when there were rainy events shortly before the measurements or there was water ponding over the soil surface.Nevertheless, SM higher than the field capacity is not addressed in the developed method.
Another possible reason for the discrepancy is the mismatch between point-scale measurements and the 1 km-resolution MODIS input.Although the study region is relatively flat with uniform vegetation cover, SM may still have great local variability [49].Therefore, point observations may not be able to fully represent the average SM status within a MODIS pixel.Comparison of estimated seasonal evolution of daily station-averaged surface SM against observations is shown in Figure 5.As MODIS images were often contaminated by clouds in the region, only days with available MODIS data at more than 13 sites were used for this analysis.Zhao et al. [32] suggested that averaging SM measurements from across the 13 sites can effectively reflect SM averaged with all 50 sites with an accuracy of R ≥ 0.99 and RMSE ≤ 0.02 m 3 •m −3 for the study region.Generally, the simulated SM from the trapezoidal approach captures the observations fairly well, indicating the reliability of the algorithm over a wide range of SM conditions.However, consistent with the results in Figure 4, the station-averaged estimates are generally lower than the observed ones at high soil water content.

Spatial Distribution of Soil Moisture
An example of the spatial distribution of estimated surface SM for 2 September 2010 is shown in Figure 6b.The MODIS images acquired on this day were generally cloud-free over the entire domain.Up-scaled regional surface SM from observed points by applying the IDW method is shown in Figure 6a.Overall, both results are in good agreement, with surface SM being lower in the southwest and middle parts and higher in the northeast parts of the region.However, a slight difference was observed for the highest SM zone.The highest SM of the trapezoidal method occurs in the northwest side of that estimated with the IDW method.This discrepancy may be attributed to a limited number of observations in the northwest area (see Figure 3), which may have significantly reduced the reliability of the IDW method in interpolating soil moisture and soil texture.On the other hand, the high SM zones in Figure 6bcorrespond to high NDVI (Figure 6c) and low LST (Figure 6d), which is a reflection of high moisture content, suggesting that the estimated SM from the trapezoidal method may be more reasonable.In addition, the estimated surface SM in the southeast part of the study site is slightly lower than the spatially interpolated one.For other days, the standard deviation of estimated surface SM agrees reasonably well with measurements (Figure 7), suggesting a good performance of the trapezoidal method in capturing the spatial variability in surface SM.

Theoretic Boundary vs. Observed Boundary
Accurate determination of the boundaries (warm and cold edges) of the Fc/LST space is the most critical procedure of the trapezoidal method, as they jointly determine the upper and lower limits of the surface SM conditions within the domain and the relative position of each pixel in the trapezoid.However, traditional ways to determine the boundaries are based on visual examination of the images being used, which suffers from great subjectivity and uncertainties [12,13].Simple regression of these boundaries based on extreme surface temperature values at each Fc interval is highly sensitive to outliers (e.g., standing water bodies, clouds, and terrain effects) and the Fc interval specified by the operator.The regressed limiting edges also depend on satellite images with varying sizes and spatial resolutions.Figure 8 shows Fc/LST scatter points on 2 September 2010 over the entire study region.For the lower boundary, the observed cold edge (the lower envelope of scatter points in Figure 8) agrees well with the theoretical one (spatially averaged Ta); whereas for the upper boundary, the observed warm edges are far below the theoretical warm edge.In addition, three different warm edges were determined by three different experienced operators based on the observed Fc/LST scatterplot.Four cases of the combinations of limiting edges were set up in the trapezoidal model to estimate surface SM.Results show that the case applied with the theoretical limiting edges gives the best surface SM estimates when compared with observations, while the other three cases applied with the observed warm edge all underestimate surface SM seriously (Figure 9 and Table 1).For the study day, the observed minimum surface SM is 0.18 m 3 •m −3 , which is significantly higher than the residual water content (~0.05 m 3 •m −3 ), suggesting a relatively wet surface condition.Under such conditions, extremely dry surfaces do not necessarily exist in the modelling domain and the underestimation of the warm edge is likely to occur because pixels suffering from antecedent soil wetness may be considered as the dry surface pixels by the visual examination method.In addition, even though based all on images, different limiting edges determined by different operators may bring further uncertainties in the trapezoidal method, which makes model output not deterministic or ambiguous in most cases.The study mentioned above suggests that determination of limiting edges based on visual examination/or simple regression, depending largely on subjectivity and spatial extent being studied, is far from satisfactory.The theoretical limiting edges adopted in this study effectively obviate those uncertainties, and therefore produce more reliable spatial patterns of surface SM estimates.Taking 2 September, 2010 as an example, statistics shown in Table 1 suggest that the theoretical limiting edges can greatly improve the accuracy of surface SM estimates compared with observed limiting edges in terms of an RMSE on an order of ~0.13 m 3 •m −3 and a mean bias on an order of ~0.1 m 3 •m −3 .

Figure 9.
Cumulative distribution curves of observed surface soil moisture and those estimated from the trapezoidal method with different boundaries edges.Cases 1, 2 and 3 correspond to different warm edges shown in Figure 7.
Table 1.Statistical analysis of surface soil moisture from observations and estimated from the trapezoidal method with different boundaries edges on 2 September 2010.Cases 1, 2 and 3 correspond to different warm edges shown in Figure 7.

Sensitivity Analysis
A local sensitivity analysis was conducted to examine how uncertainties in trapezoidal method-estimated surface SM could be apportioned to different sources of uncertainty in model input.The sensitivity to the ith forcing variable or parameter is assessed by calculating surface SM with a set of baseline parameters (SM0) and comparing this with surface SM calculated by varying the ith parameter (SM±); the sensitivity index is The variation ranges and steps of each input variable are set as 2 K for temperature variables, with a step of 0.5 K, and 20% for other variables, with a step of 5%.Data on 2 September 2010 was used for analysis because of better data quality.
As shown in Table 2, the surface SM estimate is most sensitive to changes in field capacity (θF).An increase of 20% in θF resulted in a 19.08% increase in estimated surface SM, and a decrease of 20% in θF resulted in a 19.02% decrease in surface SM estimates.In the current study, θF was underestimated in some regions and therefore resulted in underestimation of surface SM.However, another important hydrological variable, surface soil moisture availability (= / or EF in Figure 1) [17] does not require information on soil hydraulic properties.Thus, even under conditions without a priori knowledge of soil parameters, surface soil moisture availability can still be reliably retrieved by the developed trapezoidal method.
Temperature plays a fundamental role in determining the shape of the trapezoidal space.The estimated surface SM showed positive correlations with Ta but is negatively correlated with LST (Table 2).An increase of 2 K in LST and Ta resulted in a 6.06% decrease and 5.43% increase in estimated SM, while a 2 K decrease in LST and Ta could lead to a 6.06% increase and 5.22% decrease in SM estimates, respectively.
For other variables, results indicate that wind speed (u), αc_max, αs_max, and ea are negatively correlated with SM estimates, while θR, NDVImax, NDVImin and n show positive relationships with SM estimates.However, changes in none of those variables resulted in significant variability in estimated surface SM.This suggests that although there are uncertainties in determining some of those variables, e.g., αc_max, αs_max, NDVImax, NDVImin and n, it will not greatly affect the accuracy of surface SM estimates.

Comparison with Land Surface Models
Comparison of surface SM from the developed trapezoidal approach and output from four LSMs in GLDAS-1 and the one estimated with the trapezoidal method with surface SM observations are shown in Figure 10.Again, this comparison was made only on days with sufficient remote sensing data (data available at more than 13 sites).Two grid cells of the LSMs were averaged to represent the mean surface SM status over the study region (~1° × 1°).To better match the depth at which surface SM were measured (0-5 cm), LSM-based surface SM from the CLM model was obtained by averaging the uppermost two layers (0-4.5 cm) and those from other three LSMs were acquired directly from the first soil layer, i.e., 0-2 cm for Mosaic, 0-10 cm for VIC, and 0-10 cm for Noah.
Figure 10.Comparison of daily area-averaged surface soil moisture estimated from the trapezoidal method and four land surface models against soil moisture observations.Only days with MODIS data available at more than 13 sites were shown.
Results show that all four LSMs systematically underestimated the surface SM, while the remote sensing-based trapezoidal method slightly underestimated surface SM when surface SM is high and overestimated it when surface SM is low (Figure 10).Among all methods, the trapezoidal model performed best during the study period in terms of RMSE and mean bias.With respect to temporal dynamics, Noah and CLM performed better than the trapezoidal method, as evidenced by a higher R 2 and lower unbiased RMSE (RMSEunbiased) (see [50] for detail about RMSEunbiased).However, the RMSE and mean bias of Noah are ~two times and four times larger than the trapezoidal method, and those of CLM are ~three times and five times larger than the trapezoidal method, respectively, suggesting a systematical underestimation of surface SM by the two models (i.e., Noah and CLM).For the four LSMs, Noah and VIC have smaller underestimation followed by CLM, and the Mosaic model shows the largest negative mean bias.Nevertheless, this may be partially ascribed to different model structures, e.g., varying surface soil layers in different LSMs, as the shallower soil layer generally dries out faster than deeper ones.
As discussed in previous studies [51][52][53][54] reliability of LSMs depends largely on accuracy of atmospheric forcing and model parameters.The more inputs for a model, the higher uncertainties may be brought.Compared with LSMs, the remote sensing-based trapezoidal method requires much fewer inputs, e.g., precipitation and soil hydraulic conductivity required in these prognostic LSMs are not involved in the diagnostic trapezoidal approach.Alternatively, LSMs provide continuous outputs while remote sensing-based approaches are often subjected to various types of noises, e.g., cloud contamination.Considering the advantages of both approaches, a data assimilation framework using remotely estimated surface SM is advocated to improve performance of LSMs and would consequently provide better surface SM output.

Conclusions
This study develops an algorithm to retrieve surface SM based on visible, near-infrared, and thermal infrared remotely sensed information and the trapezoidal Fc/LST space.Different from other trapezoidal (or triangular) methods, theoretical solutions of the limiting edges were derived in the developed algorithm to reduce subjectivity and uncertainties in configuring the Fc/LST space arising from visual interpretation of satellite images.This feature makes the trapezoidal method applicable to regions with various land cover conditions, even over highly heterogeneous surfaces.As the extreme conditions theoretically exist and the estimation of surface SM is pixel-independent, the VI/LST space and surface SM can be accurately estimated for each pixel once the spatial distribution of meteorological variables can be well defined [13].
Performance of the proposed algorithm was tested at the SMTMN site in central Tibet with MODIS data.Results showed that the trapezoidal method is able to reproduce spatial and temporal patterns of observed surface SM with an RMSE of 0.06 m 3 •m −3 at the field scale and 0.03 m 3 •m −3 at the regional scale, respectively.In addition, a case study on 2 September 2010 demonstrates the importance of theoretical boundaries, as the surface SM estimated with the traditionally observed boundaries is significantly underestimated, which consequently resulted in the underestimation of surface SM.Compared with LSMs, the remote sensing-based trapezoidal method gives generally better surface SM estimates, while LSMs show systematic underestimation, which may be attributed to uncertainties in model inputs and parameters over the Tibet Plateau [51].Sensitivity analysis suggests that the trapezoidal method is most sensitive to field capacity and temperature but less sensitive to other meteorological observations and parameters.
Compared with microwave remote sensing-based SM estimates which have a mesoscale spatial resolution (e.g., ~40 km), our results show that there are merits of using thermal/visible satellite images to retrieve surface SM at a much finer scale due to their relatively short wavelength.On the other hand, the relatively short wavelength does not allow the sensors to penetrate through cloud covers, which makes images prone to be contaminated by clouds (particularly for humid regions where cloud covers are often thick).This would lead to non-continuous SM estimates both spatially and temporally.One solution is to combine the thermal/visible images with microwave remote sensing to make the most of both observations (i.e., finer spatial resolution for thermal infrared/visible remote sensing and spatiotemporal consistency and continuity for microwave remote sensing).This may be achieved by either using thermal infrared/visible information to downscale microwave-based SM estimates [9] or using microwave-based SM estimates to constrain SM from thermal infrared/visible information at large scales.

Figure 1 .
Figure 1.A sketch of the trapezoidal vegetation coverage/land surface temperature (Fc/LST) space.The nomenclature is given in Section 2.1.

Figure 2 .
Figure 2. Flow chart of the trapezoidal model for estimating surface soil moisture (SM).Trapezoids represent input variables or parameters, and rectangles represent intermediate variables or parameters.

Figure 3 .
Figure 3. Location of the study area and a sketch map of the Soil Moisture/Temperature Monitoring Network (SMTMN) site distribution.The black triangles represent the first 30 sites built in 2010 and the white triangles denote the 20 sites setup in 2011.The black pentagram shows the location of Naqu City.

Figure 4 .
Figure 4. Comparison of estimated surface soil moisture against observed ones at each observation site.

Figure 5 .
Figure 5.Comparison of daily station-averaged surface soil moisture with observed ones.Only days with MODIS data available at more than 13 sites are shown.

Figure 6 .
Figure 6.Spatial distribution of surface soil moisture estimated from: (a) the Inverse Distance Weighting (IDW) method; (b) the trapezoidal method for 2September 2010.The spatial patterns of concurrent normalized difference vegetation index (NDVI) and LST are shown in (c) and (d), respectively.Blank pixels in (b)-(d) indicate either water surfaces or cloudy areas.

Figure 7 .
Figure 7.Comparison of standard deviation of estimated surface SM against that of observed surface SM.Each point represents the spatial variability in surface SM (standard deviation) for each day.Only days with MODIS data available at more than 13 sites are shown.

Figure 8 .
Figure 8. Scatter plot of the Fc/LST space and its warm and cold edges over the entire study area for 2 September 2010.The theoretical warm and cold edges shown in the figure represent the mean positions of the warm and cold edges over all pixels.Cases 1-3 represent three warm edges determined from three different operators based on visual inspection of the VI/LST data cloud.

Table 2 .
Relative sensitivity of estimated surface soil moisture from the trapezoidal method to input variables at the SMTMN site on 2 September 2010.Numberswithin each bracket indicate variations in temperature variables (i.e., −2 °C-2 °C), and numbers outside each bracket show variations in other variables (i.e., −20%~20%).