Impacts of Land Cover and Seasonal Variation on Maximum Air Temperature Estimation Using MODIS Imagery

Daily maximum surface air temperature (Tamax) is a crucial factor for understanding complex land surface processes under rapid climate change. Remote detection of Tamax has widely relied on the empirical relationship between air temperature and land surface temperature (LST), a product derived from remote sensing. However, little is known about how such a relationship is affected by the high heterogeneity in landscapes and dynamics in seasonality. This study aims to advance our understanding of the roles of land cover and seasonal variation in the estimation of Tamax using the MODIS (Moderate Resolution Imaging Spectroradiometer) LST product. We developed statistical models to link Tamax and LST in the middle and lower reaches of the Yangtze River in China for five major land-cover types (i.e., forest, shrub, water, impervious surface, cropland, and grassland) and two seasons (i.e., growing season and non-growing season). Results show that the performance of modeling the Tamax-LST relationship was highly dependent on land cover and seasonal variation. Estimating Tamax over grasslands and water bodies achieved superior performance; while uncertainties were high over forested lands that contained extensive heterogeneity in species types, plant structure, and topography. We further found that all the land-cover specific models developed for the plant non-growing season outperformed the corresponding models developed for the growing season. Discrepancies in model performance mainly occurred in the vegetated areas (forest, cropland, and shrub), suggesting an important role of plant phenology in defining the statistical relationship between Tamax and LST. For impervious surfaces, the challenge of capturing the high spatial heterogeneity in urban settings using the low-resolution MODIS data made Tamax estimation a difficult task, which was especially true in the growing season.


Introduction
Surface air temperature (Ta) not only serves as a direct indicator of climate change [1], but also plays a crucial role in a range of studies for understanding the complex land surface processes, such as drought monitoring [2,3], greenhouse effect [4], hydrological simulation [5], urban heat island analysis [6][7][8][9], epidemiological modeling [10], and snowpack investigation [11].Spatially explicit Ta mapping has been traditionally relying spatial interpolation to scale up the point-based temperature measurements recorded by sample meteorological stations to the landscape level.
Remote Sens. 2017, 9, 233 2 of 14 However, the performance of such an approach is known to be highly dependent on the density of station network and the spatial heterogeneity of the geographic regions of interest [12].
Remote sensing provides a viable alternative that makes use of infrared sensors to measure land temperature through detecting thermal infrared radiation emitted from the geographic objects on land surfaces [13,14].The detected thermal radiation describes the kinetic energy of the ground objects informing land surface temperature (LST).While the LST is not equivalent to the air temperature above the surface (i.e., Ta), previous studies have confirmed a promising relationship between LST and Ta using three main categories of modeling approaches, including statistical, TVX (temperature-vegetation index), and energy-balance modeling.Specifically, (i) statistical modeling aims to establish an empirical relationship between Ta and LST through statistical regression analysis [15][16][17][18][19][20].To account for spatial heterogeneity in the land surface, ancillary data are included as supplementary explanatory variables, such as elevation, albedo, NDVI (normalized difference vegetation index), NDWI (normalized difference water index), and solar radiation [21][22][23][24][25][26][27][28][29][30]; (ii) The empirical TVX method relies on the negative linear relationship between LST and the abundance of live vegetation, as represented by a vegetation index, e.g., NDVI [31,32].The basic assumption is that Ta approaches LST under full canopy cover.Hence, Ta can be estimated by extrapolation of the LST-NDVI relationship to an NDVI value corresponding to the full canopy cover.However, recent studies have discovered inconsistent, non-linear relationships between LST and NDVI under complex land-cover conditions, which could negatively affect the TVX modeling performance (e.g., [33,34]); (iii) Energy-balance modeling is a physically-based approach directly grounded in thermodynamics.It treats Ta as the function of a range of variables describing the behavior of energy dynamics in the Earth system, such as radiation balance, and soil, sensible, and latent heat fluxes [35].While theoretically promising, it remains challenging to accurately estimate such a large number of model inputs solely through remote detection [16].
Among the three aforementioned categories of approaches, statistical modeling has been the most popular mainly due to its simplicity and relative robustness [36].However, few studies have systematically evaluated how the statistical relationship between Ta and LST varies across major land cover types.This is a crucial issue, because the Ta-LST relationship was found to be inconsistent over various landscapes, even if supplementary explanatory variables are used [24,37].Additionally, seasonal variation is likely to introduce extra uncertainties in Ta modeling.Especially at forested and agricultural sites, vegetation phenology plays a vital role affecting the actual land cover condition (e.g., leaf-on versus leaf-off).Previous efforts have discovered significantly varied modeling performance from the environments with low temperatures (or in winter seasons) versus those on hot summer days [38][39][40][41].To date, the joint impacts of land cover and seasonal variation remain to be systematically evaluated.
The main objective of this study was to improve our understanding of the statistical relationship between Ta and LST through remote detection.Here, we emphasized peak daytime air temperature (Tamax).Although the daytime surface energy flux is more complex than that at night, accurately estimating daytime air temperature was proven to more effectively inform extreme heat events and the associated human health issues (e.g., mortality) than using nighttime or daily average temperature [25,27,34,39,42].To achieve our goal, we ask two research questions in this study: (i) How do major land-cover types influence the statistical relationship between Ta and LST? (ii) How does seasonal variation affect such a relationship by changing land-cover conditions?

Study Area
Our study area is located in the middle and lower reaches of the Yangtze River (a.k.a., Changjiang River), the longest river in Asia and the third-longest in the world.The area is centered at the Poyang Lake Basin (Figure 1), which accounts for about 19% of the Yangtze River drainage area, with total coverage of 340,838 km 2 .Over the past decades, the region has experienced dramatic changes in landscape, resulting in variable land-cover types, significant environmental degradation, and substantial alteration in water and energy balance processes [43][44][45][46].The need for improved conservation of natural resources has stimulated intensive management and research activities in the area, e.g., monitoring crucial climatic variables (e.g., surface air temperature) from the recently established meteorological network [45,47].
Remote Sens. 2017, 9, 233 3 of 14 area, with total coverage of 340,838 km 2 .Over the past decades, the region has experienced dramatic changes in landscape, resulting in variable land-cover types, significant environmental degradation, and substantial alteration in water and energy balance processes [43][44][45][46].The need for improved conservation of natural resources has stimulated intensive management and research activities in the area, e.g., monitoring crucial climatic variables (e.g., surface air temperature) from the recently established meteorological network [45,47].

Air Temperature Meteorological Station Data
The meteorological time series of daily maximum air temperature observations from 1 December 2008 to 31 January 2010 were consistently collected by the China Meteorological Administration

Air Temperature Meteorological Station Data
The meteorological time series of daily maximum air temperature observations from 1 December 2008 to 31 January 2010 were consistently collected by the China Meteorological Administration (CMA).
We acquired data from 234 meteorological stations, which were evenly distributed throughout the study area (Figure 1).The locations of these stations can be identified on the MODIS (Moderate Resolution Imaging Spectroradiometer) LST images, then, LST data matching with daily maximum air temperature ground observations are extracted for study.Since MODIS LST products have a large amount of missing data due to cloud, heavy aerosols, only those observations matching with MODIS LST product data will be used for further study.Finally, our study used a total of 3916 observations, which were all taken at the height of 2 m above the ground level in accordance with the standards defined by the World Meteorological Organization (WMO).

MODIS Land Surface Temperature (LST) Data
The LST products from the MODIS system are one of the most widely used temperature products via remote detection.The system is comprised of sensors on two platforms of Terra (morning overpass at 10:30 am local) and Aqua (afternoon overpass at 1:30 pm local).Compared to Terra, Aqua has a closer overpass time when daily maximum air temperature occurs.However, previous efforts discovered superior performance of applying the LST data from MODIS/Terra to estimate air temperature [34,48].Hence, this study employed one of the Terra's daily land surface temperature and emissivity products, MOD11A1, in the succeeding analyses.MOD11A1 was produced at the spatial resolution of approximately 1 km (precisely 0.928 km), and gridded in the Sinusoidal projection.The MOD11A1 products were validated through field campaigns and radiance-based validation studies, with results indicating an accuracy of better than 1.0 K in most cases [49].
We downloaded the daily MOD11A1 products from the MODIS Land Processes Distributed Active Archive Center (LP DAAC) (http://lpdaac.usgs.gov/),with spatial and temporal coverage meeting the requirements of our research design.We further conducted data pre-processing by eliminating low-quality pixels (as indicated in the MOD11A1 Quality Control data layer) that corresponded to large errors induced by several factors, such as surface emissivity estimation error, atmospheric contamination, and viewing geometry on retrieval [49,50].We retained the 'good quality' pixels in L1B bands 31 and 32, with an average emissivity error of less than 0.01, and an average LST error of no more than 1.0 K.

Land-Cover Map
To evaluate the potential impact of land cover on the Ta-LST relationship, we adopted a relatively new land-cover product FROM-GLC (Finer Resolution Observation and Monitoring of Global Land Cover), provided by the Center for Earth System Science (CESS) at Tsinghua University [51,52].FROM-GLC was generated at the initial spatial resolution of 30 m using Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) data circa 2009 and 2010, which is consistent with the time period of our study.With the use of majority aggregation and proportion aggregation approaches, the products were delivered at multiple resolutions (e.g., 250 m, 500 m, 1 km, 5 km, and 10 km).To match the resolution of the selected MODIS LST product, the 1 km FROM-GLC map was chosen in this study, which was evaluated at an overall accuracy of 73.47% [52].Because we only focused on major land-cover types, the land-cover classes in FROM-GLC were further combined to represent forest, shrub, water, impervious surface, cropland, and grassland.To evaluate the accuracy of the final land-cover map at our study site, we further selected 200 randomly distributed sample points through field surveys and image interpretation, with the result reaching an accuracy of 82%.The percentage and area extent of each land-cover type are demonstrated in Table 1.

Model Development and Validation
Statistical analysis of the relationship between air and land surface temperature has been primarily relying on regression modeling in previous studies (e.g., [16,36,37,40,41]).Here, we adopted such approach using ordinary least squares regression (OLS) to establish the linkage between Tamax and LST.To account for the impacts of land cover and seasonal variation, models were developed for six land-cover types (i.e., forest, shrub, water, impervious surface, cropland, and grassland) and two main seasons (plant growing and non-growing seasons).For comparison purposes, we also built OLS models without distinguishing between land-cover types, and between seasons.In total, 19 models were developed, calibrated by a randomly selected 80% of the observational data from ground stations.
Model validation was conducted with the remaining 20% of the observational data.The performance of the models was evaluated using three criteria: adjusted coefficient of determination (R 2 adjusted ), mean absolute error (MAE), and root mean square error (RMSE).

Estimation of Maximum Air Temperature for Major Land-Cover Types
The scatter plots in Figures 2 and 3 illustrate the relationships between LST and the estimated Tamax across various land-cover types.Similar to previous studies [25,37,39,40], LST and Ta were found to have a relatively strong linear correlation.The general OLS model without considering land cover and seasonal variation reached an adjusted R 2 value of 0.87 (Figure 2).Compared to the general model, the models developed for specific land-cover types demonstrated varying levels of performance.More specifically, cropland, grassland, shrub, and water had better performance with higher adjusted R 2 values of 0.89, 0.92, 0.89, and 0.96 (Figure 3).Modeling forests and impervious surfaces showed inferior performance with lower adjusted R 2 values of 0.85 and 0.84 (Figure 3).The corresponding model validation results presented similar trends of performance using RMSE and MAE (Table 2).Distinguishing between land-cover types introduced an average of 4.5% of variance in adjusted R 2 .Particularly for water and grassland, LST was found to explain 92% and 96% of the variation in Tamax, indicating noticeable improvement in model performance for certain land-cover types.In contrast, forests were found to weaken the strong relationship between LST and Tamax.Compared to the general model, the forest model reduced the adjusted R 2 value by 2.0%, and considerably increased the model error by 15.2% in MAE and 19.6% in RMSE (Table 2), respectively.
Statistical significance tests showed that the effect of land cover on the Tamax-LST relationship was significant at the 0.01 level.The tests were performed by comparing the results predicted respectively using the general model and the land-cover specific models, using Paried-Samples T-TEST function in the SPSS software package (IBM, Armonk, NY, USA).

Model Development and Validation
Statistical analysis of the relationship between air and land surface temperature has been primarily relying on regression modeling in previous studies (e.g., [16,36,37,40,41]).Here, we adopted such approach using ordinary least squares regression (OLS) to establish the linkage between Tamax and LST.To account for the impacts of land cover and seasonal variation, models were developed for six land-cover types (i.e., forest, shrub, water, impervious surface, cropland, and grassland) and two main seasons (plant growing and non-growing seasons).For comparison purposes, we also built OLS models without distinguishing between land-cover types, and between seasons.In total, 19 models were developed, calibrated by a randomly selected 80% of the observational data from ground stations.
Model validation was conducted with the remaining 20% of the observational data.The performance of the models was evaluated using three criteria: adjusted coefficient of determination (R 2 adjusted), mean absolute error (MAE), and root mean square error (RMSE).

Estimation of Maximum Air Temperature for Major Land-Cover Types
The scatter plots in Figures 2 and 3 illustrate the relationships between LST and the estimated Tamax across various land-cover types.Similar to previous studies [25,37,39,40], LST and Ta were found to have a relatively strong linear correlation.The general OLS model without considering land cover and seasonal variation reached an adjusted R 2 value of 0.87 (Figure 2).Compared to the general model, the models developed for specific land-cover types demonstrated varying levels of performance.More specifically, cropland, grassland, shrub, and water had better performance with higher adjusted R 2 values of 0.89, 0.92, 0.89, and 0.96 (Figure 3).Modeling forests and impervious surfaces showed inferior performance with lower adjusted R 2 values of 0.85 and 0.84 (Figure 3).The corresponding model validation results presented similar trends of performance using RMSE and MAE (Table 2).Distinguishing between land-cover types introduced an average of 4.5% of variance in adjusted R 2 .Particularly for water and grassland, LST was found to explain 9%2 and 96% of the variation in Tamax, indicating noticeable improvement in model performance for certain land-cover types.In contrast, forests were found to weaken the strong relationship between LST and Tamax.Compared to the general model, the forest model reduced the adjusted R 2 value by 2.0%, and considerably increased the model error by 15.2% in MAE and 19.6% in RMSE (Table 2), respectively.
Statistical significance tests showed that the effect of land cover on the Tamax-LST relationship was significant at the 0.01 level.The tests were performed by comparing the results predicted respectively using the general model and the land-cover specific models, using Paried-Samples T-TEST function in the SPSS software package (IBM, Armonk, NY, USA).

Estimation of Maximum Air Temperature for Growing and Non-Growing Seasons
The scatter plots in Figure 4 illustrate the relationships between LST and the estimated daily maximum surface air temperature (Tamax) in the plant growing season across the six major land-cover types; while Figure 5 shows a similar type of scatter plots illustrating the Tamax-LST relationship in the non-growing season.Specifically, in the growing season, LST was found to explain 79, 72, 87, 75, 93, and 74 percent R 2 with Tamax across the six land-cover types of cropland, forest, grassland, shrub, water, and impervious surface, respectively.Compared to the growing season, most of non-growing season models revealed better performance with the corresponding adjusted R 2 values found to be 86, 80, 87, 86, 97, and 82 percent.The difference in model performance mainly occurred over the landscape covered by vegetation, e.g., cropland (7% increase), forest (8% increase), and shrub (9% increase).We also discovered a better model fitting over impervious surfaces (8% increase).However, seasonal variation did not demonstrate itself as an essential factor influencing the Tamax-LST relationship over water and grasslands, where the changes in adjusted R 2 were less than 5%.

Estimation of Maximum Air Temperature for Growing and Non-Growing Seasons
The scatter plots in Figure 4 illustrate the relationships between LST and the estimated daily maximum surface air temperature (Tamax) in the plant growing season across the six major land-cover types; while Figure 5 shows a similar type of scatter plots illustrating the Tamax-LST relationship in the non-growing season.Specifically, in the growing season, LST was found to explain 79, 72, 87, 75, 93, and 74 percent R 2 with Tamax across the six land-cover types of cropland, forest, grassland, shrub, water, and impervious surface, respectively.Compared to the growing season, most of non-growing season models revealed better performance with the corresponding adjusted R 2 values found to be 86, 80, 87, 86, 97, and 82 percent.The difference in model performance mainly occurred over the landscape covered by vegetation, e.g., cropland (7% increase), forest (8% increase), and shrub (9% increase).We also discovered a better model fitting over impervious surfaces (8% increase).However, seasonal variation did not demonstrate itself as an essential factor influencing the Tamax-LST relationship over water and grasslands, where the changes in adjusted R 2 were less than 5%.

Impact of Land Cover on the Estimation of Maximum Air Temperature
The modeling results for the six land-cover types confirmed the hypothesis that land surface characteristics have an impact on the Tamax-LST relationship.Similar analyses were carried out by Xu (2012) [24], Benali et al. (2012) [25], and Janatian et al. (2016) [30].The result in our study agrees with the findings of Xu ( 2012

Impact of Land Cover on the Estimation of Maximum Air Temperature
The modeling results for the six land-cover types confirmed the hypothesis that land surface characteristics have an impact on the Tamax-LST relationship.Similar analyses were carried out by Xu (2012) [24], Benali et al. (2012) [25], and Janatian et al. ( 2016) [30].The result in our study agrees with the findings of Xu ( 2012 2012) concluded that "model performance was independent of the meteorological stations' land cover".Janatian et al. ( 2016) indicated that a unique calibration equation can be adequately applied to a very large area (as large as half of Iran), despite an extreme heterogeneity of the land cover.We note that several auxiliary variables (e.g., elevation, zenith angle, reflectance, NDVI, etc.) were included in the Benali et al. (2012) and Janatian et al. (2016) models, which could have indirectly accounted for the effects by the land-cover variations.Another potential reason for contradiction is the differences in vegetation cover, and surface roughness, which may have influenced the surface-atmosphere energy exchanges.For example, Janatian et al. ( 2016) indicated that the land surface within the study area (eastern part of Iran) mainly consisted of bare or sparsely vegetated soils, which were different from our study area.In this study, such impact was found to vary across different types of terrestrial features.Particularly, estimating Tamax demonstrated superior performance over grasslands and water bodies (Figure 3).This may be explained by the fact that those landscapes are homogenous in surface cover, resulting in spatially homogenous energy exchange between land surfaces and the overlying atmosphere.
In contrast, modeling Tamax across forested areas showed noticeable inferior performance (Figure 3).This is possibly explained by the extensive heterogeneity in plant structure, live biomass, aerodynamic roughness, and leaf area index, which caused substantial variation of surface albedo in forests.Because our study area is located in the humid subtropical climatic region, highly mixed coniferous and broadleaved forest of rich species types (e.g.,China fir (Cunninghamia lanceolata), slash pine (Pinus elliotii), and Chinese red pine (Pinus massoniana)) and structures dominate the region [53][54][55], which further complicated the Tamax-LST relationship.Topography is another crucial factor.Our study area features rugged landscapes (e.g., karst), with steep slopes, deep valleys, and mixed mountains and plains.For example, 60% of the area of Jiangxi Province, centered at the study area, are mountains and hills [56].The elevation difference in some areas reaches over 2000 m [57].As mountainous areas are difficult to access, they are mainly covered by natural or planted forests (e.g., fruit trees and fast-growing timber trees).Hence, the combined heterogeneity in topography and forest structure possibly have negatively contributed to Tamax modeling in forested areas.

Impact of Seasonal Variation on the Estimation of Maximum Air Temperature
Plant surfaces typically show distinct characteristics in growing versus non-growing seasons corresponding to their phenological development.This changed plant surface albedo, and in turn affected the LST-Tamax relationship as revealed by the study (Figures 4 and 5).Their variation in adjusted R 2 was found to be 8% and 6% for the growing and non-growing season, respectively.However, such variation was 4% without considering seasonal variation (Figure 3).
Our results further indicated that the errors of estimating Tamax likely increase during the growing season for the landscapes covered by crops, forests, and shrubs (Figures 4 and 5).However, their models all achieved better performance in the non-growing season, with an average increase of 8% in adjusted R 2 .This was possibly owing to the subtropical climate of the study area, featuring a wide range of vegetation types.For example, local crops include cereal, wheat, cotton, peanut, tea, and tobacco.These areas demonstrated a sharp contrast in vegetation cover between growing and non-growing seasons.Forests in the study area consist of both natural and planted trees, such as those producing non-timber products of pineapple, oak, sugarcane, tea, and orange [58,59].The 1 km resolution data was coarse, which did not allow us to capture such heterogeneity in vegetation type and structure.From the phenology perspective, those plants exhibited various forms of growth during the growing season.In one of our earlier evaluations, we tried to divide the long growing season into several sub-seasons; however, the performance of the developed models exhibited high variation with no clear patterns.Such variation led to a high level of heterogeneity in surface characteristics in the growing season, as compared to the landscapes with reduced vegetation cover in the non-growing season.
Besides the vegetation-covered landscapes, impervious surfaces also demonstrated noticeable model improvement using data from the non-growing season, with results showing an increase of 8% in adjusted R 2 (Figures 4 and 5).This may be explained from two aspects.First, impervious surfaces were primarily located in urban areas, interspersed with vegetation within the city bounds.At the spatial resolution of 1 km, a land patch classified as impervious surface possibly contained plants, where their spectral variation in the non-growing season was much lower than that in the growing season as described in the previous paragraph [60].Second, the growing season in the study site is part of the year with a significant increase of precipitation.Researchers have argued that topography and land cover can lead to substantial variability of rainfall both spatially and temporally [61,62], which adds variation to land-atmosphere interactions across the urban areas.
Our analysis was based on a linear relationship between Tamax and LST, although several scatter plots of Figures 2-5 indicate a 'convex upward' distribution.This suggests a potential nonlinear quadratic relationship between those two variables, inferring a potential effect of within-season variation on Tamax estimation.It can be further observed that the regression lines for the non-growing season (Figure 5) are steeper than those for the growing season (Figure 4).That is, Tamax rises more rapidly when LST increases in colder seasons than that in warmer seasons.However, our emphasis was the model comparison among various land-cover types, and the nonlinear quadratic relationship only showed slightly improved model performance.We decided to use linear models in this study.

Spatial Uncertainties in the Estimation of Maximum Air Temperature
Since land cover and seasonal variation were found to affect the estimation of Tamax, we have applied the previously developed models to produce two temperature maps considering major land-cover types for two dates in the growing season (on 24 May 2010; Figure 6a) and non-growing season (on 31 December 2010; Figure 7a), respectively.Figure 6a reveals a clear temperature gradient in spring, i.e., the closer to water bodies, the lower the maximum air temperature.For example, Fujian Province in the southeast faces East China Sea to the east, and South China Sea to the south, and its estimated Tamax was much lower than that in Hubei Province, located in the Central China region.Besides ocean, inland water (e.g., the Poyang Lake in Jiangxi Province) also showed an evident impact in cooling down the neighborhood regions in the upper middle of our study area.In contrast, the spatial pattern of winter temperature was mainly controlled by the latitudinal variation that has more than five degrees across out study area (Figure 7a).
Remote Sens. 2017, 9, 233 10 of 14 Our analysis was based on a linear relationship between Tamax and LST, although several scatter plots of Figures 2-5 indicate a 'convex upward' distribution.This suggests a potential nonlinear quadratic relationship between those two variables, inferring a potential effect of within-season variation on Tamax estimation.It can be further observed that the regression lines for the non-growing season (Figure 5) are steeper than those for the growing season (Figure 4).That is, Tamax rises more rapidly when LST increases in colder seasons than that in warmer seasons.However, our emphasis was the model comparison among various land-cover types, and the nonlinear quadratic relationship only showed slightly improved model performance.We decided to use linear models in this study.

Spatial Uncertainties in the Estimation of Maximum Air Temperature
Since land cover and seasonal variation were found to affect the estimation of Tamax, we have applied the previously developed models to produce two temperature maps considering major land-cover types for two dates in the growing season (on 24 May 2010; Figure 6a) and non-growing season (on 31 December 2010; Figure 7a), respectively.Figure 6a reveals a clear temperature gradient in spring, i.e., the closer to water bodies, the lower the maximum air temperature.For example, Fujian Province in the southeast faces East China Sea to the east, and South China Sea to the south, and its estimated Tamax was much lower than that in Hubei Province, located in the Central China region.Besides ocean, inland water (e.g., the Poyang Lake in Jiangxi Province) also showed an evident impact in cooling down the neighborhood regions in the upper middle of our study area.In contrast, the spatial pattern of winter temperature was mainly controlled by the latitudinal variation that has more than five degrees across out study area (Figure 7a).To better understand the spatial uncertainty of our estimates, the same-day field observations were used as a comparison, where the radius of each circle corresponds to the level of error at individual meteorological stations (Figures 6b and 7b).The MAE in Figures 6b and 7b ranges from 2.46 to 3.00 and from 1.34 to 1.60 respectively.The land-cover map was used as a backdrop.The circles in the two figures demonstrated consistent findings, where the low-performance stations were mostly located in vegetated areas, especially those vegetated areas in the mountains.In addition to the aforementioned reason of high level of heterogeneity in species types and topography, a portion of errors may stem from the data acquisition stage.For example, as pointed out by Lee et al. (2011) [63], meteorological stations that represent the surrounding forests are in fact typically located in forest To better understand the spatial uncertainty of our estimates, the same-day field observations were used as a comparison, where the radius of each circle corresponds to the level of error at individual meteorological stations (Figures 6b and 7b).The MAE in Figures 6b and 7b ranges from 2.46 to 3.00 and from 1.34 to 1.60 respectively.The land-cover map was used as a backdrop.The circles in the two figures demonstrated consistent findings, where the low-performance stations were mostly located in vegetated areas, especially those vegetated areas in the mountains.In addition to the aforementioned reason of high level of heterogeneity in species types and topography, a portion of errors may stem from the data acquisition stage.For example, as pointed out by Lee et al. (2011) [63], meteorological stations that represent the surrounding forests are in fact typically located in forest clearings, which makes the observations less representative with high variation.Another source of error is possibly from land cover classification.Benali et al. (2012) [25] argued that land cover data with different accuracies and from various classification schemes could result in high uncertainties in temperature estimation.To mitigate the impact of error propagation from a complex classification scheme, a relatively simple and straightforward classification system may be a better option [49].In addition, we suggest that such a system needs to incorporate vegetation phenology in classification, so the results can better assist with our assessment of the Tamax-LST relationship in the plant growing season.
clearings, which makes the observations less representative with high variation.Another source of error is possibly from land cover classification.Benali et al. (2012) [25] argued that land cover data with different accuracies and from various classification schemes could result in high uncertainties in temperature estimation.To mitigate the impact of error propagation from a complex classification scheme, a relatively simple and straightforward classification system may be a better option [49].In addition, we suggest that such a system needs to incorporate vegetation phenology in classification, so the results can better assist with our assessment of the Tamax-LST relationship in the plant growing season.

Conclusions
Results in this study revealed noticeable impacts of land cover and seasonal variation on the estimation of Tamax using MODIS-derived LST.Owing to the high level of heterogeneity in plant structure and topography, the Tamax-LST relationship was found to be negatively affected in forest landscapes with high uncertainties.Estimating Tamax over grasslands and water bodies achieved superior performance, which could be explained by the spatially homogenous energy exchange between land surfaces and the overlying atmosphere.We also found that all the land-cover specific models developed for the non-growing season outperformed the corresponding models developed for the growing season.Discrepancies in model performance mainly occurred in the vegetated areas (forest, cropland, and shrub), suggesting that plant phenology plays an important role defining the statistical relationship between Tamax and LST.Because different plant species exhibit varying growth patterns over time, understanding specific species types may significantly improve model performance in the growing season.Compared to the other land-cover types, impervious surfaces in urban areas contained high spatial heterogeneity that was difficult to capture by MODIS, which led to high uncertainties in Tamax estimation.In contrast to previous studies, this research provides a more in-depth analysis emphasizing on major land-cover types and seasons in the middle-low reaches of Yangtze River featuring high variation in topography and climatic conditions.Our research suggests the necessity of improving the current Tamax-LST relationship through modeling the seasonal, sub-pixel landscape heterogeneity that typically relies on medium to coarse spatial resolution Earth observation data, although we take advantage of the sensors' short revisit intervals.

Conclusions
Results in this study revealed noticeable impacts of land cover and seasonal variation on the estimation of Tamax using MODIS-derived LST.Owing to the high level of heterogeneity in plant structure and topography, the Tamax-LST relationship was found to be negatively affected in forest landscapes with high uncertainties.Estimating Tamax over grasslands and water bodies achieved superior performance, which could be explained by the spatially homogenous energy exchange between land surfaces and the overlying atmosphere.We also found that all the land-cover specific models developed for the non-growing season outperformed the corresponding models developed for the growing season.Discrepancies in model performance mainly occurred in the vegetated areas (forest, cropland, and shrub), suggesting that plant phenology plays an important role defining the statistical relationship between Tamax and LST.Because different plant species exhibit varying growth patterns over time, understanding specific species types may significantly improve model performance in the growing season.Compared to the other land-cover types, impervious surfaces in urban areas contained high spatial heterogeneity that was difficult to capture by MODIS, which led to high uncertainties in Tamax estimation.In contrast to previous studies, this research provides a more in-depth analysis emphasizing on major land-cover types and seasons in the middle-low reaches of Yangtze River featuring high variation in topography and climatic conditions.Our research suggests the necessity of improving the current Tamax-LST relationship through modeling the seasonal, sub-pixel landscape heterogeneity that typically relies on medium to coarse spatial resolution Earth observation data, although we take advantage of the sensors' short revisit intervals.

Figure 1 .
Figure 1.Maps of the study area (a) Location map; (b) Land cover and the spatial distribution of meteorological stations map; (c) Topographic map.

Figure 1 .
Figure 1.Maps of the study area (a) Location map; (b) Land cover and the spatial distribution of meteorological stations map; (c) Topographic map.

Figure 2 .
Figure 2. Scatter plot of the daily maximum air temperature (Tamax) and surface temperature (LST) relationship without distinguishing between land covers.

Figure 2 .
Figure 2. Scatter plot of the daily maximum air temperature (Tamax) and surface temperature (LST) relationship without distinguishing between land covers.

Figure 6 .
Figure 6.Spatial patterns of the estimated Tamax (a) on 24 May 2010, and the corresponding uncertainties at the stations (b).The level of uncertainty is represented by the size of circle.

Figure 6 .
Figure 6.Spatial patterns of the estimated Tamax (a) on 24 May 2010, and the corresponding uncertainties at the stations (b).The level of uncertainty is represented by the size of circle.

Figure 7 .
Figure 7. Spatial patterns of the estimated Tamax (a) on 31 December 2010, and the corresponding uncertainties at the stations (b).The level of uncertainty is represented by the size of circle.

Figure 7 .
Figure 7. Spatial patterns of the estimated Tamax (a) on 31 December 2010, and the corresponding uncertainties at the stations (b).The level of uncertainty is represented by the size of circle.

Table 1 .
Areal coverage of individual major land cover types in the study area.

Table 2 .
Model performance for the combined and individual land-cover types.

Table 2 .
Model performance for the combined and individual land-cover types.