Empirical Estimation of Near-Surface Air Temperature in China from MODIS LST Data by Considering Physiographic Features

Spatially and temporally resolved observations of near-surface air temperatures (Ta, 1.5–2 m above ground) are essential for understanding hydrothermal circulation at the land–atmosphere interface. However, the uneven spatial distribution of meteorological stations may not effectively capture the true nature of the overall climate pattern. Several studies have attempted to retrieve spatially continuous Ta from remotely sensed and continuously monitored Land Surface Temperature (LST). However, the topographical control of the relationship between LST and Ta in regions with complex topographies and highly variable weather station densities is poorly understood. The aim of this study is to improve the accuracy of Ta estimations from the Moderate Resolution Imaging Spectroradiometer (MODIS) LST via parameterization of the physiographic variables according to the terrain relief. The performances of both Terra and Aqua MODIS LST in estimating Ta have been explored in China. The results indicated that the best agreement was found between Terra nighttime LST (LSTmodn) and the observed Ta in China. In flat terrain areas, the LSTmodn product is significantly linearly correlated with Ta (R2 > 0.80), while, in mountainous areas, the LSTmodn-Ta relationship differed significantly from simple linear correlation. By taking the physiographic features into account, including the seasonal vegetation cover (NDVI), the altitudinal gradient (RDLS), and the ambient absolute humidity (AH), the accuracy of the estimation was substantially improved. The study results indicated that the relevant environmental factors must be considered when interpreting the spatiotemporal variation of the surface energy flux over complex topography.


Introduction
The near-surface air temperature (Ta), which is traditionally recorded at 1.5-2 m above ground, is an indispensable parameter that drives the energy and water exchanges among the hydrosphere, atmosphere and biosphere.The energy exchange in terms of transpiration and/or evaporation from the canopy and soil via sensible and latent heat fluxes is mainly due to the vertical temperature and humidity gradients between the surface and the atmosphere [1][2][3].Given the ecological significance of Ta, the demand for accurate spatial data that indicate local and global environmental processes has risen greatly [4].Many studies have used spatial interpolations of Ta [5,6] and atmospheric reanalysis datasets [7,8] to monitor climate change and to derive ecological models.However, the gridded Ta data is generally limited by station coverage, especially in extensive mountainous regions, which has highly variable densities and elevation distributions of meteorological observations [9,10].Furthermore, the reanalysis climate variable is too spatially coarse to adequately represent the regional climate patterns.Remote sensing provides a high spatiotemporal precision of land surface temperature observations (LST), providing a feasible proxy for accurately retrieving Ta and thus improving knowledge of surface energy fluxes.
Given the strong correlation between the land surface temperature (LST) and Ta, recent efforts have developed an alternative method for retrieving Ta from high-resolution satellite-derived LST data [11][12][13].The LST products derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) have proven to be efficient for estimating Ta at multiple temporal scales [14,15] and across different land cover types [12,16].However, the estimation of Ta from remotely sensed LST is far from straightforward [11,17].Unlike the thermodynamic temperature (Ta), the LST is the surface radiometric temperature retrieved from thermal infrared radiance by an instantaneous view of the sensor [18,19].The difference between LST and Ta is strongly influenced by the surface characteristics and atmospheric conditions [20,21].The LST seems to be better at estimating minimum Ta than maximum Ta because the LST-Ta relationship is attenuated by solar radiation and susceptible to the more complicated heat transfer processes that occur during the day than those that occur at night [22][23][24].In contrast, Benali et al. [12] found that the daily maximum Ta could be accurately estimated from either daytime or nighttime Terra LST products with R 2 equal to 0.88 and 0.85, and root mean square error (RMSE) of 2.3 ˝C and 2.5 ˝C, respectively.Some studies have tried to improve the accuracy of the estimated Ta from MODIS LST via the inclusion of relevant environmental variables.Recondo et al. [25] retrieved daily Ta of peninsular Spain under multiple regression model by considering the Normalized Difference Vegetation Index (NDVI), water vapor, and spatiotemporal variables, increasing R 2 by 5%-19% and decreasing residual standard error by 20%-30%.Kloog et al. [17] applied a generalized additive mixed model combining NDVI, wind speed, elevation, and impervious surface with spatial smoothing technique to estimated Ta, yielding R 2 of 0.94.It is essential to sufficiently describe the spatial heterogeneity of the LST and Ta with the surface structures and climate conditions in different regions.
The LST-Ta relationship is mainly controlled by the surface energy balance, but it also depends on factors that are closely linked to energy processes [23,26].The vegetation cover representing the surface features plays a vital role in interpreting the ground thermal regime.For example, Sun and Kafatos.[27] used the NDVI-LST relationship for drought monitoring over North America, and Nieto et al. [28] used the temperature-vegetation method to obtain NDVImax, with a mean absolute error ranging between 2.8 ˝C and 4 ˝C.Atmospheric water vapor is closely coupled to the temperature through the evapotranspiration process [29], and it is vital for the accuracy of LST retrieval in the split-window algorithm [30,31].The high atmospheric water vapor content results in low LST due to less shortwave radiation reaching the surface, but its effects on Ta are positive by increasing downward radiation [26].Additionally, several studies have reported that both the LST and Ta experience complicated variability in mountainous areas with steep slopes, altitudinal gradients and heterogeneous landscapes [32][33][34].The complexity of mountain terrain tends to dominate the spatial distribution of land surface heat fluxes [35] due to the variations in solar radiation loading on different slope orientations [36][37][38].The LST appears to be well connected with the Ta on north-facing slopes, while on south-facing slopes, the LST is strongly modified by the land cover type [39].Oyler et al. [4] combined a topographic dissection index with MODIS LST to capture complex topoclimatic variations in the conterminous United States and thus yielded an accurate estimation of daily Ta with mean absolute error below 1 ˝C.To better account for the LST-Ta relationship, the physiographic factors should be considered in mountainous landscapes where weather stations are sparse.
Hence, the main objective of this study was to develop a simplified parameterization model for estimating monthly maximum, minimum, and mean air temperatures (Tmax, Tmin, and Tmean) from MODIS LST products.Because China is an expansive region with large elevation changes and climatic gradients extending from mountainous regions to flat plains and from tropical to alpine climates [40,41], the LST-Ta relationship has considerable spatial variability across China.Here, a stepwise multiple linear regression approach was explored to take into account the influences of the relevant physiographic variables on LST and Ta via the inclusion of the land surface heterogeneity (e.g., NDVI and elevation) and the atmospheric water vapor condition (absolute humidity, AH).The study is divided into three parts.First, we determined the reliable LST product via individual assessment of the linear relationship between the monthly Ta and the MODIS LST products from on board the Terra and Aqua satellites.Second, simple and multiple linear regression models were established to accurately derive Ta from LST according to the different terrain reliefs.Finally, the performances of the statistical models are validated using ground-based meteorological data.

Meteorological Data
The meteorological observation data of monthly Ta and relative humidity (RH) from 2000 to 2010 used in this study were obtained from standard meteorological ground stations at 2 m above ground.Additionally, the elevations around the weather stations were also considered.The observation data were continuously measured at 688 weather stations, located in various climate regions (Figure 1).These data were available from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn).The RH was then used to calculate the AH according to the formula provided by [42].The AH was selected in this work because it is a direct measure of the water vapor content in the surrounding air [43].The weather dataset from 2000 to 2005 was required for model construction, and the remaining data from 2006 to 2010 were used to validate the model.
where Max(H) is the maximum grid-cell elevation, Min(H) is the minimum grid-cell elevation, and P(A) is the flat plain area, A is the total area, i is the spatial window, and 500 is the benchmark mountain height in China.The RDLS value reflects the level of terrain relief relative to the surrounding environment.A value less than one indicates that the relative elevation difference does not exceed 500 m, while a value greater than one indicates a difference of more than 500 m.The elevation data were obtained from the weather stations.
To further characterize the influence of topography on the ground thermal regime, we developed a multi-scale topography index (RDLS) based on elevation.The RDLS for a specific window size reflects the altitudinal gradient of a grid cell relative to the surrounding terrain.We calculated the multi-scale RDLS across five spatial window sizes, including 6, 10, 25, 50, and 100 km.Generally, the RDLS value was concentrated in the range of 0 to 1, indicating that most of the altitude difference across China is less than 500 m.However, the value of the RDLS appeared to increase as the spatial window size increased, and its distribution became homogeneous.
As illustrated in Figures 1 and 2, the RDLS calculated using spatial window sizes less than 25 km that were capable of capturing distinct terrain patterns across China.The percentages of RDLS within the range of 0 to 1 were 71.5% and 60.3% for the 6 km and 10 km spatial resolutions, respectively, whereas the RDLS at spatial resolutions larger than 25 km were too coarse to fully describe the terrain relief between the plain and mountainous areas.The percentage of RDLS in different groups gradually decreased as the spatial resolution decreased.The percentage of the RDLS that was within the range of 0 to 1 decreased from 39.7% at the 25 km resolution to 13.1% at the 100 km resolution, and the percentage tended to distribute evenly across the entire range.Because the spatial window size has a serious impact on the RDLS value, the selection of window size should objectively account for the spatial terrain variation.The RDLS threshold to distinguish terrain relief also varied with the spatial window sizes.Therefore, the RDLS in the present study that was calculated at 6 km resolution is capable of capturing the distinctive terrain features and being in accordance with the resolution of the remotely sensed data.

Satellite Data
The satellite data were downloaded from the MODIS on board the NASA EOS Terra and Aqua satellites [44].The monthly NDVI data were from the MODIS vegetation index product (MOD13A3) at a spatial resolution of one kilometer.The land surface temperature and emissivity products (MOD11C3 and MYD11C3) retrieved from the MODIS TIR band via the generalized split-window algorithm [18,45] were used in this study.The Terra and Aqua LST products both contain the ascending and descending orbits with different overpass times, generating pairs of daytime and nighttime LST products in a 0.05 ˝ geographic climate modeling grid (5600 m at the Equator).The Terra descends (ascends) the equator at approximately 10:30 am (10:30 pm) local solar time, while Aqua descends (ascends) the equator at approximately 1:30 am (1:30 pm) [46,47].For simplicity, we renamed the daytime and nighttime LST products LSTmodd and LSTmodn for Terra and LSTmydd and LSTmydn for Aqua, respectively.Both the NDVI and LST data were reprocessed according to their quality control files; thus, only reliable pixels and good quality data were used.After comparing the four LST products with the observed Tmax, Tmean and Tmin, we ultimately chose MOD11C3 nighttime LST (LSTmodn) for the estimation of Ta.

Relief Degree of the Land Surface
Due to the influence of topographic relief, mountain temperature variability is more complicated because mountainous terrain experiences a broader range of spatiotemporal scales than flat terrain [34,48,49].The inhomogeneous surface heat flux modulates the turbulence structure, and there are strong lateral heat fluxes due to topography variability, violating the generally assumption of predominantly vertical heat transport [37,50].However, on relatively flat terrain and vegetation, the horizontal homogeneity leads to stable atmospheric conditions [51].In mountainous regions, the topography determines the temperature through the variations in solar radiation loading on different slope orientations [36][37][38].Furthermore, Ta is affected more by the wind and free-air advection from mountains than the local sensible heat fluxes [52][53][54].
To characterize the effects of the topography on the surface thermal processes, we derived a topographic index (the relief degree of the land surface, RDLS) from the elevation [55]: where Max(H) is the maximum grid-cell elevation, Min(H) is the minimum grid-cell elevation, and P(A) is the flat plain area, A is the total area, i is the spatial window, and 500 is the benchmark mountain height in China.The RDLS value reflects the level of terrain relief relative to the surrounding environment.A value less than one indicates that the relative elevation difference does not exceed 500 m, while a value greater than one indicates a difference of more than 500 m.The elevation data were obtained from the weather stations.
To further characterize the influence of topography on the ground thermal regime, we developed a multi-scale topography index (RDLS) based on elevation.The RDLS for a specific window size reflects the altitudinal gradient of a grid cell relative to the surrounding terrain.We calculated the multi-scale RDLS across five spatial window sizes, including 6, 10, 25, 50, and 100 km.Generally, the RDLS value was concentrated in the range of 0 to 1, indicating that most of the altitude difference across China is less than 500 m.However, the value of the RDLS appeared to increase as the spatial window size increased, and its distribution became homogeneous.
As illustrated in Figures 1 and 2, the RDLS calculated using spatial window sizes less than 25 km that were capable of capturing distinct terrain patterns across China.The percentages of RDLS within the range of 0 to 1 were 71.5% and 60.3% for the 6 km and 10 km spatial resolutions, respectively, whereas the RDLS at spatial resolutions larger than 25 km were too coarse to fully describe the terrain relief between the plain and mountainous areas.The percentage of RDLS in different groups gradually decreased as the spatial resolution decreased.The percentage of the RDLS that was within the range of 0 to 1 decreased from 39.7% at the 25 km resolution to 13.1% at the 100 km resolution, and the percentage tended to distribute evenly across the entire range.Because the spatial window size has a serious impact on the RDLS value, the selection of window size should objectively account for the spatial terrain variation.The RDLS threshold to distinguish terrain relief also varied with the spatial window sizes.Therefore, the RDLS in the present study that was calculated at 6 km resolution is capable of capturing the distinctive terrain features and being in accordance with the resolution of the remotely sensed data.

Relationships between the Ta and the Four LST Products over China
Both the Terra and Aqua LST products were compared with the ground-based monthly Ta to determine the optimum predictor for simulating Ta.As shown in the point density plots (Figure 3), the nighttime LST datasets (LSTmodn and LSTmydn) and the observed Ta are more linearly concentrated along the fitting line than the daytime datasets.Strong correlations were observed between the nighttime LST and Tmean and Tmin with minimal bias (R 2 > 0.90, RMSE and mean difference (MD) < 7 °C).Specifically, the LSTmodn tends to be more accurate for the estimation of Ta with lower intercepts, smaller biases and larger sample sizes than LSTmydn.For Tmax, the LSTmydd had good agreement with a slope of -1 and an intercept of less than 2 °C.This is most likely due to the fact that the Aqua overpass time (1:30 pm) is closest to the time when Tmax was recorded.However, the LSTmydd was inapplicable for retrieving Ta due to the data scarcity.As a consequence, LSTmodn seems to be best for linear estimation of Tmax, Tmean and Tmin among the LST products.

Relationships between the Ta and the Four LST Products over China
Both the Terra and Aqua LST products were compared with the ground-based monthly Ta to determine the optimum predictor for simulating Ta.As shown in the point density plots (Figure 3), the nighttime LST datasets (LSTmodn and LSTmydn) and the observed Ta are more linearly concentrated along the fitting line than the daytime datasets.Strong correlations were observed between the nighttime LST and Tmean and Tmin with minimal bias (R 2 > 0.90, RMSE and mean difference (MD) < 7 ˝C).Specifically, the LSTmodn tends to be more accurate for the estimation of Ta with lower intercepts, smaller biases and larger sample sizes than LSTmydn.For Tmax, the LSTmydd had good agreement with a slope of -1 and an intercept of less than 2 ˝C.This is most likely due to the fact that the Aqua overpass time (1:30 pm) is closest to the time when Tmax was recorded.However, the LSTmydd was inapplicable for retrieving Ta due to the data scarcity.As a consequence, LSTmodn seems to be best for linear estimation of Tmax, Tmean and Tmin among the LST products.

Relationships between the Ta and the Four LST Products over China
Both the Terra and Aqua LST products were compared with the ground-based monthly Ta to determine the optimum predictor for simulating Ta.As shown in the point density plots (Figure 3), the nighttime LST datasets (LSTmodn and LSTmydn) and the observed Ta are more linearly concentrated along the fitting line than the daytime datasets.Strong correlations were observed between the nighttime LST and Tmean and Tmin with minimal bias (R 2 > 0.90, RMSE and mean difference (MD) < 7 °C).Specifically, the LSTmodn tends to be more accurate for the estimation of Ta with lower intercepts, smaller biases and larger sample sizes than LSTmydn.For Tmax, the LSTmydd had good agreement with a slope of -1 and an intercept of less than 2 °C.This is most likely due to the fact that the Aqua overpass time (1:30 pm) is closest to the time when Tmax was recorded.However, the LSTmydd was inapplicable for retrieving Ta due to the data scarcity.As a consequence, LSTmodn seems to be best for linear estimation of Tmax, Tmean and Tmin among the LST products.

Model Building
As expected, good agreement between LSTmodn and the observed Ta (2000-2005) is observed in China (R 2 > 0.8, slope-1) (Figure 4).The RMSE ranges from 2.6 ˝C to 12.6 ˝C, and the MD varies from ´12.0 ˝C to ´1.1 ˝C.The LSTmodn tends to match Tmin best (R 2 , slope and MD-1) but shows large biases toward Tmax and Tmean.This is most likely due to the LSTmodn overpass time (10:30 pm), which is closer to the time when Tmin was recorded than when Tmax and Tmean were recorded, thus reducing the difference.However, some of the comparison scatter points are more scattered above the regression line.As indicated by the red circle in Figure 4, the linear LSTmodn-Ta relationship is significantly attenuated by this outlier dataset, revealing a triangular trend.The surface thermal dynamics are strongly influenced by the local surface features.Consequently, we identify this outlier dataset as collected from the Himalayas and the Hengduan mountain area located in southwest China where the elevation exceeds 1000 m (digital elevation model (DEM) > 1000 m) and the topographic relief is more than 600 m (RDLS > 1.2) (Figure 1).The complicated LSTmodn-Ta relationship in the Himalayas and the Hengduan mountain area is modified by the complex terrain and its related factors, such as elevation and solar radiation interception.

Model Validation
To assess the model performance, the Ta from routine meteorological stations (2006-2010) is used as the reference to validate the estimated Ta (Figure 5).The simple and multiple linear regression models exhibit good performance in the estimation of Ta in flat terrain and mountainous areas, respectively.The multiple linear regression model including NDVI, elevation, and AH can effectively interpret the LSTmodn-Ta relationship in the heterogeneous landscape of China.The Ta that are predicted over the entirety of China by combining the flat terrain area and mountainous area Based on the findings, we attempted to divide the comparison datasets into two parts, i.e., linear and nonlinear datasets, according to the altitudinal gradient.The former dataset exhibits a stable Remote Sens. 2016, 8, 629 7 of 15 and simple linear relationship (0.86 ď R 2 ď 0.95) between LSTmodn and Ta in the flat terrain areas (DEM ď 1000 m and RDLS ď 1.2); therefore, this dataset could be treated as empirical and used in applications.However, the latter dataset tends to show a triangular trend because it covers mountainous areas with relatively large local topographic reliefs (DEM > 1000 m and RDLS > 1.2).The nonlinear comparison dataset showed that the LSTmodn dramatically underestimated Tmax (MD = ´17.8˝C, slope = 0.68), resulting in the largest estimated bias (RMSE = 17.5 ˝C), followed by Tmean (RMSE = 10.4 ˝C, MD = ´9.7 ˝C) and Tmin (RMSE = 5.9 ˝C, MD = ´4.2˝C).Because the simple linear relationship is unable to fully clarify the LSTmodn-Ta relationship in the high altitude variation mountainous regions, the additional terrain-related variables should be incorporated to consider the effects of complex terrain on surface heat flux.A multiple linear regression model was built to select relevant environmental parameters into the LSTmodn-Ta relationship through stepwise regression analysis.The NDVI and LST are partially correlated with season [27,56,57], so the NDVI variable is incorporated only during the growing season from May to October.The simple and multiple linear regression models corresponding to different regions and seasons are shown in Table 1.

Model Validation
To assess the model performance, the Ta from routine meteorological stations (2006-2010) is used as the reference to validate the estimated Ta (Figure 5).The simple and multiple linear regression models exhibit good performance in the estimation of Ta in flat terrain and mountainous areas, respectively.The multiple linear regression model including NDVI, elevation, and AH can effectively interpret the LSTmodn-Ta relationship in the heterogeneous landscape of China.The Ta that are predicted over the entirety of China by combining the flat terrain area and mountainous area (multiple linear regression model) fit well with the observed Ta, yielding MD and intercepts close to 0 ˝C, as well as slopes and model efficiencies (EF) close to 1. Furthermore, the accuracy of the estimated Tmax and Tmean were significantly improved with the fitting line close to the 1:1 line.
In the flat terrain area, the monthly Ta is estimated as a function of LSTmodn and is highly consistent with the observed Ta.The frequency distribution of the estimation residual is symmetrical around zero (Figure 6).The relative frequency of the residual is largely concentrated within the range ´3 ˝C and 3 ˝C (66%-88%).The average observed Tmax, Tmean, and Tmin are 21.2 ˝C (ranging from ´8.1 ˝C to 41.0 ˝C), 15.4 ˝C (ranging from ´19.1 ˝C to 33.3 ˝C), and 10.8 ˝C (ranging from ´28.6 ˝C to 28.1 ˝C), respectively.The average simulated Tmax, Tmean, and Tmin are 21.5 ˝C (ranging from ´6.5 ˝C to 37.4 ˝C), 15.6 ˝C (ranging from ´14.3 ˝C to 32.5 ˝C), and 10.7 ˝C (ranging from ´21.0 ˝C to 28.7 ˝C), respectively, showing a perfect match with the observed Ta.However, in the southwest mountainous area, the multiple linear regression model is developed to derive Ta, and the outcomes are encouraging.As shown in Figure 5, the monthly Ta predicted by the multiple linear regression model appears to be more linearly correlated with the observed Ta than the simple linear regression model.By applying the multiple linear regression model, the correlation coefficient of R 2 increased substantially from 0.60-0.65 to 0.71-0.88.Meanwhile, the RMSE decreased from 3.49-3.87to 2.42-3.08,and the distribution of the estimation residual is more concentrated around 0 ˝C (Figure 6).Furthermore, the EF values of the multiple linear regression models are closer to 1, indicating a better description of the observed mean than that of the simple linear regression model.model appears to be more linearly correlated with the observed Ta than the simple linear regression model.By applying the multiple linear regression model, the correlation coefficient of R 2 increased substantially from 0.60-0.65 to 0.71-0.88.Meanwhile, the RMSE decreased from 3.49-3.87to 2.42-3.08,and the distribution of the estimation residual is more concentrated around 0 °C (Figure 6).Furthermore, the EF values of the multiple linear regression models are closer to 1, indicating a better description of the observed mean than that of the simple linear regression model.

The Accuracy of the Daytime and Nighttime LST for Ta Estimation
Terra nighttime LST shows the best agreement with the observed Ta in China (with the exception of Tmax), and the RMSE validation errors are 6.6 °C and 2.7 °C for Tmean and Tmin, respectively (Figure 3).The fact that the MODIS nighttime LST products result in better predicted Ta accuracies is also verified by the work of [14,53].During the nighttime, due to the radiative cooling at the surface, the air will generally be stably stratified and thus the turbulence generated by surface friction will be weakened [58,59].The earth surface behaves almost as an isothermal and homogeneous surface, simplifying the retrieval process and improving the estimation accuracy of the air temperature [24,60,61].However, in the daytime, the influences of direct solar illumination, the surface and weather conditions and the satellite itself result in more complex interactions between Ta and LST [17,24].The solar radiation is first absorbed by the ground, which then heats the near surface air.As the temperature increases and more thermal energy is concentrated at the surface, the daytime LST and Ta become increasingly decoupled [62], which is why the daytime LST is hotter than the observed Ta and their points are more scattered around the 1:1 line (Figure 3).In addition, a much larger effect from angular anisotropy on the remote sensing data products is expected for daytime LST than for nighttime LST due to the sunlight and shade effects within pixels during the daytime [15,63].
However, Mildrexler et al. [62] indicated that the daytime LST from Aqua had a strong positive correlation with Tmax globally, and Shen and Leptoukh [16] precisely derived Tmax from the Terra daytime LST over central and eastern Eurasia as the mean absolute error (MAE) varied from 2.4 °C

The Accuracy of the Daytime and Nighttime LST for Ta Estimation
Terra nighttime LST shows the best agreement with the observed Ta in China (with the exception of Tmax), and the RMSE validation errors are 6.6 ˝C and 2.7 ˝C for Tmean and Tmin, respectively (Figure 3).The fact that the MODIS nighttime LST products result in better predicted Ta accuracies is also verified by the work of [14,53].During the nighttime, due to the radiative cooling at the surface, the air will generally be stably stratified and thus the turbulence generated by surface friction will be weakened [58,59].The earth surface behaves almost as an isothermal and homogeneous surface, simplifying the retrieval process and improving the estimation accuracy of the air temperature [24,60,61].However, in the daytime, the influences of direct solar illumination, the surface and weather conditions and the satellite itself result in more complex interactions between Ta and LST [17,24].The solar radiation is first absorbed by the ground, which then heats the near surface air.As the temperature increases and more thermal energy is concentrated at the surface, the daytime LST and Ta become increasingly decoupled [62], which is why the daytime LST is hotter than the observed Ta and their points are more scattered around the 1:1 line (Figure 3).In addition, a much larger effect from angular anisotropy on the remote sensing data products is expected for daytime LST than for nighttime LST due to the sunlight and shade effects within pixels during the daytime [15,63].
However, Mildrexler et al. [62] indicated that the daytime LST from Aqua had a strong positive correlation with Tmax globally, and Shen and Leptoukh [16] precisely derived Tmax from the Terra daytime LST over central and eastern Eurasia as the mean absolute error (MAE) varied from 2.4 ˝C to 3.2 ˝C.The above studies also noted that their relationship greatly depended on the land cover type; a denser canopy resulted in more closely coupled daytime LST and Tmax.The Aqua/Terra LST overpass in the day is closer to the time when Tmax was recorded, which should be better at estimating Tmax than nighttime LST [14].However, Zeng et al. [64] found that a stronger correlation existed between Tmax and nighttime LST than between Tmax and daytime LST in estimating daily Ta over the US corn belt.Zhu et al. [65] and Lin et al. [66] found that the difference between daytime LST and Tmax was still too large with MAE being 6.21 ˝C and 8.12 ˝C for LSTmodd and LSTmydd, respectively.Furthermore, Huang et al. [67] retrieved daily Ta in central China, which showed similar results to this study that the RMSE of LSTmydd (2.6-4.2 ˝C) was larger than LSTmodn (1.9-2.5 ˝C).As further interpreted by Fu et al. [22] and Mostovoy et al. [68], the time discrepancy between the moment of the satellite overpass and the time when Tmax or Tmin occurred was not a key factor controlling the relationship between Ta and LST.

Using the Multiple Linear Regression Model to Derive Ta in Mountainous Areas
The LST and Ta are correlated due to the heat exchange between the surface and the atmosphere.This study found that combining the remote sensing LST and physiographic factors had the potential to improve the accuracy of estimating Ta in large regions.This finding was similar to the results of [12,69,70].Combining the simple and multiple linear regression models exhibited a good performance in the estimation of Ta over the entirety of China, yielding RMSE ranging from 2.1 ˝C to 3.3 ˝C and MD from ´0.13 ˝C to 0.26 ˝C (Figure 5).Zhang et al. [14] included the solar declination into the estimation of China daily Ta, resulting in a relatively large error where the MD were within 1.61 ˝C-3.32 ˝C and the RMSE from 2.08 ˝C to 4.10 ˝C.Yao et al. [71] retrieved monthly Tmean based on 72 weather stations in the southeastern Tibetan Plateau, yielding RMSE within 2.25 ˝C-3.23 ˝C, which showed the same accuracy as that in this study in mountainous areas (2.42 ˝C ď RMSE ď 3.08 ˝C).Zheng et al. [72] and Chen et al. [69] derived the monthly Ta over Northern China and the entirety of China, obtaining more accurate results with RMSE within 0.86 ˝C-1.13 ˝C and 0.80 ˝C-1.29 ˝C due to the combination of interpolation techniques.
The heat transfer process is strongly influenced by the local radiation budget and the atmospheric contribution [53,73].In mountainous regions, Ta is affected more by the wind and free-air advection from mountains than the local sensible heat fluxes; thus, it agrees less with LST [53,54].Additionally, in the high mountains, during snow-covered periods, the linear correlation between the air and ground surface temperatures becomes weaker than that in snow-free periods [74].Many studies have highlighted the importance of the vegetation cover, topography and climate conditions on the estimation of Ta from LST over complex terrain [4,17,75].To adequately interpret the relationship between LST and Ta, it is essential to take the local climate and surrounding environment into consideration.However, NDVI of alpine areas is unusable in the non-growing season because of snow and bare ground [24,76].The assumption of linear and negative slopes between LST and NDVI is, therefore, not applicable.Including NDVI in the regression model did not significantly improve the estimation accuracy, with the estimated R 2 improved trivially, from 0.73 to 0.76, 0.81 to 0.83, and 0.83 to 0.83 for Tmax, Tmean, and Tmin, respectively.By applying the multiple linear regression model in this study, the RMSE were 14%-37% smaller than those from the simple linear regression model (Figure 5).Correspondingly, the distribution of the estimation residual was more concentrated around 0 ˝C (Figure 6).Furthermore, the EF values and R 2 of the multiple linear regression models were closer to 1, indicating a better description of the observed mean than that of the simple linear regression model.
The Himalayas and Hengduan mountains located in the southwest of China are generally covered with alpine vegetation during warm growing seasons, while they are covered with dry grass residuals or snow during non-growing seasons.The seasonal vegetation greening tends to increase the spatial heterogeneity of the LST [77].Dry, bare, and low-density soil has relatively low thermal inertia, resulting in large temperature variations and high LSTs [78,79], while dense vegetation induces more transpiration and blocks incident shortwave radiation, efficiently preventing the surface temperature from increasing with the air temperature during the daytime, which cools the surface [27,36].Meanwhile, the climate conditions on the Tibetan Plateau also play major roles in the land surface-atmosphere interactions.Tanaka et al. [80,81] found that the ground surface and surface atmosphere layer on the Tibetan Plateau were very dry until the monsoon progressed.In summer, the strong activity of the convection cloud and mesoscale convection systems on the southern Plateau at a latitude of 33 ˝N, accompanied by the remarkable mountain uplift, resulted in heavy rainfall and high water vapor content in the mountain areas of southwest China [82,83].This area is rather close to a relatively large local topographic relief area (DEM > 1000 m and RDLS > 1.2).In addition, the Tibetan Plateau warming in summer induces an upward transfer of sensible heat flux, which then destabilizes the atmosphere and promotes convection, leading to more precipitation [84].This high atmospheric water content in the complex topography contributes to the large error in the LST determined from spectral emissivity [85,86].

Conclusions
In this study, we assessed the potential of MODIS LST for the estimation of monthly Ta through linear regression analysis of LST and observed the Ta from 688 meteorological stations in China.The Terra nighttime LST data (LSTmodn) was proven to be accurate enough for linear estimation of the monthly Tmean and Tmin with strong correlation (R 2 > 0.90) and minimal bias (RMSE and MD < 7 ˝C) among the other LST products despite their overpass time discrepancies.Because the LSTmodn and Ta are well coupled in the flat terrain areas (DEM ď 1000 m and RDLS ď 1.2), the Ta can be accurately estimated as a function of LSTmodn, yielding MD less than 0.5 ˝C, model efficiencies (EF) close to 1, and RMSE ranging from 2.08 ˝C to 3.31 ˝C.
However, in the mountainous areas of southwest China, a complex topography leads to a complicated thermal mechanism.To clarify the complex LSTmodn-Ta relationship, the physiographic factors closely linked to energy processes were taken into account.This study develops a feasibility proxy to retrieve Ta that includes the influences of the seasonal vegetation cover (NDVI) and the ambient absolute humidity (AH) by stepwise regression analysis according to the different topographical reliefs.The multiple linear model can effectively improve the estimated accuracy by significantly increasing the R 2 from 0.60-0.65 to 0.71-0.88 and decreasing the RMSE from 3.49 ˝C-3.87 ˝C to 2.42 ˝C-3.08 ˝C in the mountainous areas (DEM > 1000 m and RDLS > 1.2).This result suggested that a large terrain relief can significantly influence the spatial pattern of the climate and thus modify the surface energy flux.This finding becomes critically helpful in interpreting the LSTmodn-Ta relationship by considering the complicated hydrothermal circulation in complex terrain.Therefore, further studies should consider the local physiographic features contributing to the land surface-atmosphere interaction mechanism by reducing the bias of Ta determined from satellite LST.

Figure 1 .
Figure 1.Spatial distribution of the meteorological stations and the relief degree of the land surface (RDLS) in China.

Figure 1 .
Figure 1.Spatial distribution of the meteorological stations and the relief degree of the land surface (RDLS) in China.

Figure 2 .
Figure 2. The frequency distribution of multi-scale RDLS over China.

Figure 3 .
Figure 3.The point density plots of the land surface temperature (LST) products (LSTmodd, LSTmodn, LSTmydd, and LSTmydn) and the observed air temperature (Ta) datasets (Tmax, Tmean, and Tmin).The point density from low to high is expressed by the color ramp from blue to red.RMSE

Figure 2 .
Figure 2. The frequency distribution of multi-scale RDLS over China.

Figure 1 .
Figure 1.Spatial distribution of the meteorological stations and the relief degree of the land surface (RDLS) in China.

Figure 2 .
Figure 2. The frequency distribution of multi-scale RDLS over China.

Figure 3 .Figure 3 .
Figure 3.The point density plots of the land surface temperature (LST) products (LSTmodd, LSTmodn, LSTmydd, and LSTmydn) and the observed air temperature (Ta) datasets (Tmax, Tmean, and Tmin).The point density from low to high is expressed by the color ramp from blue to red.RMSEFigure 3. The point density plots of the land surface temperature (LST) products (LSTmodd, LSTmodn, LSTmydd, and LSTmydn) and the observed air temperature (Ta) datasets (Tmax, Tmean, and Tmin).The point density from low to high is expressed by the color ramp from blue to red.RMSE is the Root-Mean-Square Error, and MD is the Mean of the Differences between LST and Ta; the same abbreviations are used throughout the work.

15 Figure 4 .
Figure 4.The linear relationship between LSTmodn and monthly observed Ta in different regions.

Figure 4 .
Figure 4.The linear relationship between LSTmodn and monthly observed Ta in different regions.

Figure 5 .
Figure 5.The linear relationship between the monthly observed and the estimated Ta in different regions.The Ta predicted over the entirety of China is combined by the flat terrain area and mountainous area (multiple linear regression model).EF is the model efficiency.

Figure 5 .
Figure 5.The linear relationship between the monthly observed and the estimated Ta in different regions.The Ta predicted over the entirety of China is combined by the flat terrain area and mountainous area (multiple linear regression model).EF is the model efficiency.

Figure 6 .
Figure 6.The frequency distribution of the estimation residual in different regions.

Figure 6 .
Figure 6.The frequency distribution of the estimation residual in different regions.

Table 1 .
Regression models for deriving air temperature (Ta) from Terra nighttime land surface temperature (LSTmodn) with the relevant environmental factors.