Estimating Subpixel Surface Heat Fluxes through Applying Temperature-Sharpening Methods to MODIS Data

: Using high-resolution satellite data to perform routine (i.e., daily to weekly) monitoring of surface evapotranspiration, evapotranspiration (ET) (or LE, i.e., latent heat ﬂux) has not been feasible because of the low frequency of satellite coverage over regions of interest (i.e., approximately every two weeks). Cloud cover further reduces the number of useable observations, and the utility of these data for routine ET or LE monitoring is limited. Moderate-resolution satellite imagery is available multiple times per day; however, the spatial resolution of these data is too coarse to enable the estimation of ET from individual agricultural ﬁelds or variations in ET or LE. The objective of this study is to combine high-resolution satellite data collected in the visible and near-infrared (VNIR) bands with data from the MODIS thermal-infrared (TIR) bands to estimate subpixel surface LE. Two temperature-sharpening methods, the disaggregation procedure for radiometric surface temperature (DisTrad) and the geographically-weighted regression (GWR)-based downscaling algorithm, were used to obtain accurate subpixel land surface temperature (LST) within the Zhangye oasis in China, where the surface is heterogeneous. The downscaled LSTs were validated using observations collected during the HiWATER-MUSOEXE (Multi-Scale Observation Experiment on Evapotranspiration) project. In addition, a remote sensing-based energy balance model was used to compare subpixel MODIS LST-based turbulent heat ﬂuxes estimates with those obtained using the two LST downscaling approaches. The footprint validation results showed that the direct use of the MODIS LST approach does not consider LST heterogeneity at all, leading to signiﬁcant errors (i.e., the root mean square error is 73.15 W · m − 2 ) in LE, whereas the errors in the LE estimates obtained using DisTrad and GWR were 45.84 W · m − 2 and 47.38 W · m − 2 , respectively. Furthermore, additional analysis showed that the ability of DisTrad and GWR to capture subpixel LST variations depends on the value of Shannon’s diversity index (SHDI) and the surface type within the ﬂux contribution source area.


Introduction
As a major component of the surface energy balance and the water budget, evapotranspiration (ET) has substantial effects on global climate change, water management and crop yields [1][2][3]. ET is the water transferred from the land surface to the atmosphere, and the latent heat accompanying ET is LE, where L is the latent heat of vaporization. Satellite-based remote sensing has been identified as a more suitable mean of mapping the spatial distribution of ET or LE [4]. Nevertheless, the land surface is highly heterogeneous in terms of its geometric and physical aspects, and the transfer or exchange processes between the land surface and the atmosphere are likely to be nonlinear. Models and methods developed assuming a homogeneous land surface may confront serious problems in the estimation of ET or LE using medium-or coarse-resolution data [5]. Moderate Resolution Imaging Spectrometer (MODIS) data have a daily return time (one overpass per day at varying view angles). At present, such data are ideal for the simulative monitoring of regional-scale changes in surface energy and water fluxes. Studies have shown that a pixel resolution of 500 m or less is necessary for accurate assessment of surface energy balance and ET or LE changes [1,6]. However, the spatial resolution of the MODIS thermal-infrared sensor (i.e.,~1000 m) is low, and the influence of the spatial scale (mixed-pixel) effect is substantial, especially where complex patterns of land use exist. Due to differences in the geometry of the underlying surface and land cover, it is difficult to use the mean or a single parameter to represent entire pixels. Moreover, farms and fields in many regions of the world are significantly smaller than individual pixels; in particular, in developing countries, fields may be less than 0.1 ha (30 m × 30 m) in size [6]. Hence, coarse-resolution data cannot be used to estimate ET or LE from individual fields. Developing a model of heat fluxes for use with mixed pixels and using the spatial information of the study area to get more reliable results may represent a solution to this problem. Middle-or high-resolution data (e.g., from the Land remote-sensing satellite Thematic Mapper (Landsat TM/ETM+), the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), and the Environment and Disaster Monitoring and Forecasting Satellite Constellation (Chinese HJ-1A/HJ-1B)) have high resolution but a low frequency of satellite coverage over regions of interest (e.g., approximately every 2 weeks for Landsat TM and ASTER; 4 days for the Chinese HJ-1A/HJ-1B). Moreover, cloud cover further reduces the number of useable observations of surface conditions resulting in high-resolution satellite imagery of a region typically being available once a month, which is not very useful for routine ET or LE monitoring [7].
To estimate ET or LE using surface energy balance models (e.g., the two-source model (TSM) [8] or the Surface Energy BALance (SEBAL) model [9]) and remote sensing data, the sensible heat flux (H) is the most important component, and it is also complex to calculate [10]. H is the transfer of turbulent heat between the surface and the atmosphere. It is driven by temperature differences and controlled by resistances that depend on local atmospheric conditions and land cover properties. In particular, the land surface temperature (LST) and the dynamic parameters (i.e., the zero plane displacement height and the roughness length) are more influential than other input parameters, because the LST provides an indication of the net effect of land-atmosphere interactions, which results in the process of evapotranspiration [2]. Sensitivity analysis has shown that uncertainty in LST has a significant impact on the estimation of H, and a 3 K change in LST yielded a change in H of up to 75% in TSM. A similar change in the difference between surface and air temperatures led to changes in H of up to 45% in SEBAL [11]. However, middle-or coarse-resolution satellite cannot provide accurate subpixel LST estimates. Furthermore, in order to drive a model with remote sensing data, the zero plane displacement height and the roughness length are often calculated using a vegetation index [12]. The vegetation index is directly related to the dynamic parameters. Hence, it is a promising scientific theme, both theoretically and practically, to model and monitor heat fluxes using high temporal resolution LST estimates derived from middle-or coarse-resolution data (e.g., MODIS) and high spatial resolution vegetation index values derived from middle-or high-resolution data (e.g., that obtained by the, Chinese HJ-1A/HJ-1B satellites). It is also significant in expanding the scope of application of satellite data and improving the utilization of the data.
Within the last two decades, an alternative technique, called LST downscaling, has been developed to derive high-resolution LST data using the high-resolution visible/NIR bands that are available from either the same or different satellite systems [7,13]. The most widely used LST downscaling approaches are the disaggregation procedure for radiometric surface temperature (the DisTrad algorithm) [1] and its descendant, the algorithm for sharpening thermal imagery (the TsHARP model) [14]. These two approaches are classified as statistical models, and they are based on the assumption that the relationship between the normalized difference vegetation index (NDVI) and LST is scale invariant. The performance of these algorithms has been evaluated over different land cover types, including agricultural areas [15], urban areas [16], and other areas, some of which display high heterogeneity [17,18]. Recently, many improvements have been made to DisTrad and TsHARP so that they include a greater number of independent variables (e.g., the normalized difference built-up index (NDBI), the surface albedo, and the surface reflectance) [19][20][21]. Moreover, many studies have questioned the hypothesis that the relationship between LST and NDVI is stationary in space [19,22,23]. In fact, the relationship between LST and other land surface parameters is geographically variable. Geographically weighted regression (GWR) is an extension of traditional standard regression techniques that permits local rather than global parameter estimation. It can be used to investigate nonstationary relationships between the dependent variable and the independent variables using geographically varying regression coefficients. Duan et al. [23] proposed a GWR-based LST downscaling algorithm that is based on the assumption that the relationship between LST and other land surface parameters and the digital elevation model (DEM) is geographically variable. They also find that the GWR-based algorithm outperforms the TsHARP algorithm, and the newly added parameters can effectively describe the variation in LST. The DisTrad model, which is based on NDVI segmentation, is generally considered to be more accurate [24], but the effectiveness of the GWR-based algorithm and the DisTrad model have not been compared. The purpose of this work is to investigate the reliability of downscaling algorithms for use with MODIS LST and the estimation of subpixel LE based on downscaled LST values for use in watershed studies. We downscale the MODIS LST products by using the GWR-based algorithm and the DisTrad model, and the performance of the two models is compared with ground-based observations of surface temperature. In addition, a remote sensing-based energy balance model is used to study the effect of the disaggregation of LSTs on the estimation of LE. Moreover, the relative ability of the two downscaling methods to describe the spatial heterogeneity of mixed pixels is assessed by comparing the turbulent heat fluxes estimated from the MODIS LST to the turbulent heat fluxes estimated from the disaggregated temperatures. The differences in the heat flux estimates obtained using the two downscaling methods are analyzed and discussed. This paper is organized as follows. Section 2 describes the study area and datasets. Section 3 presents a methodology for estimating LE and the temperature-sharpening methods, which are based on statistical relationships. The results and an error analysis are presented in Section 4. The uncertainties associated with this study and further analysis are discussed in Section 5. Some conclusions are drawn in Section 6.

Study Area
Our study was conducted within an agricultural oasis, the Zhangye oasis, which is located in the central reaches of the Heihe River Basin (HRB) and in the arid region of Gansu Province in northwestern China (100.10 • E-100.66 • E, 38.68 • N-39.15 • N). The Zhangye oasis has a semi-arid climate. The annual precipitation is less than 200 mm, but the potential evapotranspiration is approximately 1200-1800 mm yearly, and the annual temperature is approximately 7 • C [23,25]. Its terrain is flat, and its elevation is approximately 1480 m. A large portion of the Gobi Desert and the alpine vegetation of the Qilian Mountains are located near the study area ( Figure 1). Within the oasis, the land surface is mainly covered by maize, and the growing season lasts from May to September. The artificial oasis is highly heterogeneous, which impacts the thermal-dynamic and hydraulic features. Consequently, monitoring and managing the water resources within the oasis during the crop growing season are urgently needed [26,27]. The HRB has long served as a test bed for integrated watershed studies, as well as land surface and hydrological experiments [28]. One major objective of the Heihe Watershed Allied Telemetry Experimental Research (HiWATER), namely the Multi-Scale Observation Experiment on Evapotranspiration (HiWATER-MUSOEXE), is to capture the strong land surface heterogeneities and associated uncertainties within a watershed [25,27,29]. Two nested observation matrices that covered 30 km × 30 km and 5.5 km × 5.5 km were in place within the study area from May to September 2012. Within the larger observation matrix, four eddy covariance (EC) systems and one superstation (EC15) which is very high and equipped with a lot of observation system were installed in the oasis-desert ecosystem [25]. Within the oasis, a total of 17 automatic meteorological stations (AMSs) were distributed, including one in a vegetable field (EC01), one in a residential area (EC04), one in an orchard (EC17), and 14 in maize fields ( Figure 1).
Remote Sens. 2017, 9, 836 4 of 29 watershed studies, as well as land surface and hydrological experiments [28]. One major objective of the Heihe Watershed Allied Telemetry Experimental Research (HiWATER), namely the Multi-Scale Observation Experiment on Evapotranspiration (HiWATER-MUSOEXE), is to capture the strong land surface heterogeneities and associated uncertainties within a watershed [25,27,29]. Two nested observation matrices that covered 30 km × 30 km and 5.5 km × 5.5 km were in place within the study area from May to September 2012. Within the larger observation matrix, four eddy covariance (EC) systems and one superstation (EC15) which is very high and equipped with a lot of observation system were installed in the oasis-desert ecosystem [25]. Within the oasis, a total of 17 automatic meteorological stations (AMSs) were distributed, including one in a vegetable field (EC01), one in a residential area (EC04), one in an orchard (EC17), and 14 in maize fields ( Figure 1).

Datasets
Various types of remote sensing and ancillary data were used in this study, including highresolution satellite data (HJ-1B images), low-resolution satellite data (Terra/MODIS images), Digital Earth data (digital elevation models or DEMs), land use classification data [3], and in situ data describing radiation, heat fluxes and meteorological measurements.

Remote Sensing Data
Three kinds of remote sensing data were used in this study, specifically high spatial resolution data from the Chinese HJ-1B satellite, low spatial resolution MODIS data and the ASTER Global Digital Elevation Map (GDEM), which was used as auxiliary data. The HJ-1 A/B satellites, which were launched in September 2008, carry a total of four CCD cameras (i.e., CCD1 and CCD2 for HJ-1

Datasets
Various types of remote sensing and ancillary data were used in this study, including high-resolution satellite data (HJ-1B images), low-resolution satellite data (Terra/MODIS images), Digital Earth data (digital elevation models or DEMs), land use classification data [3], and in situ data describing radiation, heat fluxes and meteorological measurements.

Remote Sensing Data
Three kinds of remote sensing data were used in this study, specifically high spatial resolution data from the Chinese HJ-1B satellite, low spatial resolution MODIS data and the ASTER Global Digital Elevation Map (GDEM), which was used as auxiliary data. The HJ-1 A/B satellites, which were launched in September 2008, carry a total of four CCD cameras (i.e., CCD1 and CCD2 for HJ-1 A and similar for HJ-1 B). These satellites have quasi-sun-synchronous orbits, an altitude of 650 km, and a swath width of 700 km [30]. Apart from the first band of the HJ-1 CCDs (B1 = 475 nm), which is ≈10 nm wider than that of Landsat TM/ETM+ (B1 = 485 nm), the other three bands are identical (B2 = 560 nm, B3 = 660 nm and B4 = 830 nm). The spatial resolution of the four CCD sensors matches the first four bands of Landsat TM/ETM+ (see Table 1). The side rotation (±30 • ) feature of the HJ-1 A/B CCDs confers an advantage over Landsat, in that the revisit time is shortened from 16 days to 48 h or less [31,32]. The images produced by the HJ-1 A/B CCDs are stable, and the performances of each band are balanced [33]. After correction (i.e., geometric correction, radiometric calibration and atmospheric correction), the quality of the images is very similar to that of Landsat-5 TM [34]. The potential for high-quality ET or LE monitoring using the fine spatial resolution of the HJ-1 CCD sensors has already been demonstrated [24,35]. The data from the HJ-1 CCD sensors are freely available from the China Center for Resources Satellite Data and Application (CRESDA) website (http://218.247.138.121/DSSPlatform/index.html). We used the HJ-1B image data that cover the HRB region in 2012. Because many algorithms for estimating evapotranspiration using surface parameters require clear-sky conditions, we combined available information on data quality with visual interpretation to select satellite images without clouds. The selected images made up eleven groups, where a group consists of the matching IRS and CCD data covering the study area. The data making up these groups were collected on 19 June, 30 June, 8 July, 27 July, 3 August, 15 August, 22 August, 29 August, 2 September, 13 September, and 14 September.
The HJ-1B satellite data covering the HRB were pre-processed using methods including geometric correction, radiometric calibration, and atmospheric correction [24]. To estimate the regional LE, various surface variables, including the downward shortwave radiation [36], the downward longwave radiation [24], the surface albedo [37,38], the land surface temperature [24], the emissivity [39], the leaf area index (LAI) and the fractional vegetation coverage (FVC), were retrieved by applying the lumped method to the HJ-1B data (CCD/IRS) [24].
The MODIS sensor was launched on 18 December 1999 aboard the Terra satellite and has been in operation for nearly 18 years. It provides global data that are used to monitor long-term regional land surface heat and water vapor fluxes and changes in the Earth system [40][41][42]. The MODIS data have high temporal resolution (i.e., at least once every day) and high spectral resolution (i.e., 36 bands in total). The parameters used in the estimation of heat fluxes, such as NDVI, albedo, LSTs and emissivity, can be freely downloaded (https://ladsweb.nascom.nasa.gov/data/search.html). The Remote Sens. 2017, 9,836 6 of 30 MODIS products that cover the study area and contain LST-related parameters were downloaded and corrected for geometry and projection effects.
The ASTER GDEM is a product generated from optical data collected by the ASTER instrument onboard NASA's Terra satellite. This dataset is the only DEM that covers~99% of the entire land surface at high resolution. In this paper, we mainly use the ASTER GDEM v2, which has an accuracy of 17 m at the 95% confidence level [43]. The three kinds of satellite data were registered and matched with each other in the same projection.

Data from the HiWATER Experiment
As mentioned above, HiWATER was designed as a comprehensive ecohydrological experiment, and it was implemented in the HRB. MUSOEXE, a component of HiWATER, was implemented to reveal the spatial heterogeneities of heat and water fluxes by constructing two flux matrices [44]. In this study, EC data from 19 towers within the two matrices were used. The details of these EC systems are shown in Table 2, and the distribution of the stations is shown in Figure 1. EC16 and HHZ were excluded because they lacked surface net radiation (R n ) and soil heat flux observational data. Note: EBR represents energy balance closure, and the calculation method can be found in [28].
The ground observation data include the H and latent heat flux (LE). The consistency of all of the EC systems and the quality of the data were ensured using the following steps. First, the consistency of the EC instruments was tested during an inter-comparison campaign in the flat, open Gobi desert, over a surface covered by coarse grain sand and small pebbles with withered sparse scrub vegetation. This campaign was completed before HiWATER-MUSOEXE was conducted [29]. Second, all of the raw data were acquired at 10 Hz and processed. The processing included spike removal and corrections for density (WPL correction) [45,46]. Third, a four-step procedure was performed to control the quality of the EC data; for more details, see [24,29]. The soil heat flux (G) was measured by using three soil heat plates at a depth of 6 cm at each site, and the surface G was calculated based the soil temperature and moisture above the plates [47]. Energy imbalances are common in ground flux observations. To force the energy balance, the conserving Bowen ratio (H/LE) and the residual closure technique are often used. Computing the LE as a residual variable may be a better method for producing energy balance closure under conditions with large LE [48]. Thus, the residual closure method was applied because the "oasis effect" was distinctly observed in the desert-oasis system on clear days during the summer [24]. The data obtained during the HiWATER experiment can be freely obtained from the Heihe Plan Data Management Center (http://www.heihedata.org/).

Model of Remote Sensing-Based LE Estimation
The surface energy balance describes the transfer of energy between the land surface and the atmosphere. The land surface energy balance equation is the basis of remote sensing-based LE estimation, and it is commonly expressed as follows: H is the flux caused by the temperature gradient between the land surface and the air, whereas LE is the change in energy due to evapotranspiration or condensation. Atmospheric turbulence at the boundary layer is the motivation for the heat transfer between the land surface and the atmosphere. LE can be derived as the residual of the energy balance equation after R n , G and H are known. The models that permit estimation of R n , G and H using remote sensing data are introduced below.
R n is the main energy source of land surface heat and mass transfer and exchange, such as evaporation. The estimation of R n is of great importance in the remote sensing-based estimation of land surface heat fluxes. It is the difference between incoming and outgoing radiation, and it is calculated as follows: where R n is the surface net radiation (W·m −2 ), R s ↓ is the downward shortwave radiation (W·m −2 ), R s ↑ is the shortwave radiation reflected by the surface (W·m −2 ), R L ↓ is the downward atmospheric longwave radiation (W·m −2 ), and R L ↑ is the upward longwave radiation emitted by the surface (W·m −2 ). α is the surface albedo, ε s is the emissivity of the land surface, R L ↓= ε a σT 4 a , ε a is the effective emissivity of the atmosphere, σ = 5.67 × 10 −8 W·m −2 ·K −4 is the Stefan-Boltzmann constant, T a is the air temperature, and T s is the surface radiation temperature.
G is the heat exchange between the surface and the soil deeper underground, and it is related to the exchanges of water and heat within the soil. Generally, heat is transferred from depth to the surface at night. Studies have found that soil heat flux is correlated with R n and can be derived using land surface parameters [24,49]. Because canopies exert a significant influence on G, the fractional canopy coverage f c is used to determine the ratio of G to R n as follows [24]: where Γ c = 0.05 and Γ s = 0.315 represent the ratio of G to R n for full canopy coverage and bare soil.
H is the turbulent heat transfer between the surface and the atmosphere that is derived using the temperature gradient between them. According to gradient diffusion theory: where ρ is the air density (kg·m −3 ), C p is the specific heat of air at a constant pressure (J·kg −1 ·K −1 ), T aero is the radiative surface temperature (K), T a is the air temperature at a reference height (K), and r ac is the aerodynamic resistance (s −1 ·m), which is determined using the zero plane displacement height and the roughness length. r ac was calculated based on the Monin-Obukhov similarity theory (MOST) using a stability correction function [12,50,51]. Because remote sensing-based methods cannot obtain T aero , the value of T aero is usually replaced by the radiative surface temperature (T s ). It should be noted that T s is not strictly equal to T aero in Equation (4). The difference between these terms for homogeneous and full-coverage vegetation is approximately 1-2 K [52], and it can reach 10 K in sparsely vegetative areas [53]. The method of Chen that corrects for this discrepancy by adding an "excess resistance" r bh was used in this study [54]. The modified H can be specified as follows:

Temperature-Sharpening Methods Based on Statistical Relationships
The estimation of H is the most important and difficult step in the remote sensing-based assessment of the surface energy balance, in which radiative surface temperature and dynamic parameters play important roles. In this study, high-resolution vegetation index values from the HJ-1B data were used to downscale the coarser MODIS LST values. Hence, we can correct the errors in the heat flux caused by mixing by obtaining an "effective surface temperature" for use in the sensible heat flux. Two different LST downscaling methods are introduced in the following sections.

The DisTrad Downscaling Method
The DisTrad method is a statistical temperature-sharpening method that assumes that the negative correlation between NDVI and LST is invariant [1]. This method can combine high spatial resolution HJ-1B data and high temporal resolution MODIS data to estimate ET or LE on regional scales and over long periods [55]. The flow chart describing this method is shown in Figure 2 and is summarized as follows: (1) The HJ-1B NDVI data with a resolution of 30 m were aggregated to the 990 m resolution of the NDVI data. We first selected a subset of pixels from the scene where NDVI is most uniform at a pixel resolution of 990 m. Using the original NDVI data (NDVI 30 ), the coefficient of variation (CV; the standard deviation divided by the mean) is calculated among the 33 × 33 pixels (1089 pixels in total) that make up each 990 m NDVI pixel (NDVI 990 ). (2) Next, the NDVI 990 is divided into several classes (0 ≤ NDVI 990 < 0.2, 0.2 ≤ NDVI 990 < 0.5 and 0.5 ≤ NDVI 990 ). A fraction (25%) of the pixels having the lowest CV were then selected from each class. Lower CV values correspond to more homogeneous land surface values. (3) The relationship between NDVI 990 and T 990 can be established using the selected pixels.
whereT 990 (NDV I 990 ) represents the LST value estimated using NDVI on the 990 m scale. (4) The surface temperatures calculated using Equation (6) differ from the MODIS LST values (T 990 ) due to the influence of soil moisture and cannot be ignored. Based on the correlation between NDVI and LST is invariant, the simulated high-resolution LST data (T 30 (NDV I 30 )) can be calculated according to Equation (7) as follows: where ∆T 990 = T 990 −T 990 is the spatial variability due to the effects of soil moisture.

The DisTrad Downscaling Method
The DisTrad method is a statistical temperature-sharpening method that assumes that the negative correlation between NDVI and LST is invariant [1]. This method can combine high spatial resolution HJ-1B data and high temporal resolution MODIS data to estimate ET or LE on regional scales and over long periods [55]. The flow chart describing this method is shown in Figure 2 and is summarized as follows:  (1). The HJ-1B NDVI data with a resolution of 30 m were aggregated to the 990 m resolution of the NDVI data. We first selected a subset of pixels from the scene where NDVI is most uniform at a pixel resolution of 990 m. Using the original NDVI data (NDVI30), the coefficient of variation (CV; the standard deviation divided by the mean) is calculated among the 33 × 33 pixels (1089 pixels in total) that make up each 990 m NDVI pixel (NDVI990). (2). Next, the NDVI990 is divided into several classes (0 ≤ NDVI990 < 0.2, 0.2 ≤ NDVI990 < 0.5 and 0.5 ≤ NDVI990). A fraction (25%) of the pixels having the lowest CV were then selected from each class.
Lower CV values correspond to more homogeneous land surface values. (3). The relationship between NDVI990 and T990 can be established using the selected pixels.

GWR-based LST Downscaling Algorithm
GWR is an extension of traditional standard regression techniques, such as ordinary least squares, that permits local (rather than global) parameter estimation [56]. It can be used to investigate nonstationary relationships between the dependent variable and the independent variables using geographically varying regression coefficients. The main objective of the GWR-based downscaling algorithm is to address the spatially heterogeneous relationships between LST and the other land surface parameters. According to the surface characteristics of the HRB and previous studies of this area [19,23], NDVI, albedo and the digital elevation model (DEM) data were selected as auxiliary variables for use in LST downscaling. DEM is recognized as an important factor in characterizing variations in LST [23]. The surface albedo was selected because it influences the surface radiation balance and is mostly negatively related to soil moisture over barren agricultural areas [19]. Because longitude and latitude were the input parameters in the GWR model, these two parameters were not explicitly selected as auxiliary variables in the LST downscaling process. The specific steps followed in the GWR-based LST downscaling algorithm are shown in Figure 3. They are summarized as follows: (1) Perform spatial aggregation of the NDVI, DEM and albedo values through pixel averaging.
NDVI 30 , DEM 30 and albedo 30 denote these variables at the pixel resolution of HJ-1B, whereas NDVI 990 , DEM 990 and albedo 990 represent the aggregated variables at the pixel resolution of MODIS.
(2) Use GWR to establish a nonstationary relationship between T 990 and NDVI 990 , albedo 990 , and DEM 990 , which can be expressed as:

Accuracy Assessment of the LST Downscaling Results
The performance of the DisTrad downscaling method and the GWR-based LST downscaling algorithm can be evaluated from qualitative and quantitative perspectives. The qualitative analysis

Accuracy Assessment of the LST Downscaling Results
The performance of the DisTrad downscaling method and the GWR-based LST downscaling algorithm can be evaluated from qualitative and quantitative perspectives. The qualitative analysis was performed by visual examination of the results. The quantitative analyses were performed using the root mean square error (RMSE), the mean bias error (MBE), and the coefficient of determination (R 2 ). Figure 4 displays the spatial pattern of the LST data downscaled using DisTrad and GWR on 19 June 2012. For comparison, the spatial distribution of the MODIS LST is also shown. In comparison with the MODIS LST, both the DisTrad-and GWR-downscaled LST data effectively enhance the spatial resolution while preserving the characteristics and spatial distribution of the MODIS LST. Moreover, the results of the two models provide a more detailed description of the spatial variability of the thermal characteristics of the surface. However, a smoothing effect is observed in the LST downscaled using GWR when compared to the DisTrad-estimated LST because the regression process is based on the minimum mean square error (MMSE). Methods based on the MMSE tend to underestimate high values and overestimate low values; thus, the predicted parameters are frequently smoothed [23,57].  The ground-based land surface temperature (TS) was calculated using the Stefan-Boltzmann law from the AMS measurements of the longwave radiation fluxes [58]: where ↑ and ↓ are the in situ surface upwelling and atmospheric downwelling longwave radiation, respectively, and is the surface broadband emissivity, which is regarded as the value of the HJ-1B pixel corresponding to the AMS. Figure 5 shows the comparison of the LSTs obtained using the two downscaling methods with the ground-based measurements of land surface temperatures. As seen in Figure 5l, the R 2 , MBE and RMSE of the DisTrad downscaled LSTs are 0.83, 0.62 K and 2.59 K, respectively, and these values reflect better results than those associated with the GWR downscaled LSTs (0.82, 0.75 K, 2.66 K). On the whole, the differences between the two methods are not remarkable. However, the method that displays better performance varies depending on the date examined. For example, the GWR downscaled LSTs are better than the DisTrad downscaled LSTs on 19 June, 30 June, 13 September and 14 September, whereas DisTrad is superior to GWR on 8 July, 27 July, 3 August, 15 August, 22 August, 29 August and 2 September. The reason for these differences between the two techniques may be related to changes in land surface parameters (i.e., NDVI and albedo); the study period covers the entire life cycle of plants, extending from when they begin to grow to when they grow rapidly and eventually wither. The MBE difference between the DisTrad downscaled LSTs and the GWR downscaled LSTs ranges from −0.49 K~0.2 K, and the RMSE difference is −0.47 K~0.42 K. The results of the two methods vary substantially on 8 July and 29 The ground-based land surface temperature (T S ) was calculated using the Stefan-Boltzmann law from the AMS measurements of the longwave radiation fluxes [58]: where L ↑ and L ↓ are the in situ surface upwelling and atmospheric downwelling longwave radiation, respectively, and ε S is the surface broadband emissivity, which is regarded as the value of the HJ-1B pixel corresponding to the AMS. Figure 5 shows the comparison of the LSTs obtained using the two downscaling methods with the ground-based measurements of land surface temperatures. As seen in Figure 5l, of the LSTs. Therefore, the accuracy of the GWR downscaled LSTs is better at this time. As seen in Figure 6, the percentage contribution of albedo is higher on 3 August and 22 August (i.e., 26.03% and 26.89%, respectively) because irrigation occurred on these two days. Irrigation increases soil moisture and cools the land surface. Because albedo is mostly negatively related with soil moisture within barren agricultural areas [19], it displays a larger influence over the results.   To determine the reasons for the above phenomenon, we calculated the contributions of the different independent variables to the LSTs on each date separately by calculating the normalized regression coefficients [19]. The results are shown in Figure 6. The ability of the individual independent variables to explain the LSTs determines the accuracy of the two methods. NDVI explains only 49.39% and 47.09% of the LSTs on 19 June and 30 June, respectively, which is lower than the combination of albedo and DEM. Therefore, the method of joining albedo and DEM used in GWR is better than the DisTrad method. However, as the vegetation grows (see Figure 6 where the mean FVC starts to increase), NDVI plays an increasingly important role. NDVI explains over 53.67% of the LSTs on 8 July, 27 July, 3 August, 15 August, 22 August and 29 August, whereas the total contribution from albedo and DEM is only 15.20% to 46.33%. Thus, on these days, the downscaled LSTs produced by the DisTrad method are more accurate. Afterwards, while the mature crops (in our study area, maize is the main crop) gradually wither, the albedo and DEM explain over 66.61% of the LSTs. Therefore, the accuracy of the GWR downscaled LSTs is better at this time. As seen in Figure 6, the percentage contribution of albedo is higher on 3 August and 22 August (i.e., 26.03% and 26.89%, respectively) because irrigation occurred on these two days. Irrigation increases soil moisture and cools the land surface. Because albedo is mostly negatively related with soil moisture within barren agricultural areas [19], it displays a larger influence over the results. To further compare the downscaling results of these two methods for different objects, the downscaling results were verified separately according to the dominant surface type at each site. Here, we use the stations SSW and GB to validate the downscaling results for bare land; station EC04 was used to validate the downscaling results for buildings; station SD was used to validate the downscaling results for swampland; and the other stations were used to validate the downscaling results for vegetated areas. The assessments for the different land cover types are shown in Table 3. Overall, the results of the DisTrad downscaling method for areas covered with swampland, buildings and bare land are inferior to those of GWR, whereas the results in the vegetated areas are the opposite. The reason may be that NDVI is most suitable for vegetation monitoring within areas with moderate vegetation coverage [59], and it is considered unsuitable for studying areas with low vegetation coverage, such as urban build-up areas, where the detection capacity is only 30% [60]. The DisTrad downscaling results represent overestimates in areas covered with swampland (MBE > 0), To further compare the downscaling results of these two methods for different objects, the downscaling results were verified separately according to the dominant surface type at each site. Here, we use the stations SSW and GB to validate the downscaling results for bare land; station EC04 was used to validate the downscaling results for buildings; station SD was used to validate the downscaling results for swampland; and the other stations were used to validate the downscaling results for vegetated areas. The assessments for the different land cover types are shown in Table 3. Overall, the results of the DisTrad downscaling method for areas covered with swampland, buildings and bare land are inferior to those of GWR, whereas the results in the vegetated areas are the opposite. The reason may be that NDVI is most suitable for vegetation monitoring within areas with moderate vegetation coverage [59], and it is considered unsuitable for studying areas with low vegetation coverage, such as urban build-up areas, where the detection capacity is only 30% [60]. The DisTrad downscaling results represent overestimates in areas covered with swampland (MBE > 0), whereas they are underestimates (MBE < 0) in areas covered by buildings and bare land. On the other hand, GWR displays higher accuracy (|MBE| < 0.1 K) in both regions. Two possible reasons can explain this phenomenon. First, swampland contain water bodies. However, the DisTrad downscaling model eliminates water bodies when constructing the regression function, and the presence of the secondary term can cause considerable interference in the downscaling results for areas covered by water, buildings and bare land [14]. Furthermore, the GWR-based algorithm includes a spatially nonstationary relationship between LST and the other land surface parameters (NDVI, albedo and DEM) that is based on geographically varying regression coefficients [24].

Validation of Heat Fluxes Calculated Using MODIS LSTs
Before comparing the two downscaled LSTs in calculating surface heat fluxes, we must verify the LE model by comparing the remote sensing-estimated fluxes with the ground-based measurements. In this section, the MODIS LSTs, which have a resolution of 990 m, are resampled to a spacing of 30 m using the nearest neighbor method. The resampled LSTs and the surface VNIR variable, which have resolutions of 30 m, are applied to the LE algorithm. Figure 7 provides the turbulent heat flux results (which are termed the MODIS-LST heat fluxes) calculated using the MODIS LST data on 19 June 2012. The spatial distribution of the turbulent heat fluxes is related to the underlying land cover. The LE values are lower in areas covered by buildings and uncultivated land, including land patches in the desert areas, barren areas, and Gobi regions, and they are higher in the water-covered and agricultural areas within the oasis. Comparing Figure 7c,d, we see that the LE values (H values) of the barren areas in the southwest are higher (lower) than those of the desert areas in the southeast. Because the barren areas in the southwest border the Qilian Mountains, this result is due to snow melt and the downward movement of water. Comparing the four maps, we see that Figure 7c,d display obvious "boxy" anomalies, which are due to our direct use of the MODIS LSTs.
In recent years, many researchers have employed the footprint model in validating remote sensing-based LE estimates. When combined with the flux contribution source areas, this model solves the problems associated with mismatches between observed values and remote sensing-based estimates and leads to better results [61][62][63]. In this study, we employ a user-friendly tool that is based on the Eulerian analysis flux footprint model proposed by Kormann and Meixner (2001) to calculate the footprints of the function parameters [63,64]. Based on the relative weights of the pixels in which the source area is located, we determine the remote sensing-based estimates with the same spatial representation range as the observed values. The results of the footprint validation are shown in Figure 8 and Table 4. The R 2 , MBE, and RMSE values of LE are 0.75, −59.78 W·m −2 , and 73.15 W·m −2 , respectively, and those of H are 0.63, 38.90 W·m −2 , and 77.26 W·m −2 , respectively. The accuracy of LE is impacted by R n , surface G, and H, because it is calculated as a residual term. All of the input parameters may contribute to the error in LE; these errors are discussed in Section 4.4.
in the water-covered and agricultural areas within the oasis. Comparing Figure 7c,d, we see that the LE values (H values) of the barren areas in the southwest are higher (lower) than those of the desert areas in the southeast. Because the barren areas in the southwest border the Qilian Mountains, this result is due to snow melt and the downward movement of water. Comparing the four maps, we see that Figure 7c,d display obvious "boxy" anomalies, which are due to our direct use of the MODIS LSTs. In recent years, many researchers have employed the footprint model in validating remote sensing-based LE estimates. When combined with the flux contribution source areas, this model solves the problems associated with mismatches between observed values and remote sensing-based estimates and leads to better results [61][62][63]. In this study, we employ a user-friendly tool that is based on the Eulerian analysis flux footprint model proposed by Kormann and Meixner (2001) to calculate the footprints of the function parameters [63,64]. Based on the relative weights of the pixels in which the source area is located, we determine the remote sensing-based estimates with the same spatial representation range as the observed values. The results of the footprint validation are shown in Figure 8 and   As shown in Figure 8, the majority of the LE values are large because June, July, August, and September constitute the growing season, when soil moisture is adequate and ET is strong. ET leads to substantial reductions in air temperature; thus, the surface and air temperature gradient is small, leading to low values of H. The points with small LE values are influenced by uncultivated land. Within our study area, the barren areas, desert areas and Gobi regions make up the uncultivated land. The points in the scatterplot with large H values represent desert areas [24], where the H values reach approximately 300 W·m −2 . We also note that some points in the H scatterplot are less than zero, possibly due to the oasis effect or irrigation.
As shown in Figure 8, the majority of the LE values are large because June, July, August, and September constitute the growing season, when soil moisture is adequate and ET is strong. ET leads to substantial reductions in air temperature; thus, the surface and air temperature gradient is small, leading to low values of H. The points with small LE values are influenced by uncultivated land. Within our study area, the barren areas, desert areas and Gobi regions make up the uncultivated land. The points in the scatterplot with large H values represent desert areas [24], where the H values reach approximately 300 W m −2 . We also note that some points in the H scatterplot are less than zero, possibly due to the oasis effect or irrigation.

Compare of Heat Fluxes Calculated Using Different LST Datasets
To determine which of the LST downscaling methods is more suitable for calculating the subpixel surface heat fluxes, the downscaled LSTs obtained using the DisTrad and GWR methods were used to calculate the surface heat fluxes (here, we call these heat fluxes the DisTrad-LST and GWR-LST heat fluxes, respectively). These three methods were evaluated using (1) validation of the heat fluxes calculated with different LST downscaling methods using ground-based observational data and (2) quantitative analysis based on the spatial distribution and scatterplots of the energy balance components (i.e., H and LE).

Compare of Heat Fluxes Calculated Using Different LST Datasets
To determine which of the LST downscaling methods is more suitable for calculating the subpixel surface heat fluxes, the downscaled LSTs obtained using the DisTrad and GWR methods were used to calculate the surface heat fluxes (here, we call these heat fluxes the DisTrad-LST and GWR-LST heat fluxes, respectively). These three methods were evaluated using (1) validation of the heat fluxes calculated with different LST downscaling methods using ground-based observational data and (2) quantitative analysis based on the spatial distribution and scatterplots of the energy balance components (i.e., H and LE). Table 5 provides the in situ validation results of H and LE calculated using the DisTrad downscaled LSTs and the GWR downscaled LSTs. Compared to the validation results of the MODIS-LST heat fluxes in Table 4 and Figure 8, the DisTrad-LST heat fluxes display better retrieval accuracy than the GWR-LST heat flux, and the GWR-LST heat fluxes are generally better than the MODIS-LST heat fluxes. Also, the GWR-LST heat fluxes and the DisTrad-LST heat fluxes are both better than the MODIS-LST heat fluxes on all of the days examined; moreover, the MBE and RMSE values decreased and the R 2 values increased on most days. Using the data collected on 19 June as an example, the spatial distributions of H and LE calculated using the GWR-LST and DisTrad-LST are shown in Figure 9. Based on Figure 9, Figure 7c,d, the results obtained using the GWR-LSTs and DisTrad-LSTs are synoptically smoother than the MODIS-LST results; no "boxy" anomalies like those shown in Figure 7c,d are present. Because the GWR-LST heat fluxes and DisTrad-LST heat fluxes consider the heterogeneity of surface temperatures that could be a significant source of error in the LE models, the differences among the three heat fluxes result mainly from the differences between the sharpened LSTs and resampled LSTs of the subpixels at a resolution of 30 m. For example, the boundary between the oasis and uncultivated land becomes a belt of intermediate H and LE values because these mixed pixels include uncultivated land and vegetation. However, the temperatures of the mixed pixels in the coarse-resolution image (which shows the MODIS LSTs) are replaced by a single value. This result neglects the significant differences in the subpixel temperatures.   We also note that the difference in MBE between the heat fluxes determined from the downscaled LSTs and the MODIS-LSTs is larger at EC04 and EC11 than EC15. Figure 10 shows that the classes and temperatures of 10 × 10 subpixels at 30 m correspond to the pixels within the MODIS pixels, which have a resolution of 990 m, at the EC tower. After temperature sharpening, the distribution of temperatures at a resolution of 30 m agrees with the classification. EC04 shows the temperature sharpening resulting from different methods on 19 June. The temperature retrieved at the MODIS resolution was 305.24 K. Compared with the in situ measurement of 314.04 K, the temperature was underestimated at a resolution of 990 m. The mean value of LE calculated using the Inhomogeneity is an inherent property of the Earth's surface parameters. Because of its complexity, this heterogeneity can seriously affect the accurate acquisition of land surface heat fluxes, and it is one of the main factors that cause various surface parameters to change with scale [65,66]. Landscape inhomogeneity can be classified using two conditions: nonlinear vegetation density variations between subpixels (e.g., different types of vegetation mixed with each other) and coarse pixels containing different landscapes (e.g., vegetation or bare soil mixed with areas covered by buildings or water). To evaluate the effects of LST, stations EC04, EC15, and EC11, which display a typical severely heterogeneous surface, a uniform surface, and a weakly heterogeneous surface, respectively, were selected to analyze the heat fluxes calculated using the results of different LST downscaling methods. EC04 is used as an example because its land cover displays a complex pattern. Table 6 compares the turbulent heat fluxes calculated using the MODIS-LST, GWR-LST and DisTrad-LST methods. For EC04, the LE calculated using the GWR-LST is more consistent with the in situ measurements than the LE calculated using the MODIS-LST; the MBE and RMSE are much smaller, the R 2 is larger, and the MBE improves by 30.67 W·m −2 . However, the LE calculated using the DisTrad-LST is more accurate than the LE calculated using the GWR-LST. For EC15, the differences between the three methods were small; however, the LE calculated using the GWR-LST is more accurate than the LE calculated using the MODIS-LST and DisTrad-LST for EC11. The reason is related to the land cover types within the flux contribution source area, which we discuss in Section 4.3.2. We also note that the difference in MBE between the heat fluxes determined from the downscaled LSTs and the MODIS-LSTs is larger at EC04 and EC11 than EC15. Figure 10 shows that the classes and temperatures of 10 × 10 subpixels at 30 m correspond to the pixels within the MODIS pixels, which have a resolution of 990 m, at the EC tower. After temperature sharpening, the distribution of temperatures at a resolution of 30 m agrees with the classification. EC04 shows the temperature sharpening resulting from different methods on 19 June. The temperature retrieved at the MODIS resolution was 305.24 K. Compared with the in situ measurement of 314.04 K, the temperature was underestimated at a resolution of 990 m. The mean value of LE calculated using the MODIS-LST of the 10 × 10 subpixels was 342.22 W·m −2 , and the H was 158.96 W·m −2 . However, when substituting the downscaled LST into the LE model, the DisTrad-LST LE was 364.79 W·m −2 , and the GWR-LST LE was 372.08 W·m −2 , respectively, indicating that the errors caused by surface temperature nonuniformity cannot be ignored in the second landscape inhomogeneity condition.

Validation of the DisTrad-LST and GWR-LST Heat Fluxes
The land surface of EC15 is uniform and consists of pure pixels covering maize fields. Consequently, the temperature distribution at a resolution of 30 m was very homogeneous, and the ranges of the LST obtained using both downscaling methods displayed a range of 1.5 K. The temperature observed through in situ measurements was 301.99 K. For comparison, the temperature retrieved at the MODIS resolution was 302.94 K. The mean LE calculated using the MODIS-LST, DisTrad-LST and GWR-LST were 281.14 W·m −2 , 287.05 W·m −2 and 287.95 W·m −2 , respectively. For homogeneous areas, the differences among the three methods were not large (within 8 W·m −2 ). The weakly heterogeneous land surface at EC11 contained barley and maize, and it displays the first condition of heterogeneity (i.e., nonlinear vegetation density variations among subpixels). The mean LE values calculated using the MODIS-LST, DisTrad-LST and GWR-LST were 585.45 W·m −2 , 602.82 W·m −2 and 605.41 W·m −2 , respectively. The improvements in H and LE by the temperature sharpening were not obvious as the improvements at EC04, which contains a totally different landscape. At EC11, the range of temperature-sharpening results is 1.6 K. The mean temperatures retrieved using the MODIS-LST, DisTrad-LST and GWR-LST were 300.66 K, 299.40 K and 299.69 K, respectively, and the observed ground-based temperature was 297.83 K. Hence, the LE differences between the DisTrad-LST LE and GWR-LST LE are small. However, compared to MODIS-LST, temperature sharpening can still decrease the heterogeneity that results from thermal dynamics. The land surface of EC15 is uniform and consists of pure pixels covering maize fields. Consequently, the temperature distribution at a resolution of 30 m was very homogeneous, and the ranges of the LST obtained using both downscaling methods displayed a range of 1.5 K. The temperature observed through in situ measurements was 301.99 K. For comparison, the temperature retrieved at the MODIS resolution was 302.94 K. The mean LE calculated using the MODIS-LST, DisTrad-LST and GWR-LST were 281.14 W m −2 , 287.05 W m −2 and 287.95 W m −2 , respectively. For homogeneous areas, the differences among the three methods were not large (within 8 W m −2 ). The weakly heterogeneous land surface at EC11 contained barley and maize, and it displays the first condition of heterogeneity (i.e., nonlinear vegetation density variations among subpixels). The mean LE values calculated using the MODIS-LST, DisTrad-LST and GWR-LST were 585.45 W m −2 , 602.82 W m −2 and 605.41 W m −2 , respectively. The improvements in H and LE by the temperature sharpening were not obvious as the improvements at EC04, which contains a totally different landscape. At EC11, the range of temperature-sharpening results is 1.6 K. The mean temperatures retrieved using the MODIS-LST, DisTrad-LST and GWR-LST were 300.66 K, 299.40 K and 299.69 K, respectively, and the observed ground-based temperature was 297.83 K. Hence, the LE differences between the DisTrad-LST LE and GWR-LST LE are small. However, compared to MODIS-LST, temperature sharpening can still decrease the heterogeneity that results from thermal dynamics.

Comparison of the DisTrad-LST and GWR-LST Heat Fluxes
Based on Table 5, the MBE ranges of H and LE associated with the DisTrad-LST method and the GWR-LST method are −13.65 W·m −2~6 .99 W·m −2 and −7.23 W·m −2~1 6.32 W·m −2 , respectively. The GWR-LST heat fluxes are better than the DisTrad-LST heat fluxes on 19 June, 30 June, 13 September and 14 September, whereas the DisTrad-LST heat fluxes are superior to the GWR-LST heat fluxes on 8 July, 27 July, 3 August, 15 August, 22 August, 29 August and 2 September. The results obtained here are consistent with the results of the LST comparison because LST is the only difference between the two methods used to calculate the flux. To achieve an improved understanding of the spatial distribution of the differences between the three heat fluxes, we performed a scatterplot analysis by selecting 1000 × 1000 subpixels on 19 June as an example. The results are shown in Figure 11. Comparing the downscaled LST heat fluxes with those obtained from the MODIS-LST, it can be seen that the bias of LE estimation inside the oasis, Gobi, and desert areas is relatively small (pixels with a relatively large bias are covered with a red triangle on the left side of Figure 11, and the corresponding spatial distributions are marked in red on the right side of Figure 11). In those areas, the surface structure is very similar, resulting in a more uniform temperature distribution. However, where substantial discontinuities in land cover type, surface structure, or soil moisture occur (dry and wet areas, as well as areas with dense and sparse vegetation, appear together in many situations, such as at the boundaries between the oasis, the desert and the Gobi areas), the errors could be significant because of discontinuities in the LE model and its parameters. This phenomenon can be noted in Figure 11, where the pixels marked in red are mainly located in the boundary areas. Thus, in complex areas, we cannot directly use the MODIS-LST to calculate the regional LE, because the surface temperatures of these regions are either ignored or exaggerated in the low-resolution MODIS mixed pixels; thus, we cannot use a single value instead of the subpixel values. Comparing Figure 11b,d, we find that the red points obtained with the DisTrad-LST are much greater in number than those obtained with the GWR-LST (these points mainly indicate building areas within the oasis). The reason is that the regression relationship used in the DisTrad downscaling method is based on hierarchical pure pixels. Moreover, it assumes that the LST-NDVI relationship is stationary in space, but that is not the case [19,22]. The DisTrad downscaling method also tends to overestimate the LST of subpixels within cooler landscapes and to underestimate the LST of subpixels for warmer areas [1]. discontinuities in land cover type, surface structure, or soil moisture occur (dry and wet areas, as well as areas with dense and sparse vegetation, appear together in many situations, such as at the boundaries between the oasis, the desert and the Gobi areas), the errors could be significant because of discontinuities in the LE model and its parameters. This phenomenon can be noted in Figure 11, where the pixels marked in red are mainly located in the boundary areas. Thus, in complex areas, we cannot directly use the MODIS-LST to calculate the regional LE, because the surface temperatures of these regions are either ignored or exaggerated in the low-resolution MODIS mixed pixels; thus, we cannot use a single value instead of the subpixel values. Comparing Figure 11b,d, we find that the red points obtained with the DisTrad-LST are much greater in number than those obtained with the GWR-LST (these points mainly indicate building areas within the oasis). The reason is that the regression relationship used in the DisTrad downscaling method is based on hierarchical pure pixels. Moreover, it assumes that the LST-NDVI relationship is stationary in space, but that is not the case [19,22]. The DisTrad downscaling method also tends to overestimate the LST of subpixels within cooler landscapes and to underestimate the LST of subpixels for warmer areas [1]. The differences between the GWR-LST LE and the DisTrad-LST LE occurred mainly at the boundaries between oasis and water areas and uncultivated land (see Figure 11f). To analyze these differences further, we calculated the differences between the two methods at each site. Table 7 shows the information for sites that display differences between the two methods that are greater than 40 W m −2 . Table 8 shows the site information for small differences between the two methods. The difference between GWR-LST LE and DisTrad-LST LE is −63.51 W m −2~7 9.72 W m −2 . The biggest difference occurred in SD on 14 September, and the DisTrad-LST are better than the GWR-LST. However, in The differences between the GWR-LST LE and the DisTrad-LST LE occurred mainly at the boundaries between oasis and water areas and uncultivated land (see Figure 11f). To analyze these differences further, we calculated the differences between the two methods at each site. Table 7 shows the information for sites that display differences between the two methods that are greater than 40 W·m −2 . Table 8 shows the site information for small differences between the two methods. The difference between GWR-LST LE and DisTrad-LST LE is −63.51 W·m −2~7 9.72 W·m −2 . The biggest difference occurred in SD on 14 September, and the DisTrad-LST are better than the GWR-LST. However, in Section 4.1, we verify that the LST of DisTrad in areas covered by swampland are higher than the observed values. The reason is that we use the footprint model to validate the heat fluxes, and the verification result is determined by the contribution source area. Even at the same station, which of the two methods is better changed over time. For example, at EC11, the GWR-LST LE are better than the DisTrad-LST LE on 19 June and 30 June, but the opposite pattern occurs on 13 September. The reason may be related to the underlying surface type of the flux contribution source area. According to the calculation principle of the footprint model [63,64], the flux contribution source area is closely related to the wind speed and wind direction, and the wind speed and wind direction at the same station on different dates are likely to be different. To confirm this, we calculated Shannon's diversity index (SHDI) of the flux contribution source area (see Table 7). SHDI reflects the heterogeneity of a landscape, and it is sensitive to unbalanced distributions of various patches. In our paper, patches means different types of land use. SHDI also reflects the complexity and variability of each landscape patch type [67]. It was calculated using the following formula: where P i is the area ratio of the i th patch type in the entire landscape, and m is the number of patch types. SHDI = 0 indicates that the entire landscape consists of only one patch type. The main patch type in our study area include maize, spring wheat, barley, other crop, buildings and uncultivated land. The values of SHDI are shown in Tables 7 and 8. Compared to the measured values, the GWR-LST LE is always better than the DisTrad-LST LE when the SHDI value is high (>0.74). In such cases, the footprint area contains three or more different types of objects. Combined with specific surface types, we find that this situation belongs to the second condition of heterogeneity. At EC11 on 30 June, EC14 on 13 September, and EC09 on 8 July, the footprint area contains buildings (the proportions of buildings were 7.34%, 2.56%, 1.52%, respectively) and three types of vegetation. The SHDI values at EC03 on 3 August and EC01 on 30 June are lower than those of the site described earlier, because the footprint area is composed of only two types of vegetation and buildings (the proportions of buildings were 4.06% and 7.87%, respectively). At EC11 on 19 June, SD on 8 July and SD on 27 July, SHDI is high, but its underlying surface did not contain buildings, and the GWR-LST LE is still better than the DisTrad-LST LE because the footprint area at EC11 on 30 June contains bare land and SD on 8 July and 27 July contains water. When the underlying surface is very complex (i.e., it contains buildings or bare land and water), the probability that the DisTrad method will choose such a pixel to establish the regression equation is very small or even zero. Thus, the downscaled LSTs obtained using the DisTrad method for these pixels are unstable. However, the GWR-based algorithm includes a spatially nonstationary relationship between the LST and other land surface parameters (NDVI, albedo and DEM) based on geographically varying regression coefficients. In contrast, the GWR-based method takes into account a greater number of independent variables that can effectively express the distribution of surface thermal dynamics. The DisTrad-LST LE is better than GWR-LST LE when the SHDI value lies between 0.59 and 0.75. In this case, the footprint area contains a mixture of three or fewer different types of vegetation, and in vegetated areas, the LSTs obtained using the DisTrad downscaling method are superior to those produced by the GWR downscaling method (see Table 3). At EC11 on 13 September, the footprint area contains two types of vegetation which is maize and other crop. However, at EC15 on 13 September, SHDI = 0 and the footprint area contains only one type of vegetation. In this case, the DisTrad-LST LE is much better than the GWR-LST LE. Because the observation system at EC15 was a superstation with a 40-m-tall tower, the sensor may receive radiation from the tower, resulting in an overestimated albedo and eventually leading to the LSTs being overestimated. For example, at EC15 on 8 July and 13 September, the GWR-LSTs are significantly overestimated compared to the EC-LSTs and the DisTrad-LSTs. Based on Equations (2)-(4), R n and G are underestimated, whereas H is overestimated. Thus, in conjunction with Equation (1), the GWR-LST LE is significantly lower than the other two. However, in most cases, when the SHDI value is low (less than 0.2) and the underlying surface of the footprint area contains a single land cover type, the difference between the two methods is less than 10 W·m −2 . Table 8. The site information for small differences between the two methods.

Error Analysis
Many sources of error are unavoidable. Examples of these unavoidable sources of error include the calibration coefficients of the CCDs onboard the HJ-1B satellite and drifts in the IRS because of aging of the instruments; mismatches in the timescales of the EC data, which are averaged over time, and the remotely sensed turbulent heat fluxes, which are instantaneous; and the geometric corrections cause a half-pixel bias equal to or less than the deviation of the artificially subjective interpretation. LE is calculated as a residual item in the energy balance equation, whereas the errors of other three energy balance components accumulate in the LE term.

Errors from R n and G
Since G is calculated from R n , this paper mainly analyzes the errors associated with R n . According to a sensitivity analysis of Equation (2), R n is highly sensitive to R s ↓ and R L ↓ , whereas the albedo and ε s are less influential. The station validation results of R s ↓ and R L ↓ are shown in Table 8 On overall, the R 2 , MBE and RMSE values of R L ↓ were 0.73, 0.28 W·m −2 , and 21.24 W·m −2 , respectively. As seen in Table 9, the accuracies at SD and SSW were low. The low accuracies at the SD station potentially resulted from high humidity, which resulted in low at-nadir brightness temperatures and low retrieved R L ↓ . The models use surface temperature when estimating R L ↓ to approximate or substitute the near-surface temperature [24], and the inaccuracy of the SSW LST may influence the R L ↓ results.  (4), it can be seen that the four energy balance components are related to LST. Because LE is calculated as a residual term in the energy balance equations, the sensitivity of R n , G and H to LST was analyzed first. However, except for LST, land surface variables, including canopy height, FVC, and LAI, and meteorological variables, including wind speed, air temperature, and relative humidity, are also the important factors for a sensitivity analysis of H. Figure 12 presents a sensitivity analysis for R n , G and H. In this case, LST is 304 K, and it ranges from 297.5-311 K with a step size of 0.5 K; moreover, LAI ranges from 0.1-2.5 with a step size of 0.1. The canopy height is 0.69 m, and it ranges from 0.1-1.3 m. The wind speed ranges from 0.34-4.25 m·s −1 , the air temperature T a = 297.11 K, FVC = 0.54, and the relative humidity RH = 35.43%. In addition, the land cover is maize.
The air pressure is stable over short periods and has little effect on the LE results. Although H is sensitive to meteorological variables, such as wind speed (see Figure 12) and air temperature, the meteorological data used in the LE model are derived from ground observations; thus, the meteorological factors are relatively accurate. As shown in Figure 12, LAI, LST, and the canopy height are influential variables. H is sensitive to LAI when it is less than 1. The momentum roughness length increases as LAI increases, and the turbulent exchange is enhanced. However, when LAI is greater than 1, the plant canopy is regarded as a continuum that is not a sensitive variable. Because our study area is dominated by agriculture and the study period extended from June to September, the crops in the middle reaches of the HRB grew quickly; thus, LAI was usually greater than 1. Hence, LST and canopy height are the main sources of error.
(1) Errors in LST. As shown in Figure 12, R n and G are not very sensitive to LST, but H is more sensitive to LST. A 1 K bias in LST would result in an error in H of 18.86%, whereas the corresponding errors in R n and G are both only −0.96%. Therefore, in this paper, the LST is mainly affected by the impact of H on the role of LE. However, as we mentioned in Section 3.1, H depends on the temperature gradient between the surface and the atmosphere. Hence, the sensitivity of the LST is unstable and depends on the strength of the turbulence. The strength of the turbulence determines the mass and energy transport and the resistance to heat transfer, which influences the sensitivity to the LST. A weaker turbulence corresponds to a weaker LST sensitivity and vice versa.
(2) Errors in canopy height. The canopy height was obtained from a phenophase and classification map. Thus, the accuracy of the canopy height depends mainly on the plant growth state and the classification accuracy. Even within the same region, the canopy height of a crop can differ due to differences in seeding times and soil attributes, such as soil moisture and fertilization. For example, the orchard land use type was present at EC17; however, on our land use type classification map, the land use at EC17 was other crops, which includes orchards and vegetables. Thus, it was difficult to set the canopy height. In our study area, the canopy heights of orchards were approximately 4 m, most of the other crops were vegetables, which had a canopy height of 0.2 m. Hence, a value of 0.2 m would lead to an overestimate in LE. The LE estimates with incorrect canopy heights and the correct orchard canopy heights at EC17 are shown in Table 10. The LE-I relative error indicates the relative error between the LE obtained using the incorrect canopy height and the observed LE and is expressed as ((LE obtained using the incorrect canopy height)-(EC-LE))/(EC-LE) × 100%. Moreover, the LE-C relative error is the relative error between the LE from correct canopy height and the observed LE and is expressed as ((LE obtained using the correct canopy height)-(EC-LE))/(EC-LE) × 100%. Days with large LST biases were removed, and the bias between the model results and the ground-based observations decreased.

Discussion
Because HJ-1B has a single thermal channel, and single-channel LST retrievals may be unstable [1], we have not used the LST retrieved by HJ-1B to perform comparisons. As shown in Section 4.3, the direct use of MODIS-LST calculation results is ill-founded and causes obvious "boxy" anomalies in estimates of LE. In contrast, the temperature-sharpening algorithms used to estimate the GWR-LSTs and the DisTrad-LSTs are capable of decreasing the effects of the heterogeneity of the LST, which is consistent with results from previous studies [1,24]. However, the following problems require further analysis.
(1) In this paper, LST is derived from MODIS and NDVI uses 30 m information from the VNIR bands of HJ-1B. When LST and NDVI are derived from different sensors, inter-sensor calibration should be completed; the comparison of the thermal images across sensors is complex, due to the differences in the instantaneous field of view, the swath width, the orbital geometry, the spectral response functions and the overpass times [68].
(2) We use a one-source energy balance model when estimating LE. As discussed in Section 4.4.2, the accuracy of the estimated canopy height greatly influences the estimated turbulent heat flux. The canopy height is known a priori based on phenophase classifications and influences the accuracy of the surface roughness calculation. Multi-source remote-sensing data, such as active microwave and LiDAR data, could be used to obtain canopy heights [24], which would decrease the dependence on the accuracy of the classification. In addition, to correct the discrepancies between the remotely sensed radiative surface temperatures and the aerodynamic temperatures at the source of heat transport, a brief and well-performed parameterization scheme (under a uniform and flat plant surface) including excess resistance was used to calculate the aerodynamic resistance to heat transfer [24]. Since the surface of the HRB is very heterogeneous, multiple parameterization methods should be compared to select the optimum method.
(3) The energy balance closure has a significant influence on the evaluation of the model-estimated heat flux values. In this paper, the EC energy balance closure ratio was greater than 0.75. Studies have shown that extension of the averaging time [69], uncaptured low-frequency eddies [70], and the lack of an accurate accounting of heat storage terms [71], among other potential causes, may lead to these energy imbalances. The residual closure technique and the conserving Bowen ratio are often used to force the energy balance. We chose the former method because the conserving Bowen ratio method yields irrational sensible heat fluxes, due to the small or negative Bowen ratios that occur in oasis and desert environments [24]. The energy balance closure is problematic at times for turbulent flux systems and tends to be associated with significant discrepancies in LE [72].

Conclusions
In this paper, by using remote sensing data from multiple sources (i.e., MODIS and HJ-1B) and combining the advantages of different data types in terms of their spatial and temporal resolution, we enhance the accuracy of the estimated LE. The DisTrad downscaling method and the GWR-based downscaling algorithm have been shown to be effective in obtaining the surface temperatures within MODIS subpixels. However, the two methods have different advantages that depend on the descriptive ability of the independent variables at different times.
Compared to the turbulent heat fluxes calculated using different LST datasets, the direct use of MODIS LST approach does not consider LST heterogeneity at all, which causes significant errors in the latent heat flux (i.e., the mean bias error is −59.78 W·m −2 ), especially in the boundary areas and for mixed pixels. The errors obtained using DisTrad-LST LE and GWR-LST LE were −23.17 W·m −2 and −26.52 W·m −2 , respectively, and are less than that of MODIS LST based on the footprint validation results. However, this error is non-negligible. Additionally, it also shows that both LST downscaling methods can decreases the thermodynamic uncertainty caused by the land surface. As we discussed, although there are a lot of uncertainties in this paper, the differences between the turbulence heat fluxes calculated using two downscaling methods is only LST. Compared to DisTrad-LST LE and GWR-LST LE, the differences between the two methods occur in areas where the pattern of the surface types within the flux contribution source area is complex. The SHDI value can reflect the landscape heterogeneity, but in actual applications, the specific type of the underlying surface must be considered to distinguish between the advantages and disadvantages of the two methods. That is, when the SHDI value is high (>0.74) and the footprint area contains buildings or bare land and water, the GWR-LST LE outperforms the DisTrad-LST LE. On the other hand, when the SHDI value is between 0.59-0.75 and the footprint area contains three or fewer different types of vegetation, DisTrad-LST LE is better. When the SHDI value is low (<0.2) and the underlying surface of the footprint area contains a single type of vegetation, the difference between the two methods is less than 10 W·m −2 . As a sensitive variable of the ET model, canopy height is mainly determined by classification in this paper, and the application of classification at a 30 m resolution can improve the accuracy of the canopy height.
The HJ-1B satellite data have considerable advantages because of their high spatial resolution, and the MODIS data have fine temporal resolution. Because both satellites are still in operation and the data can be freely accessed, the long-term data are promising for applications involving the monitoring of energy budgets using a combination of data from the two satellites.