A climatology of surface–air temperature difference over the Tibetan Plateau: Results from multi‐source reanalyses

The Tibetan Plateau (TP), known as earth's “Third Pole,” influences regional and even global weather and climate systems through its mechanical and thermal‐dynamical forcing. Near‐surface (2 m) air temperature (Ta) and surface (skin) temperature (Ts) are two crucial parameters of land–atmosphere interactions and climate change. Their difference (ΔT = Ts − Ta) determines the heating source over the TP that drives the Asian summer monsoon. This study focuses on climatology, inter‐annual variability, and long‐term trend of ΔT over the TP in the last four decades (1979–2018), based on four latest reanalysis datasets including ERA‐Interim, ERA5, MERRA2, and JRA55, along with observational data. We show that ΔT‐based different datasets display fairly different climatology in terms of seasonality, spatial distribution, and long‐term trend. ΔT exhibits a clear seasonality with negative value in winter and positive ones in summer despite different strengths and timings presented by the reanalyses. Along with global warming, all reanalyses except JRA55 exhibit obvious downwards trends of ΔT in a spatially non‐uniform way. The median ΔT among the four reanalyses features uniform decreases in all seasons, with the most distinct area on the northern TP, as well as the largest and least decreases in autumn and spring, respectively. Further analysis shows that the differences in ΔT are most likely associated with discrepancies in radiation forcing, snow cover, wind speed, and boundary layer height within the reanalyses. The present findings highlight the difficulty for the state‐of‐the‐art reanalyses to represent the climate change over the TP and point to possible factors behind the deficiencies identified.

atmosphere interactions . Previous studies have attempted to employ ΔT to conduct remote sensing retrieval of energy fluxes on the TP  and to assess the capability of climate model (Koven et al., 2013). Compared with other regions at same Northern Hemisphere latitudinal zone, the TP acts as a huge heat source to the atmosphere (Duan and Wu, 2008;Yang et al., 2011), mainly through sensible heat flux (SHF). Driven by the heating, normally termed as sensible heat driven air-pump , the water vapour and air mass in the low-level atmosphere around the plateau are pumped to the TP during spring and summer. ΔT is also a direct and main contributor to SHF which plays a crucial role in the seasonal transition, onset, and maintenance of the Asian summer monsoon, affecting precipitation in China . Therefore, deep understanding of the spatial patterns and variations of ΔT over the TP is of great importance for elucidating their influence on SHF and heat source.
The entire TP has poor coverage of meteorological stations, most of which are concentrated on the central and eastern TP and low-altitude areas . This makes it difficult and speculative to comprehend key processes crucial to regional climate and land-atmosphere interactions . Meanwhile, for routine T s meteorological observation, China Meteorological Administration widely adopted automatic observation to replace previous manual observation after 2005. As the difference in criterion exists between the two methods, T s recorded by automatic stations is generally higher in winter when snow cover exists than by manual observation (Liao et al., 2019). T s is measured at the 0-cm level, mainly reflecting the bare land skin temperature, but the central and eastern TP is actually covered by vegetation (Zhu et al., 2012). Thus, there remains a large uncertainty in T s from meteorological stations. And the observed ΔT do not reflect negative value and the wellknown cold source over the TP during winter and may overestimate the summer heat source (Liao et al., 2019). This severely hampers our understanding of the regional ΔT change. These factors all motivate us to adopt alternative dataset which can provide much better geographic coverage and supporting climate variables.
The reanalysis products become increasingly acceptable and have been in vogue to be used by a greater number of scientists since it can best describe long-term and large-scale thermo-dynamical state of the atmosphere (Gao et al., 2019). By combining in-site observations and remote sensing data into an atmospheric model via data assimilation, model processes are constrained as much as possible to the real atmospheric state. This improved reliability enables attractive applications for reanalysis data in a series of scientific studies (Hinkelman, 2019). The prevalent reanalysis data with high-resolution performs well in reproducing the TP climate variables. For instance, MERRA has a high correlation with observations of surface meteorological variables (Wang and Zeng, 2012). ERA-Interim data shows relatively small cold bias in T a (Wang et al., 2017) and better performance in water cycle (Wang and Zeng, 2012), as well as SHF on the TP .
Compared to the individual variable, either T a or T s , the spatial-temporal patterns of ΔT are more complicated (Feng and Zou, 2019). In recent years, numerous studies suggested that the SHF over the TP has undergone a weakening trend (e.g., Duan and Wu, 2008;Wang et al., 2013;Duan et al., 2018), which is attributed to the decreased surface wind speed (Duan and Wu, 2008). However, there has been relatively little systematic research on another important contributor ΔT over the TP, in particular, including its seasonality, inter-annual variability, and long-term trend. The objectives of this study are to identify these based on latest multi-source reanalysis datasets, and to further investigate possible causes for the discrepancies in ΔT among the reanalyses.
2 | DATA AND EVALUATION STRATEGY

| Reanalysis datasets
Four reanalysis products containing T s and 2-m T a (hereafter referred to simply as T a ) are comprehensively used in the study (Table 1), since they are the latest generation and state-of the-art reanalyses. Monthly T a and T s are obtained from the Europe Centre for Medium-Range Weather Forecasting (ECMWF) ERA-Interim (hereafter ERAI) with a spatial resolution of~79 km or~0.75 for the period 1979-2018. The ERAI is an updated version of the ERA40 reanalysis and utilizes the data assimilation system of the Integrated Forecast System (Cy31r2). Besides, ERAI employs a four-dimensional variational analysis (4D-Var) and has improved variational bias correction of satellite-derived radiance data as well as humidity analyses and model physics (Dee et al., 2011).
Recently, ECMWF has developed its fifth-generation reanalysis product, namely ERA5 (the successor to ERAI) with a spatial resolution of~31 km or~0.25 (Hersbach et al., 2018). Relative to the former ERAI, substantial changes have been made for ERA5 in terms of higher temporal and spatial resolutions. ERA5 integrates ample amounts of reprocessed observations (satellite, ozone, aircraft, and surface pressure data) into global estimates using advanced earth system model and data assimilation system of the IFS (Cy41r2).
In a reanalysis, a skin layer represents the vegetation layer, the top layer of the bare soil, or the top layer of the snow pack in its land surface model. T s indicates the temperature at the interface between the land and atmosphere. In the ECMWF reanalysis family T s is physically acquired from the surface-atmosphere energy balance of the land surface model, which connects the surface with the lowest atmospheric level through dry static energy, moisture, and thermal contact with a single four-layer soil profile (or one layer if snow exists; Dee et al., 2011). Another reanalysis is the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA2), which is released by the NASA Global Modeling and Assimilation Office. MERRA2 has a native spatial resolution of 0.5 × 0.625 and is available from 1980 onwards. Modern hyperspectral radiance and microwave observations, as well as GPS-Radio Occultation datasets, are incorporated into the NASA's Goddard Earth Observing System version 5 (GEOS-5) via the GSI assimilation system (Gelaro et al., 2017). Thanks to these advances, many aspects in MERRA2 have been updated, including the representation of the cryosphere, T a , and global mean radiative fluxes (Hinkelman, 2019). MEERA2 T s is from the catchment land surface model, a component of the GEOS-5, and assimilation system (Draper et al., 2018).
In addition, we also make use of the Japanese 55-year Reanalysis (JRA55; 1.25 × 1.25 grid) conducted by the Japan Meteorological Agency (Kobayashi et al., 2015). T s is represented by the brightness temperature in JRA55, which is computed from surface upwards longwave radiation assuming that surface acts as a black body. Similar to ERAI, JRA55 incorporates snow observations over the TP region.
To investigate the differences in ΔT variations among reanalyses, downward shortwave and longwave radiations (DSR and DLR), snow depth, wind speed at 10 m above ground level (WS), and boundary layer height (BLH) are extracted from the corresponding reanalysis datasets. Following Orsolini et al. (2019), snow depth is converted to snow cover fraction (SCF).

| Observational datasets
Two sets of ancillary observational data of 2-m T a are also utilized. Observational CN05.1 gridded data (hereafter CN05, a spatial resolution of 0.25 × 0.25 ) covering the TP is used, which is interpolated based on 2,416 stations in China (Wu and Gao, 2013). It has been widely used in the validation and evaluation of climate simulations and climate change (Wang et al., 2016;Guo et al., 2018). For comparison purpose, another monthly T a data is obtained from the gridded observations CRU (Climate Research Union) with a spatial resolution of 0.5 × 0.5 produced by University of East Anglia (Mitchell and Jones, 2005). Likewise, the gridded temperatures are interpolated from quality-controlled meteorological station data based on a presumed correlation decay distance. Furthermore, we also resort to the monthly satellite land surface temperature (LST) datathat is MODIS/Terra MOD11C3 Collection 6 covering the time period from February 2000 to December 2018 at 0.05 spatial resolution (Wan et al., 2015)-as a reference data to manifest the performance of the reanalysis products in T s .

| Evaluation strategy
Since the time spans of these reanalysis datasets are somewhat different, this paper simply focuses on the common time period 1979-2018 (40 years) except for MERRA2 starting from 1980 and for CN05 and CRU being used until 2016 and 2017, respectively. Due to different spatial resolutions (Table 1), for the spatial distribution and trend of the median ΔT among the four reanalyses, we interpolate T s and T a into a common resolution of 1.5 × 1.5 based on the bilinear interpolation method, although this may obscure some of the added value of high-resolution data and introduce some errors because of complex climate and underlying surfaces.
The performances of the four reanalysis datasets in T a and T s are evaluated prior to the ΔT analysing. To examine the statistical difference of the two variables, significance test for the difference ΔT is implemented using the twotailed Student's t test (Wilks, 2006). Linear regression trends based on the least-squares method are also applied for regional averaged winter (December-February, DJF), spring (March-May, MAM), summer (June-August, JJA), and autumn (September-November, SON) climate variables. The significance test for the linear trend is performed using the two-tailed Student's t test (Santer et al., 2000).
The median of the multi-reanalysis products is computed, which can reduce the influence by extremely large and small values . Relationships (r) between climate variables (DSR, DLR, SCF, WS, and BLH) and ΔT within the same reanalysis product are detected using the Pearson correlation analysis. We also use the seasonality and long-term trends of climate variables, together with the correlations to explain the differences in ΔT among the reanalyses.

| Evaluation of temperature (T a and T s )
We first compare climatological annual cycle of T a from observations (CRU and CN05) and reanalyses ( Figure 2).
All reanalyses reproduce reasonably well the annual changes of T a . T a in ERAI is slightly higher than the observations, consistent with the result of Wang et al. (2017). Yet other three reanalyses have a negatively bias partly due to the fact that most of meteorological stations record warm temperature since they are concentrated in inhabited valleys.
Intercomparisons among the reanalyses show that ERA5 produces relatively cold temperatures, particularly in winter (Figure 2a and Table S1, Supporting Information). Topography may play a vital influence on nearsurface temperature. There are evident elevation differences among these reanalyses (Figure 1). A noted fact is that much finer resolution model can better depict topography and likely generates high terrain, leading to a cold temperature. Nevertheless, coarse resolution modelling generates relatively low and flat terrain, resulting in a warm temperature. Between the ECMWF family reanalyses, ERAI has a better performance in reproducing snow depth over the TP (Orsolini et al., 2019). It is due to the fact that snow cover from Interactive Multisensor Snow and Ice Mapping System (IMS) was considered by ERAI, even in high-altitude regions, including the TP, whereas IMS data above 1,500 m was not assimilated in the production of ERA5, resulting in overestimated SCF and days on the TP, nearly 3 times of observed data (Orsolini et al., 2019). Consequently, ERA5 most likely simulates much stronger albedo effect of snow, forming cold temperature. This feature is more obvious in cold season (Table S1). Compared with CN05, MERRA2 and JRA55 have small biases and RMSE, ERAI performs best in T a except for DJF series.
When it comes to T s , though appreciable differences exist, the annual change is consistently reproduced by the four reanalyses. Warmest T s is produced by MERRA2 across the whole year (Figure 2b), coldest by JRA55 during summer, ERA5 still tends to yield cold T s with exception of summer when it shares nearly similar T s with ERAI (Table S1). These differences may be attributable to numerous factors, such as diverse land surface schemes, atmospheric models, and observational data assimilating sources (Table 1). Land surface schemes which emphasize different land surface processes, such as vegetation, snow melting process, and soil water-heat transfer, simulate diversified partition of surface energy (Yang et al., 2009;Wang et al., 2014;. Additionally, atmospheric models' resolution and physics design, together with the use of data assimilation system, have large impacts on atmospheric circulations and surface energy and water budgets. In comparison with MODIS LST, all reanalysis data generate warmer T s in JJA but colder ones in the rest of the seasons. JRA55 forms much closer T s with MODIS data in DJF and JJA, while MERRA2 shows the best performance in MAM (Table S1). Considering above comparisons, no reanalysis clearly stands out as the best performing one for all seasons and variables, and each of the reanalyses displays different biases. It is highlighted that combined utilization of reanalysis data with various data sources is beneficial for understanding the TP ΔT changes. Figure 2c shows the climatological annual cycle of ΔT over the TP. Overall, highest ΔT is attained by MERRA2, lowest ΔT by ERA5 during winter (reaching up to −3.4 C) and spring and by ERAI during transitional seasons. Although there are some differences in the magnitudes of ΔT among the four reanalyses, the variations are basically consistent. However, the timing of the peak has a disagreement across the reanalyses. ΔT in MERRA2 and JRA55 peak in June, 1 month earlier than the ECMWF family of reanalyses. These may reflect inconsistencies in the activity and onset of summer monsoon among these reanalyses. The median of the four reanalyses shows that ΔT has strong seasonal variability, with maximum in JJA (1.7 C) and minimum in DJF (−1.8 C).

| Spatial-temporal distribution of ΔT
In terms of the spatial distribution (Figure 3), CN05 and all the reanalyses indicate that the TP has a largescale negative T a and ΔT in DJF, suggesting an obvious cold source over the TP in winter. From the ECMWF datasets, substantial negative ΔT occur over high mountains, such as the Karakorum Mountains, the Himalayas, the mountains in the southeastern TP, and the Qilian Mountains. More details are generated by ERA5 because of its higher resolution. Unlike the ECMWF reanalyses, MERRA2 produced large ΔT mainly appear along the Himalayas and in the Qilian Mountains. In spite of sharing approximately the similar spatial resolution, the result of ΔT between ERAI and MERRA2 is quite diverse. JRA55 shows large ΔT appearing on the northwestern TP. The median exhibits the strongest cold source (Figure 3f), as indicated by the most negative ΔT.
In MAM, T a gradually increases, the spatial patterns of ΔT manifested by the four reanalyses are quite different. Positive ΔT in the ECMWF begin to appear over the western TP and the Qaidam basin, while negative ones still exist in the southeastern TP and the Himalayas with decreased strengths. MERRA2 shows entirely positive ΔT on the TP. As for JRA55, positive ΔT is located on most of the TP. The median resembles the spatial pattern of ECMWF reanalyses (Figure 3l), but the regional mean ΔT quickly turns into positive value (Table S2), signifying a shift to a heat source for the TP.
As T a turns entirely positive in JJA, the four reanalyses reveal similar pattern with positive ΔT over the vast majority of the TP. Large ΔT are located in the western and northern TP, and some part in the southern TP. MERRA2 has totally positive ΔT and the largest value, while ERAI and JRA55 possess the same and smallest ΔT (Table S2). The median result indicates higher ΔT appears in the western TP than in the eastern TP, yet smallest ΔT in the southeastern TP where monsoon precipitation occurs frequently, reducing T s and ΔT ( Figure 3r). The median ΔT is the highest among four seasons, denoting the heat source reaches its greatest strength in JJA.
As solar radiation weakens (T a decreases accordingly) and the monsoon withdraws, ΔT decreases rapidly from JJA to SON and then quickly transforms into negative values, with a regional averaged value of −0.3 C. It means that the heat source reduces its strength. Discrepancies of inter-reanalysis are appreciable, which is comparable to those in MAM. The median displays negative ΔT in the southeastern TP and positive value in the northwestern TP (Figure 3x). Figure 4 shows the spatial distribution of seasonal trends in T a and ΔT over the TP. CN05 reveals significant increasing trends in T a in all four seasons, with highest in DJF and smallest in JJA. ΔT change varies spatially among these reanalyses. The two ECMWF reanalyses show a large area of decreasing trend in ΔT, most notably in the northeastern TP, while somewhat increasing trend appears in the western TP, the Himalayas, and the southeastern TP. Comparatively high decreasing trend is represented by ERAI, especially in DJF and SON. With exception of DJF, MERRA2 displays a spatial pattern marked by an increasing trend in the western TP but a decreasing trend in the eastern TP. Areas with negative and positive trends account for approximately half of the TP, respectively. On the contrary, JRA55 presents distinct increasing trends in ΔT particularly during DJF apart from the northern TP. Overall, ΔT in the median presents widespread downwards trends across the four seasons with the most prominent area occurring in the northern TP, while somewhat upwards trends are observed in the western and southeastern TP. The time series of seasonal average ΔT for the whole TP from the reanalysis products exhibit marked discrepancies ( Figure 5). The ECMWF reanalyses, especially ERAI, present significantly decreasing trends in ΔT.

| Inter-annual variation of ΔT
MERRA2 gives almost invariable trends in DJF and MAM but slightly decreasing trends in the remaining seasons. As to JRA55, there are significantly increasing trends, especially in DJF. Moreover, ΔT time series of the reanalyses exhibit large fluctuations in DJF and MAM, indicating large year-to-year variability. Similar decreasing trends are also reported by NCEP reanalysis data . Taken together, the median ΔT has decreased significantly in four seasons from 1979 to 2018 (Table S2), in accord with the decline of SHF (Wang et al., 2013;Duan et al., 2018) and the reduction of heating ability in recent decades (Yang et al., 2011). The largest decrease occurs in SON at a rate of 0.045 CÁdecade −1 , while the least decrease appears in MAM at a rate of 0.037 CÁdecade −1 . Several complex interactions for the TP rapid warming have been detected in recent decades. For instance, the TP warming is due to CO 2 increase (Chen et al., 2003). Locally, radiation balance caused by snow albedo feedback imposes great impacts on T a and T s (e.g., Chen et al., 2017;Gao et al., 2019;Guo et al., 2020). WS decline lessens the Bowen ratio and generates less surface sensible heating, accounting for the TP warming . In this section, downward radiation forcing, SCF, WS, and BLH are used here to investigate the causes of the differences in ΔT among the reanalyses.

| Downward radiation forcing
Radiation forcings (DSR and DLR) are expected to have dominant influences on ΔT as they are the control factors of temperature. Figure 6a shows scatterplot of downward radiations versus ΔT. Both DSR and DLR have significantly positive correlations with ΔT. This indicates that strong downward radiation forcing tend to generate a positive ΔT, vice versa. Furthermore, DLR has a larger impact on ΔT than DSR in terms of higher correlation coefficient except for MERRA2, mainly because there is similar seasonal evolution between DLR and ΔT, that is, highest in JJA and lowest in DJF (Figure 6d). MERRA2 has stronger DSR and reaches its peak in June, which can impose great influence on ΔT, such as large positive ΔT and advanced peak time. Over time, DSR has a significant decreasing trend and DLR has an increasing trend in JJA. The increases of DLR in the reanalyses except for JRA55 are weaker than the decreases of DSR (Table 2), thus causing the decreases in ΔT. However, it is exactly the opposite for JRA55. In DJF and MAM, DSR and DLR display slightly positive trends over the past 40 years but exert complicated influences on ΔT because of their intensities and snow cover.

| Snow cover
A prior study suggested that the albedo effect is governed by SCF (Xu and Dirmeyer, 2012). To examine the effect  of snow cover on ΔT, we plot the correlation map between SCF and ΔT over the TP on monthly time scale (Figure 7a). There are negative correlations between monthly SCF and ΔT. It means large (small) SCF appears to yield negative (positive) ΔT. Climatologically, the large ΔT generally occur at high altitudes where large SCF appear (Figure 7b). During cold season, high snow cover tends to reflect solar radiation, leading to a lower T s and a more negative ΔT, and vice versa. Among the four reanalyses, ERA5 generates the highest SCF over the TP, followed by JRA55, ERAI, and then MERRA2, in accord with the result of Orsolini et al. (2019). That helps to explain ERA5 has large negative ΔT in DJF and MAM. In the case of MERRA2, the smallest SCF is more likely to lead to less negative ΔT, and the rapid response of ΔT to SCF is partly the cause for ΔT's dramatic transition from negative in DJF to positive in MAM. For JRA55, the response of ΔT to SCF is slow from snow season to snow free season (JJA). Snow cover has been reported to have reduced over the past decades , characterized by greater decreasing magnitudes in DJF and MAM than in JJA (Table 2). High-elevation warming is ascribed to snow cove/ice-albedo feedback. Accordingly, some increases of ΔT in high-elevation areas revealed by ECMWF reanalyses are likely associated with the decrease of SCF in DJF and MAM. For JRA55, notably reduced SCF in DJF together with increased radiation forcing act to cause a large increase in ΔT. Large year-toyear variability of the snow cover among the reanalyses may contribute to the great diversity of ΔT.

| Wind speed
Wind speed (WS) is a dynamic factor affecting the exchange of turbulent heat fluxes. Since WS can affect land surface energy balance through evapotranspiration and surface heat flux, WS change is constrained by processes that have a scale beyond the atmospheric boundary layer . WS also exhibits a strengthening-weakening seasonal variation from DJF to SON, with a maximum value in MAM (Figure 8b), meaning that the peak time of WS is ahead of ΔT. Among the reanalyses, MERRA2 has the largest WS, while JRA55 has the smallest WS. At the monthly scale, negative correlations suggest that high WS is likely to produce small ΔT (Figure 8a). In terms of WS change, the reanalyses reach a consistent result that shows increases in JJA and SON, which are in agreement with the decreases of ΔT. Between the ECMWF reanalyses, having nearly similar variation tendencies of DSR and DLR in SON, the decreases of ΔT are very different, partly due to the obvious difference of WS variation (Table 2).

| Boundary layer height
The development of planetary boundary layer is strongly relevant to convection activity and complex surface forcing containing SHF, frictional drag, evapotranspiration, and terrain-induced flow modification (Stull, 1988). The effective heat capacity, resolved by the air column with mixed heat, influences T a to some degree. This heat is determined by a surface heat flux and runs through the PBL. The shallower BLH, the lower effective heat capacity, and the more sensitive T a is to forcing (Esau et al., 2012). BLH is missing in JRA55, and thus the correlations are calculated using the other three reanalysis products (Figure 9a). Positive correlations between monthly BLH and ΔT imply large BLH results in positive ΔT. The largest ΔT is generated by MERRA2 in JJA and F I G U R E 7 As Figure 6, but for snow cover fraction (SCF) [Colour figure can be viewed at wileyonlinelibrary.com] the smallest by ERA5 in DJF (Figure 9b), corresponding to their respective largest BLH and smallest BLH. A strongly amplified temperature response can be obtained in shallow boundary layers (Davy and Esau, 2016). In a warming climate, obvious decreasing trends in BLH are found in ERAI in the four seasons, which amplifies T a warming response and causes the decreases of ΔT. In JJA, BLH have reduced significantly in all reanalyses, contributing to the ΔT decreases (Table 2).

| SUMMARY AND DISCUSSION
The study utilizes T s and T a data from four recent and state-of-the-art reanalysis datasets (ERAI, ERA5, MERRA2, and JRA55) to reveal spatiotemporal variability of ΔT over the TP from 1979 to 2018. The TP ΔT exhibits a strong seasonal cycle with negative value in DJF and positive value in JJA, suggesting the seasonal transition from heat sink to heat source. However, there are significant differences in terms of peak time and strength of ΔT among the reanalyses.
Under climate warming, apart from JRA55, the reanalyses present decreasing trends of ΔT in a spatially non-uniform way over the studied period. Generally, the median of the four reanalyses exhibits significantly negative trends across four seasons with the most pronounced area in the northern TP and the largest decrease in SON as well as the least decrease in MAM.
The decreases in ΔT during DJF and SON suggest the strength of the TP cold sources have enhanced and increased near-surface temperature inversion have suppressed air convection and mixing. These results of decreasing ΔT are in favour of the decline of SHF, implying that the TP heat sources in MAM and JJA have reduced their strength, which may exert significant influences on the monsoon onset and intensity.
The differences in the spatiotemporal variability of ΔT among the reanalyses are relevant with radiation forcing, SCF, WS, and BLH at different intensities and time scales (mainly seasonality and inter-annual variability). Overall, ΔT is positively correlated to DSR and DLR as well as BLH, but negatively correlated to SCF and F I G U R E 8 As Figure 6, but for 10-m wind speed (WS) [Colour figure can be viewed at wileyonlinelibrary.com] F I G U R E 9 As Figure 6, but for boundary layer height (BLH). BLH is not available in JRA55 [Colour figure can be viewed at wileyonlinelibrary.com] WS. Additionally, other causes such as precipitation, vegetation condition, and soil moisture may also play important roles. T s is response to wetting and drying cycles of soil, which are heavily modulated by precipitation. The vegetation is normally denser in wetter region. The vegetation types and their description in the reanalyses are remarkably different (Xie et al., 2019). ERA5 and MERRA2 adopt seasonally varying vegetation index to represent vegetation dynamics, but do not allow vegetation to change annually. Moreover, the existing reanalyses do not depict vegetation changes explicitly in the TP. These aspects should be kept in mind in future research.
Given high heterogeneity among reanalysis products covering the climatology of ΔT at the seasonal and interannual time scales, using the reanalyses to represent climate change is not yet mature and should be refined in the future. Meanwhile, the different usage of data in assimilation system such as snow data in the ECMWF reanalyses, and the diverse physical parameterizations, that is, land surface and boundary layer schemes, allow highlighting the importance of incorporating more satellite data and station observations into assimilation system and optimizing model parameterizations. Although interrelationships between influencing factors and ΔT are performed based on the identical reanalysis source to alleviate errors from different datasets, inherent biases are still inevitable. Therefore, a quality-controlled database containing various variables corresponding to T s and T a should be generated in future research to verify the results and further explore their difference and coupling. What's more, regional climate models-based research is required to enhance understanding of the regional details, as topographic features and their crucial effects on regional climate are not well resolved in coarse-resolution GCMs/reanalyses (Gao and Chen, 2017;Giorgi, 2019), as well as to estimate the impacts of dynamic vegetation and soil moisture.