Global Land Surface Temperature Change (2003–2017) and Its Relationship with Climate Drivers: AIRS, MODIS, and ERA5-Land Based Analysis

Land surface temperature (LST) plays a critical role in the water cycle and energy balance at global and regional scales. Large-scale LST estimates can be obtained from satellite observations and reanalysis data. In this study, we first investigate the long-term changes of LST during 2003–2017 on a per-pixel basis using three different datasets derived from (i) the Atmospheric Infrared Sounder (AIRS) onboard Aqua satellite, (ii) the Moderate Resolution Imaging Spectroradiometer (MODIS) also aboard Aqua, and (iii) the recently released ERA5-Land reanalysis data. It was found that the spatio-temporal patterns of these data agree very well. All three products globally showed an uptrend in the annual average LST during 2003–2017 but with considerable spatial variations. The strongest increase was found over the region north of 45° N, particularly over Asian Russia, whereas a slight decrease was observed over Australia. The regression analysis indicated that precipitation (P), incoming surface solar radiation (SW↓), and incoming surface longwave radiation (LW↓) can together explain the inter-annual LST variations over most regions, except over tropical forests, where the inter-annual LST variation is low. Spatially, the LST changes during 2003–2017 over the region north of 45° N were mainly influenced by LW↓, while P and SW↓ played a more important role over other regions. A detailed look at Asian Russia and the Amazon rainforest at a monthly time scale showed that warming in Asian Russia is dominated by LST increases in February–April, which are closely related with the simultaneously increasing LW↓ and clouds. Over the southern Amazon, the most apparent LST increase is found in the dry season (August–September), primarily affected by decreasing P. In addition, increasing SW↓ associated with decreasing atmospheric aerosols was another factor found to cause LST increases. This study shows a high level of consistency among LST trends derived from satellite and reanalysis products, thus providing more robust characteristics of the spatio-temporal LST changes during 2003–2017. Furthermore, the major climatic drivers of LST changes during 2003–2017 were identified over different regions, which might help us predict the LST in response to changing climate in the future.


Introduction
Land surface temperature (LST) can be recognized as the heat emissions on the Earth's surface and refers to the radiative skin temperature on the land surface [1]. LST is a vital indicator for studying environmental changes, hydrological conditions, and the Earth's energy balance [2][3][4]. LST measurements can be obtained from (i) the ground observations of in situ stations or networks [5], (ii) retrievals of remote sensing data [6], and (iii) reanalysis data based on surface energy balance model estimations [7]. Although in situ measurements have high accuracy and temporal resolution, their point-scale representativeness and the sparse distribution of meteorological stations are major limitations for both their research and applications at continental to global scales [8]. Furthermore, the accuracy of interpolated LST products is based on the density of ground observations. future. Section 2 introduces the various datasets used. Section 3 introduces the methods employed to analyze those datasets. In Section 4, the results of the analysis are presented, and Section 5 provides discussions of the implications of the results, as well as the added value of the study with respect to existing studies. Finally, in Section 6, the conclusions of the study are presented.
Both MODIS and AIRS LST products used in this study were obtained from the Aqua satellite with local equatorial crossing times around 1:30 a.m. (descending) and 1:30 p.m. (ascending).
The MODIS instrument has 36 channels with a swath of 2330 km, including 20 visible and near-infrared bands and 16 thermal infrared bands ranging from 3 to 15 µm [38]. Eight thermal infrared bands belonging to the atmospheric window channels were designed to retrieve land surface temperatures via the day/night algorithm [39] and split-window algorithm [40]. New refinements were applied to the above two algorithms, and radiancebased validation of MODIS collection 6 (C6) was implemented at five bare soil sites in North Africa [41]. The comparison studies indicated that the improved C6 is better than C5, especially in bare areas, with errors of less than 1 K at most sites.
A swath of 1650 km and 2378 channels of thermal infrared are characteristics of AIRS that distinguish it from MODIS [42]. Experiments to evaluate the uncertainty of AIRS LST products were carried out, and the results showed that the root mean square error (RMSE) of the daytime data was between 1.5 and 3 K in different sites and had a positive bias of longwave spectral emissivity [43]. Neural network technology was applied to the latest retrieval algorithm of AIRS Version-6-more accurate initial emissivity predictions and a reduction in the impact of cloud errors are some significant advancements of this method [44]. The AIRS-Only retrievals have comparable performance to AIRS/Advanced Microwave Sounding Unit (AMSU) retrievals [45] but with longer time series publicly available. Hence, the 1 • AIRS-Only monthly LST dataset (AIRS3STM) without AMSU observations was utilized in this study.
The maximum temperature is more suitable than the mean daily average temperature for assessing thermal changes in large-scale climate systems, such as defining temperature trends [46]. Moreover, the maximum daily LST is very sensitive to the energy changes of the Earth [47]. We, therefore, used the MODIS and AIRS LST retrievals from the afternoon overpass (1:30 PM), which is very close to the maximum temperature during a day.
The reanalysis LST in this study was from ECMWF's ERA5-Land's skin temperature (hereafter ERA5-Land LST) database based on the HTESSEL model. This database uses the CY45r1 version of the Integrated Forecasting System. ERA5-Land uses an hourly output frequency, and data are available from 1981 to present. To match the satellite LST retrievals, for each grid cell, we used the ERA5-Land LST from the hour corresponding to the maximum temperature within a day for the climatology time series of ERA5-Land. Johannsen et al. [48] evaluated the ERA5 LST using the Satellite Application Facility on Land Surface Analysis over the Iberian Peninsula, and found that the standard deviation error of the ERA5 daily maximum LST is around 2-3 K.
The long-term monthly precipitation datasets of the Climatic Research Unit (CRU) were used in this study. Evaluation studies have shown the reliability of CRU precipitation products both globally and regionally for both in the long and short term [49,50]. Along these lines, CRU precipitation has been used in previous studies for trend analysis, characterizing the spatio-temporal distribution of precipitation and assessing the impact of precipitation on other climate variables [51][52][53]. Similarly, the radiation variables of the Clouds and Earth's Radiant Energy System (EBAF-Surface Ed4.0) were verified in detail in, for instance, the study of Loeb et al. [54]. Furthermore, defining the trends of radiation allows us to investigate the Earth's climate from a thermodynamic point of view [55]. This study also used AIRS CF, MODIS AOD, and AIRS CO 2 data to define atmospheric conditions. The works of Wu et al. [56] shows the AIRS cloud property retrievals to agree with Lidar observations. Moreover, the MODIS AODs were in good agreement with the sunphotometer measurements, with a correlation coefficient of 0.85 and an RMSE of 0.15 in the Yangtze River Delta in China [57]. The accuracy of AIRS CO 2 is 1-2 ppm, according to comparative work [58]. Regarding the surface conditions, the LAI and land cover data from MODIS were verified and compared at different scales [59,60] and can describe the vegetation and surface types to a certain extent. The above data are valuable to study environmental changes.
All datasets were downloaded at monthly intervals from January 2003 to December 2017, except for CO 2 , LAI, and land cover. Monthly CO 2 datasets are only available up to February 2017, so we only used CO 2 data from January 2003 to February 2017. The LAI data used here have a temporal resolution of 8 days and cover the period of 2003-2017. The land cover data have an annual temporal resolution. To enable a direct spatial comparison, relatively high-resolution products, such as MODIS LST, ERA5-Land LST, CRU P, and MODIS LAI products, were converted into 1 • data using the aggregation method, and the coarse-resolution AIRS CO 2 data were resampled as 1 • data via bilinear interpolation from their original resolutions. To convert land cover data from a 0.05 • to 1 • spatial resolution, the most common land cover type within each 1 • grid cell was used.

Identification of the Spatio-Temporal Consistency of LST Trends Derived from AIRS, MODIS, and ERA5-Land
First, for each 1 • grid cell over the land, we calculated the annual LST from 2003 through 2017 by taking the average of the monthly data from the AIRS, MODIS, and ERA5-Land products. Then, we conducted a trend analysis retrieved by linear regression during 2003-2017 over each 1 • grid cell and obtained the spatial distribution of the linear trend of annual LST from the AIRS, MODIS, and ERA5-Land products. Two statistical methods were used to examine their consistency: (1) a scatter plot of the spatial distribution of the linear trend was made for any two of these three LST products, and the corresponding correlation coefficient (R) was calculated. (2) For each continent, the changes in LST were summarized as boxplots based on the 5th, 25th, 50th, 75th, and 95th percentile values. The boxplots derived from AIRS, MODIS, and ERA5-Land products were also compared.

Identification of the Spatio-Temporal Characteristics of Climatic Variables Associated with LST Trends
At the global scale, the trend analysis was first performed on the annual P, SW↓, and LW↓ during 2003-2017 over each 1 • grid cell. Then, we applied the linear regression method [61] to identify the coefficient of determination (R 2 ) between the annual LST and these three climatic variables (i.e., annual P, SW↓, and LW↓), both individually and combined, during 2003-2017 over each 1 • grid cell. Based on the global results, regions with different major climatic drivers were selected for more detailed temporal analyses.
Over these selected regions, the trend analysis was first conducted on the regional average LST for each month from January to December during 2003-2017 to identify the months with the strongest changes in LST. During these months, we examined the time series of regional average LST, P, SW↓, LW↓, cloud fraction (CF), and aerosol optical depth (AOD) during 2003-2017 by calculating the linear trends of each variable and the correlation coefficient (R) among these variables [62].

Identification of the Spatio-Temporal Consistency of LST Trends Derived from AIRS, MODIS, and ERA5-Land
The spatial distributions of trends in annual LST during 2003-2017 derived from AIRS, MODIS, and ERA5-Land are shown in Figure 1a-c, respectively. Their spatial patterns appear similar, with most regions showing increases in LST trends. Quantitatively, the agreement in the spatial distribution of the LST trends between the two satellite products (AIRS and MODIS) was shown to be very high, with the correlation coefficient (R) value reaching 0.92 (Figure 1d). This R value was 0.74 between AIRS and ERA5-Land ( Figure 1e) and 0.73 between MODIS and ERA5-Land (Figure 1f). Apparent differences between the satellite and reanalysis LST trends were observed over eastern Brazil, where the increasing trends in LST were stronger in the satellite retrievals, as well as over central Africa and northern Australia, where the positive trends in LST were stronger in the reanalysis product. The similar result from these two satellites was mainly due to the similarities in their retrieval process and the fact that they are both thermal infraredbased observations. Thus, similar limiting factors would affect their qualities. Conversely, the difference between satellite and reanalysis datasets may be the result of a number of factors which include the land model parameterizations and sensitivities.
The boxplots shown in Figure 1g compare the trends of LST derived from the different products over different regions. Globally, the average change rate of the LST derived from AIRS, MODIS, and ERA5-Land was about 0.03, 0.02, and 0.04 K year −1 , respectively. For the continents, larger increases in rates of change greater than 0.1 K year −1 were observed over North America, Europe, and Russia (particularly over Asian Russia), whereas Oceania experienced an overall decreasing trend in LST, possibly due to Australia's contribution during 2003-2017.
According to Figure 2a,b, the overall satellite data and reanalysis data trends differ little over the global land (within 0.05 K/year). However, the satellite results were clearly underestimated over central Africa and overestimated over northeastern Brazil. Based on the spatial pattern of the trends in annual LAI (See Figure 2c), the positive bias of the trends is consistent with a decline in LAI, and the negative deviation corresponds to an increase in LAI. In northeastern Brazil, South America is mostly shown to be warming, and vegetation is also rapidly decreasing. This decline in LAI and the appearance of LST warming may indicate that dramatic environmental changes are being experienced in northeastern Brazil.

Identification of the Spatio-Temporal Characteristics of Climatic Variables Associated with LST Trends
The spatial distributions of trends in annual P, SW↓, and LW↓ during 2003-2017 are shown in Figure 3a-c. Large negative trends of P are visible in the tropical rainforests and savanna of South America, while SW↓ increased dramatically in central South America. Thus, a distinct negative correlation between SW↓ and P can be observed. In some areas of high latitudes, the opposite relationship is shown-that is, P increases while SW↓ decreases. LW↓ is overall on the rise globally, but is decreasing in some regions, such as the tropical forests of Brazil, northern Australia, and Central Asia. Among them, tropical forests in South America showed the most obvious pattern of decreasing trends in LW↓.  Figure 3e-g. From these three analyses, we examined the different driving factors in the northern and southern hemispheres. The P and SW↓ largely explain the LST changes in the southern hemisphere (such as the south of Africa, the east and south of South America, and most of Australia), while LW↓ can explain the high latitudes in the northern hemisphere. Figure 3h also shows the results of the stepwise regression. Using the high determination coefficient (R 2 ) of the model based on stepwise multiple linear regression of the three variables, we determined that the three variables can largely explain (R 2 > 0.8) the changes in LST over the major global land areas. However, regionally insignificant areas were also found in the tropical forests and barren land in Africa and South America (Figure 3d).
Based on the results, the distinct directions of the LST trends in northeast Brazil and Australia can be clearly explained by changes in P and SW↓. Reduced P and increased SW↓ were found to contribute to dramatic warming in northeastern Brazil, while the opposite trend was found over Australia. The decline in P and LW↓ in some areas of the warming Amazon, with an LST rate of change of about 0.05 K/year, exceeded −3 mm/year and −1.2 W m −2 /year, respectively, and the increase in SW↓ also exceeded 0.8 W m −2 /year, as shown in Figure 3a-c. The P, SW↓, and LW↓ all changed drastically; however, the results of the regression show an insignificant R 2 , suggesting that none of these factors are major drivers at inter-annual scales in the Amazon. For Asian Russia, intense warming was accompanied by an increase in P and, at the same time, a significant effect of LW↓ based on the R 2 over 0.8 between LST and LW↓. The result dominated by LW↓ also attracted our attention in Asian Russia. The increase in downward radiation is closely related to changes in atmospheric conditions, and more variables such as clouds, AOD, and CO 2 need to be considered for analysis. It is thus worth exploring the LST changes and potential factors based on a more detailed temporal analysis in the two areas. The black boxes in Figure 3h show the two regions that we focused on (Asian Russia: 50 • N-80 • N, 50 • E-180 • E; the Amazon: 0 • -15 • S, 75 • W-52 • W).
We computed the LST trends of the zonal averages (Asian Russia and Amazon) on a monthly basis during 2003-2017. The LST trends (left axis) and LST seasonal cycles (right axis) are shown in Figure 4 for the three sets of LST products. Regarding the LST trend and seasonal cycles, the three sets of data in the two regions all present highly consistent results in the peak months, with relatively large positive trends according to the red box in Figure 4, while differences in the LST trends in the Asian Russia results can be found for May and July between the satellite data and reanalysis data. In Asian Russia, the LST warmed abnormally (up to about 0.2 K/year) in February, March, and April, which are generally colder months within the seasonal cycle. Over the Amazon, positive trends occurred in most months, although the values were small (<0.1 K/year), confirming that this region, overall, experiences warming all year round. Based on the above analysis, Asian Russia's February to April and the Amazon's August to September data were selected as further study months to explore the relationship between LST changes and various climate drivers.
The P, SW↓, LW↓, CF, and AOD were taken into consideration to analyze the LST changes, for which the selected months were averaged for each year and compared with the LST variabilities for the same period ( Figure 5). The time series changes are analyzed together with the correlation coefficients in Table 2.
In Asian Russia, the three sets of LST showed similar variability, although their absolute values varied (Figure 5a). At the same time, the significant warming trends are very obvious (e.g., the rate of change of AIRS LST is 0.21 K/year), especially in recent years. In Figure 3e,f, the spatial patterns of P and SW↓ were found not to be significantly linked to the warming trends observed in the LST over Russia. Since the increased P of 0.17 mm/year and the decreased SW↓ of −0.21 W/ m −2 could not be the drivers increasing the LST, they were further determined not to be warming factors over Asian Russia. On the other hand, the rate of change of LW↓ reached 0.84 W/ m −2 /year and has recently been increasing year by year, as shown in Figure 5d. It is obvious that the variabilities of LST and LW↓ are consistent with each other, reaching up to a significant positive correlation (R = 0.75), as seen in Table 2. According to Figure 5e,f, the rates of change in CF reaching 0.03/year, and the slightly decreased AOD (−0.005/year), are evidence of increased LW↓, with the correlation between CF and LW↓ reaching as high as 0.79. The relationship between CF and LST is also very clear, reaching a high correlation coefficient of 0.89. Moreover, the increase in P (R = 0.64 between P and CF) and the decrease in SW↓ (R = −0.44 between SW↓ and CF) can also be reasonably explained by the change in CF.   . This special dominant factor mechanism indicates that there may be no direct connection between AOD and LST dynamics. It is worth noting that the change of LW↓ here was found to be related to AOD to a certain extent (R = 0.47 between AOD and LW↓). The low AOD values of 2011-2014 correspond to a trough of LW↓; the increase in LW↓ in the years that followed was also accompanied by a slight increase in AOD. Therefore, AOD plays an important role in the local energy balance and further affects changes in the LST.
As a greenhouse gas, CO 2 can cause environmental warming [63]. Because the time period of the CO 2 data used in this study does not match that of the previous analysis, CO 2 is analyzed separately, and the trends are inter-annual and seasonal. The inter-annual and seasonal CO 2 concentration trends are increasing on a global scale (see Figure 6) and not just anomalously over certain areas. Winter in Asian Russia has increased relatively significantly based on Figure 6b, which provides results obtained from AIRS satellite data from 2003-2016. Therefore, meteorological factors may be the most significant factors that affect changes in LST after considering the previous analysis. The continuous increase in CO 2 is also worthy of attention.

The Pattern of Spatio-Temporal LST Trends Derived from AIRS, MODIS, and ERA5-Land
This study observed that both the satellite and reanalysis data show consistent LST changes. Anomalous rises in temperature were mainly found to appear in Asian Russia and parts of South America, while anomalous temperature decreases mainly appeared in Australia. The overall results of these signals are consistent with previous studies using four daily representative observations of MODIS from Aqua and Terra satellites to evaluate the global LST changes during 2003-2016 [64]. The greatest fluctuations in warming and cooling on land occurred in Russia and eastern Australia, respectively, based on an analysis of global MODIS LST changes during 2001-2012 [65]. Zhou et al. [66] demonstrated that the trend of global skin temperature over land was approximately 0.036 K/year during June 2007 to April 2018 based on the Interferometer Atmospheric Sounding Instrument (IASI) used on the Meteorological Operational Satellites; their results also found that the decadal warming trend is more pronounced in the northern hemisphere. In addition, the fastest warming in recent decades (> 0.4 K/decade) was observed in the northern mid-latitude region based on land surface air temperature data [67].
It is worth noting that the satellite and ERA5-Land data show differences in the warming values in eastern South America and local opposite trends over Africa (see Figure 1). The availability of satellite observations is highly dependent on cloudiness, which caused differences between satellite and reanalysis LST trends over Antarctica [16], demonstrating that cloud detection uncertainty is still a hindrance to infrared remote sensing [68]. Moreover, uncertainty may be caused by the lack of a full understanding of the land surface parameters [11]. Furthermore, the use of satellite vegetation products in the assimilation system was shown to have a positive impact on near-real-time issues [69]. In addition, Zhou et al. [17] found that in about 70% of global desert areas, LST trends in the reanalysis products are generally underestimated, which suggested a limitation of the performance of the reanalysis products in reproducing LST dynamics. Although ERA5-Land does not consider the inter-annual dynamics of vegetation, the correlation coefficient is 0.74 between AIRS and ERA5-Land and 0.73 between MODIS and ERA5-Land (see Figure 1e,f), which illustrates that climatic factors play a very important role in LST changes.
Deforestation also caused a significant warming in the average temperature trends in tropical regions due to the effects of reduced evapotranspiration on temperature [70]. Recently, the forests and Cerrado have been browning in eastern South America [71], which may be related to deforestation expansion. Therefore, a decline in vegetation is also related to a dramatic increase in LST. For this reason, the regulation of vegetation on LST will contribute to monitoring the forest conservation.

Potential Climatic Variables Associated with LST Trends
Based on the spatial patterns of global changes in LST and the other climatic variables, including precipitation, shortwave radiation, and longwave radiation, the link between LST and these variables was assessed. Based on the regression method, radiation and precipitation can explain the LST changes of most global land areas at inter-annual scales where the differences can also be observed along the latitudes. Longwave radiation mainly affects the LST at high latitudes over the northern hemisphere, while shortwave radiation forcing and precipitation can largely explain the LST changes in the rest of the northern hemisphere and southern hemisphere. Radiation was also found to be the main factor that explains the inter-annual variability of the monthly average daily maximum temperature over land [72]. Therefore, enhanced radiative forcing is an important factor affecting the LST, which is consistent with our conclusions.
Additionally, Song et al. [73] studied the effects of vegetation and fine particulate matter (PM2.5) on MODIS LST changes globally during 2001-2016. The authors found that the relationship between LST and vegetation in the northern hemisphere is, overall, positively correlated, while that of the southern hemisphere is negatively correlated. Over typical areas in eastern South America and the Mediterranean region, our results also clearly reflect this rule based on trends in annual LST and LAI. It was found that PM2.5 affects radiative forcing, causing LST to increase with a small increase in PM2.5 concentrations but to decrease with a marked increase in PM2.5. As an important variable of human activities affecting the environment, PM2.5 is worth considering, to study its impact on LST changes. More detailed analyses are discussed in the section on regional changes in LST and climatic variables.
Interestingly, the increased precipitation in Asian Russia did not result in cooling effects similar to those in Australia. Australia experienced extreme precipitation in 2010/2011 via La Niña events [74], which explains Australia's cooling trend. According to our results, in Asian Russia, the close relationship between LST, longwave radiation, and clouds agrees well with the findings of Jing et al. [75], in which the longwave scattering of clouds increased downward with a longwave radiation flux, thereby increasing surface temperature. Longwave radiation was affected by the variation in cloud cover during the research period of this paper. Based on a large number of daytime observations from Russia's meteorological station, cloud cover increased significantly in the spring of 2001-2010 compared to the spring of 1990-2000 [76]. In the Arctic, there is less solar radiation in winter and spring, and anthropogenic pollution is relatively high, thereby increasing the emissivity of the clouds, which increases the longwave radiation [77]. The warming of the northern hemisphere's cryosphere is related to the radiative surface albedo feedback caused by changes in ice cover [78,79]. Bekryaev et al. [80] found that the longwave radiative forcing related to cloud and water vapor content in winter and spring is the reason for polar feedback amplification, which cannot be explained solely by the ice-albedo feedback mechanisms. As temperature rises, convective motion increases [75]; precipitation and water vapor are also increasing in the Arctic Circle [81]. As humidity and clouds in the atmosphere of the Arctic increase, longwave radiation towards the ground likewise increases, contributing to the warming [82]. Therefore, widespread and seasonal warming is accompanied by a feedback mechanism. Warming has brought about the melting of ice and snow, an unstable atmosphere, and more water vapor, clouds, and rainy weather. The current discovery of increased vegetation in the northern hemisphere was also explained by rising temperatures and increased moisture. The reduction in surface albedo and the increase in the greenhouse effect are exacerbating warming. In this process, clouds play a key role. The warming of permafrost in Russia has the potential to influence global climate change because the thawing of frozen soil releases soil organic carbon [83]. At longer timescales in the past, Russia's winter became warmer, and the greenhouse gas effect has been considered as an important factor for this warming [84]. Based on the results of the AIRS CO 2 data, CO 2 continues to increase globally. At the time scale studied in this paper, CO 2 may not be the main influencing factor for the most intense warming in Asian Russia. With such a large area of warming in February to April, the impact of clouds is much greater than air pollution and CO 2 , which is a reasonable explanation.
In the Amazon region of the southern hemisphere, stepwise multiple linear regression analysis between LST and precipitation, shortwave radiation, and longwave radiation did not show a high degree of significance at inter-annual scales. The strongest increase was found to occur during the dry season. The Amazon is an area with complex hydrometeorological conditions. Thus, inter-annual changes in the area may be determined by different driving factors in different years. Precipitation plays a major role, especially in dry years, followed by shortwave radiation in non-drought years. Extreme drought events over the Amazon basin have occurred frequently in this area in recent decades [85], and thermal anomalies can be detected by LST changes in the Amazon's drought event [86]. A recent increasing trend was also reported in the vapor pressure deficit (VPD) over tropical South America during dry months and reduced cloud cover, leading to increased incoming solar radiation in the southeast Amazon [87]. However, during the time period of this study, the decline in AOD more significantly affected the increase in shortwave radiation than the clouds in the Amazon region. AOD was also found to noticeably decline in its inter-annual and seasonal variability in Brazil during 2001-2014, which may be related to the decrease in deforestation fires since 2004 [88]. Although this study depends solely on statistical methods and analyses, which have known limitations, we are confident that the results obtained can serve as a starting point to further investigate the underlying mechanisms of these variabilities in long-term LST.

Conclusions
This study investigated the trends in annual average LST during 2003-2017 over each 1 • grid cell at the global scale using two satellite products (AIRS and MODIS) and one reanalysis product (ERA5-Land). All three products showed an increasing global trend in the annual average LST during 2003-2017, with the strongest increases found over the region north of 45 • N, particularly over Asian Russia, with a slight decrease over Australia. The spatial distribution of LST trends from these three products agrees very well; the correlation coefficient is 0.92 between AIRS and MODIS, 0.74 between AIRS and ERA5-Land, and 0.73 between MODIS and ERA5-Land. The major reason for the differences between the satellite observations and reanalysis LST trends appears to be the consideration of inter-annual vegetation dynamics.
Although the reanalysis product did not consider the inter-annual vegetation dynamics, the high correlation coefficient with satellite-based LST trends suggests the importance of climate factors in driving LST changes. Regression analysis showed that the major climate factors for LST changes during 2003-2017 varied over different regions. The interannual LST variations over the region north of 45 • N were mainly influenced by the incoming surface longwave radiation (LW↓), while precipitation (P) and incoming surface solar radiation (SW↓) played a more important role over other vegetated regions, except over tropical forests, where the inter-annual LST variation was observed to be very low, so no dominant climate factor was identified at an inter-annual scale.
Further analysis over Asian Russian and the southern Amazon at a monthly scale showed that the major climate factors for LST changes also vary temporally. Over Asian Russia, LST increases were predominantly observed during the winter to spring months (February-April), which was closely related with a simultaneous increase in LW↓ and clouds. Over the southern Amazon, the strongest LST increase was found in the dry season (August-September) and was primarily affected by a decrease in P. The increase in SW↓ associated with decreasing atmospheric aerosols was another contributing factor. The results of this study demonstrate that using both satellite observations and reanalysis can more comprehensively characterize the spatio-temporal dynamics of LST and allows us to better understand the major driving factors of LST changes over different regions.