The performance of CORDEX-EA-II simulations in simulating seasonal temperature and elevation-dependent warming over the Tibetan Plateau

To explore the driving mechanisms of elevation-dependent warming (EDW) over the Tibetan Plateau (TP), the output from a suite of numerical experiments with different cumulus parameterization schemes (CPs) under the Coordinated Regional Climate Downscaling Experiments-East Asia (CORDEX-EA-II) project is examined. Results show that all experiments can broadly capture the observed temperature distributions over the TP with consistent cold biases, and the spread in temperature simulations commonly increases with elevation with the maximum located around 4000–5000 m. Such disagreements among the temperature simulations could to a large extent be explained by their spreads in the surface albedo feedback (SAF). All the experiments reproduce the observed EDW below 5000 m in winter but fail to capture the observed EDW above 4500 m in spring. Further analysis suggests that the simulated EDW during winter is mainly caused by the SAF, and the clear-sky downward longwave radiation (LWclr) plays a secondary role in shaping EDW. The models’ inability in simulating EDW during spring is closely related to the SAF and the surface cloud radiative forcing (CRFs). Furthermore, the magnitude and structure of the simulated EDW are sensitive to the choice of CPs. Different CPs generate diverse snow cover fractions, which can modulate the simulated SAF and its effect on EDW. Also, the CPs show great influence on the LWclr via altering the low-level air temperature. Additionally, the mechanism for different temperature changes among the experiments varies with altitudes during summer and autumn, as the diverse temperature changes appear to be caused by the LWclr for the low altitudes while by the SAF for the middle-high altitudes.


Introductions
The Tibetan Plateau (TP), the largest and highest plateau in the world, exerts a profound impact on the global and regional climate (Duan and Wu 2005;Sato and Kimura 2007;Wu 1984), and is regarded as an early warning indicator of the global warming (Yao et al. 2000). In particular, the TP holds the largest ice outside the polar regions, providing freshwater resources for over 1.4 billion people downstream (Chen 2015;Immerzeel et al. 2010;Yao et al. 2012). Therefore, a deeper understanding of climate processes taking place over the TP has important implications for the hydrological regime response, and ecosystem services.
Numerous studies based on in-situ observations confirm significant warming over the TP during recent decades (Duan and Xiao 2015;Liu and Chen 2000;Rangwala et al. 2009), and most of them demonstrate that the warming rate 1 3 is amplified with elevation (Duan and Xiao 2015;New et al. 2002;Qin et al. 2009), which is referred to as elevationdependent warming (EDW). However, the available meteorological stations over the TP are unevenly distributed and mainly located over relatively low altitudes (You et al. 2020), which limits the availability of climate information for the TP and our understanding of relevant processes causing the EDW, especially for the high elevation regions.
Climate models, which overcome the inadequacies of observations and provide enough variables to unveil the physical mechanisms driving EDW, have been widely used in temperature simulations and projections over the TP (Gao et al. 2015;Ma et al. 2019;Palazzi et al. 2017Palazzi et al. , 2019Rangwala et al. 2010;Su et al. 2013;Wang et al. 2016). However, global climate models (GCMs) with coarse resolution do not adequately resolve the terrain feature and local climate, limiting their ability in simulating the climate over the complex terrain (Giorgi and Bates 1989;Salathé et al. 2008;Su et al. 2013). Regional climate models (RCMs) on the other hand, with improved representation of fine-scale topography and reasonable description of interaction among regional earth system components (Fu et al. 2005), show a clear advantage over the GCMs in simulating temperature over the TP Zhang et al. 2019), though large cold biases exist (Gou et al. 2019;Meng et al. 2018;Wang et al. 2013).
So far, only a few RCM studies have focused on the characteristics of EDW. For example, based on the results from the Weather Research and Forecasting Model (WRF) at 30 km resolution, Gao et al (2015Gao et al ( , 2018 show that a clear EDW signal occurs below 5000 m whereas no EDW exist above 5000 m over the TP, and this phenomenon will continue to exist in the future. Similar conclusions were drawn by Guo et al. (2016) based on two RCMs projections that showed the EDW over the TP does not expand to 4800 m in the future. They also highlight that the amplitude and location of the maximum warming rate depend on the intensity of warming in the future. Additionally, a recent study reveals that the relationship between temperature change and elevation varies widely across the RCMs over the Himalayan region (Nengker et al. 2018).
EDW is the response to various processes that have feedbacks at different scales, and a variety of mechanisms have been proposed to explain it (Rangwala and Miller 2012;Rangwala et al. 2010;You et al. 2020). Surface albedo feedback (SAF) is considered as the main driver for the EDW, as evidenced by the maximum warming occurring near the annual 0℃ isotherm and by strong correlations between simulated patterns of warming and albedo reductions (Giorgi et al. 1997;Minder et al. 2018;Pepin and Lundquist 2008;Yan et al. 2016). The SAF is more responsible for the change in the daily maximum temperature (Tmax) than the daily minimum temperature (Tmin) . Another important contributor to the EDW is the downward longwave radiation, especially for Tmin during winter (Palazzi et al. 2017(Palazzi et al. , 2019Rangwala et al. 2009Rangwala et al. , 2013. Furthermore, other mechanisms, such as cloud feedback (Duan and Wu 2006;Yan et al. 2016), aerosols radiative forcing (Rangwala et al. 2010), and free tropospheric warming (Bradley et al. 2004;Rupp et al. 2017), are also used to explain the EDW. However, it remains difficult to untangle the causal correlations among the relevant variables and to quantify the main process leading to the EDW (Pepin et al. 2015). Besides, the characteristics of the simulated EDW and the relative role of driving factors may depend on the adopted RCM configurations (Minder et al. 2018). To assess the impact of model configurations on the simulated EDW and clarify the mechanisms driving EDW over the TP, further investigations based on the high-resolution simulations are still needed.
In an RCM, the cumulus parameterization scheme (CP) and land surface model (LSM) are considered to play crucial roles in simulating temperature . Many sensitivity analyses have been conducted to investigate the sensitivity of EDW to LSM (Minder et al. 2016(Minder et al. , 2018. While for CPs, most of the studies mainly focus on the precipitation over the TP (Ali et al. 2015;Gu et al. 2020;Ou et al. 2020), and the sensitivity of simulated EDW to CPs is largely unknown. CPs strongly affect the surface radiation balance and temperature via altering the timing and partitioning of precipitation between convective and non-convective precipitation (Lynn et al. 2004). Besides, CPs determine the amount of precipitation and some of them occur in the form of snow over the TP, which influences the surface radiation budget and ultimately near-surface air temperature simulation. All of these imply that different CPs may be an important source of uncertainty in EDW and their effects need to be systematically evaluated. Here, five sensitivity experiments have been conducted by the WRF model (Skamarock et al. 2008) with different CPs under the framework of Coordinated Regional Climate Downscaling Experiment-East Asia (CORDEX-EA-II). In this study, we focus on the comparations of WRF performances in simulating the temperature and EDW over TP, intending to investigate the major contributing factors to EDW, and to assess the sensitivity of simulated EDW to the CPs in WRF.
This paper is structured as follows. Section 2 briefly describes the model, datasets, and methodology used in the analysis. Section 3 presents the simulations of temperature climatology and EDW in the five experiments and diagnoses the main process responsible for inter-model differences and the simulated EDW. The sensitivity of EDW to the CP in WRF is also examined. Finally, Sect. 5 summarizes the main findings and conclusions.

Model and sensitivity experiments
To assess the impact of CPs on the temperature simulation and EDW over the TP, five CPs, including the Betts-Miller-Janjic (BM, Betts and Miller 1986;Janjić 1994), the Kain-Fritsch (KF, Kain and Fritsch 1990), the KF scheme with an alternative trigger function weighted by moisture advection (KFMT, Ma and Tan 2009), the Grell-3D ensemble (GR3D5, Grell and Dévényi 2002), the updated version of the Simplified Arakawa-Schubert (SAS, Han and Pan 2011), are selected for sensitivity experiments in which the other components are fixed. All the experiments are driven by ERA-Interim reanalysis (Dee et al. 2011(Dee et al. ) for 22 years (1988(Dee et al. -2009, and adopt the same domain with a quasi-regular grid spacing of about 25 km covering the CORDEX-EA-II region (Fig. 1). The configurations of these simulations are summarized in Table 1, and a detailed description of the model setup can be found in Niu et al. (2020).

Observational gridded datasets for the air temperature over the TP
Three variables including daily mean temperature (T2m), Tmax, and Tmin from the gridded observational dataset of CN05 (Wu and Gao 2013) are used to evaluate the performance of RCMs over the TP. This gridded dataset with 0.25° × 0.25° resolution is constructed using the anomaly interpolation method (New et al. 2002) from 2416 weather monitoring stations across China. The climatology is first interpolated by thin-plate smoothing splines, in which the fitted surface are functions of latitude, longitude, and elevation. Also, the elevation is a co-predictor and a topographic correction is applied during the interpolation. Then, a gridded daily anomaly using the angular distance weighting method is added to the aforementioned climatology to obtain the final dataset. CN05 has been widely used to investigate regional climate change information and validate high-resolution climate models (Wu et al. 2017;Shi et al. 2018). Note that the sparse in-situ stations over the western TP ( Fig. 1b) results in large uncertainty in the gridded data there and probably exaggerates model bias (Lucas-Picher et al. 2011).
For the convenience of comparison, both the observations and model results are interpolated into 0.25° × 0.25° resolutions over the TP using bilinear interpolation. To reduce the influence of elevation difference on temperature, a lapse rate (0.65 °C/100 m) correction is applied to the simulated temperature with help of the CN05 and model elevations when evaluating the models' performances.

Methodology
In order to quantify the main factors contributing to the inter-model differences and EDW over the TP, the changes in surface energy budget are diagnosed following the methods of Lu and Cai (2009). Based on the land surface energy balance equation and Stefan-Boltzmann law, the land surface temperature change can be computed as follows: where σ and Ts represent the Stefan-Boltzmann constant (5.67 × 10 − 8 Wm −2 K −4 ) and land surface temperature; LW ↑ , and LW ↓ , and SW ↓ refer to the surface upward and downward longwave and downward shortwave, respectively; Q is the ground heat storage; H and LE denote the sensible and latent heat flux. Additionally, the overbar stands for the unperturbed mean state and Δ means the difference between the perturbed and unperturbed mean state.
The surface cloud radiative forcing (CRFs), calculated by the difference between the net total-sky radiation and clear-sky radiation at the surface, could stand for the effects of cloud on climate (Cess and Potter 1988). Since changes in the CRFs include parts of surface albedo feedback (SAF) and other feedbacks, Lu and Cai (2009) exclude the SAF from the changes in CRFs as follows: where (•) cld = (•) − (•) clr , in which (•) and (•) clr indicate radiations for the total sky and clear sky condition, respectively; − means the surface albedo of the unperturbed mean climate state. As the radiation emitted from the surface is not directly related to the cloud, ΔLW ↑cld in Eq. (2) is close to zero. Therefore, using Eq. (2), the Eq. (1) can be rewritten as: The six terms on the right-hand side of Eq. (3), which represent the SAF, the difference in surface CRFs, the non-SAF-induced change in clear-sky downward shortwave fluxes (SW clr ), the change in downward clear-sky LW radiation (LW clr ), surface heat storage (Q), as well as surface sensible/latent fluxes (turbulent heat flux, THF), are analyzed to detect the main factors related to the simulated temperature spread and EDW. A detailed description of this equation can be found in Lu and Cai (2009). Note that we only calculate the energy components of Eq. (3) over the land surface.

The climatological temperature
Similar distributions are observed for T2m, Tmax, and Tmin in the observational dataset, and the temperature distribution follows the topography of the TP, with the highest temperature in the Qaidam Basin and southeast TP and lowest temperature over the northwest corner above 5000 m (Fig. 2), which reflects the dominant role of orography on the temperature distribution over the TP. All the experiments reasonably capture the spatial pattern of the temperature over the TP with spatial correlation coefficients above 0.75 for all seasons, and most of them show better capacity in simulating temperature distribution in autumn than other . Generally, cold bias is observed in all the experiments with its magnitude varying spatially, and the largest bias reaching up to − 10 °C occurs over the southwestern TP during winter and over the northwestern TP during other seasons ( Fig. 3 and Supplementary Fig. 1). On the other hand, the experiments depict a slight warm bias over the scattered areas of eastern TP. The Taylor diagram also indicates that there exist large inter-model differences among the temperature simulations over the TP during summer. The spreads in the simulated temperature over the TP, represented by the standard deviation of five simulations, are present in Fig. 5. The disagreements among the simulations show a strong elevation dependence in summer and autumn, evidenced by the more pronounced inter-model difference over the western TP than the eastern TP. While the largest spread in spring and winter occurs along the southern edge of TP and southeast part of TP. In addition to the spatial variation, the model spread due to CPs also shows a seasonal dependence, with the largest magnitude in summer and lowest in winter. Moreover, the simulated Tmin exhibits a larger inter-model spread compared to Tmax, indicating stronger sensitivity of simulated Tmin to the choice of CP over the TP. The relationship between the spread and elevation in the simulations is also investigated (Fig. 6). The spreads in the simulations mostly increase with elevation and reach their maximum around 4000-5000 m, demonstrating that the simulated temperature at high elevation is more sensitive to the choice of CP. Furthermore, the relationship between spread and elevation exhibits a seasonal dependence for T2m and Tmin over the high altitudes, as the discrepancy still increases with elevation above 5000 m during summer and autumn, while decreases in other seasons.

The structure and seasonal dependence of EDW
Accelerated warming is found over the TP during the early-2000s global warming hiatus (Duan and Xiao 2015;Ma et al. 2017). To investigate the structure of temperature change with elevations, we calculate the difference in temperature between two periods 2000-2009 and 1990-1999 as the temperature change and plot it with the elevation divided into 500 m altitudinal bands. As is shown in Fig. 7, the relationship between temperature change and elevation varies across seasons. In winter, the warming rate generally increases with elevation up to 5000 m and rapidly decreases with elevation further, suggesting an EDW signal below 5000 m over the TP in winter. This is also confirmed by the satellite observational results (Qin et al. 2009). While for other seasons, the EDW signal can only be observed at certain altitudes, such as higher altitudes (above 4000 m) during spring and autumn. Moreover, the profiles of EDW in Tmax and Tmin seem to resemble that in T2m, except for Tmin in summer which presents a clear EDW signal (Supplementary Figs. 2 and 3). Compared to the warm season (summer and autumn), the cold season (winter and spring) exhibits higher spatial variability of T2m change for the entire profiles, with its magnitude increasing with elevations except for regions below 2500 m. These diverse responses in the grid points within the same elevation bins may be caused by the local factors, such as slope, land cover, exposure, etc. (Pepin et al. 2015), and imply that the EDW profiles over the TP depend on sub-regions, which is evidenced by the contrasting EDW patterns in three mountain ranges over the TP .
Most of the experiments can reasonably reproduce the seasonal elevation profiles of temperature change below 4500 m over the TP with reduced magnitudes. However, substantial differences in the profiles between the observations and simulations occur over the higher altitudes except in winter. All the simulations fail to capture the observed EDW above 4500 m in spring, and they generate an EDW in summer for the elevations around 4500-5500 m where little elevation dependency is found in the observations. Such diversity is also documented by Gao et al. (2018), who showed that gridded observations generate EDW above 5000 m while neither RCM nor ERA-Interim suggests EDW over there. These results imply that whether the EDW exists over the high altitudes in the TP is still under debate and there is an urgent need to extend observations upwards to the high altitudes.
The structure and magnitude of simulated EDW are sensitive to the choice of CPs, especially during summer and autumn. For instance, KF simulates enhanced warming below 4500 m during autumn, yielding EDW variations reaching up to 0.5 °C, while other CPs produce small EDW variations with fairly uniform temperature change. The largest difference among the simulations is located over the high altitude (around 5500 m) during summer, reaching up to 1.3 °C, and over around 4500 m during other seasons with the value exceeding 0.6 °C. With help of these temperature and other associated variables in the experiments, it is possible to investigate the controlling factors leading to the EDW quantitatively and assess the influence of CP on the simulated EDW over the TP. Before this is done, it is useful to explore the reasons behind the biases of the simulations and the differences among them.

Biases in the simulations and their differences
Cold bias over the mountainous regions seems to be a persistent feature in RCMs as noted by previous studies Giorgi et al. 2004;Mishra 2015), and such bias over the TP may be partly explained by the overestimated precipitation ( Supplementary Fig. 4). As the resolution in the RCM is not high enough to fully resolve the steep terrain and associated local processes (e.g., valley wind, orographic drag), which leads to excess water vapor transport towards the TP and finally cold biases through snow-moisture-evaporation feedback (Haslinger et al. 2013;Lin et al. 2018). Besides, models without definite consideration of vegetation canopy in their surface-albedo calculation also overestimate the albedo of surfaces with 100% snow cover and SAF (Qu and Hall et al. 2007), which can exacerbate cold biases (Minder et al. 2016(Minder et al. , 2018. On the other hand, the Noah LSM will consider the snow melting process when the snow surface temperature becomes above freezing. Therefore, the underestimated temperature in WRF over the TP is accompanied by an expansive and long-existence of snow cover, which in turn leads to cold biases over the TP. These underestimations can be improved by applying satellite-observed albedo (Meng et al. 2018) or more realistic vegetation parameters  in the model, highlighting the necessity of correctly treating land surface processes over the TP. Furthermore, the convection-permitting models, which explicitly resolve convective heat and moisture and remove   the need for CPs (Fosser et al. 2015(Fosser et al. , 2017, have become important tools for improving the simulation of climate over the TP Gao et al. 2020;Zhou et al. 2021).
Because of the lower air density at high altitudes, the surface air temperature over the TP is mainly controlled by the surface temperature (Ts). Hence, the surface energy budget method proposed by Lu and Cai (2009) is applied in our study to diagnose the main factors leading to the differences among the simulations over the TP. We decompose the temperature differences using Eq. (3) and display these discrepancies of surface energy components and their sum in Fig. 8. Here, we take the climatological mean during 1980-2009 from the multi-model mean as the base state and the deviation of the five experiments from the based state as the perturbations. Due to the inhomogeneous distribution of spread, the TP is separated into two sub-regions, eastern TP (east of 90° E) and western TP (west of 90° E). The largest term contributing to the simulated Ts diversity is the SAF for the two sub-regions during all seasons, and its maximum contribution occurs in summer over the western TP and in spring over the eastern TP, with the value reaching up to 4.69 °C and 2.46 °C, respectively. The contribution of THF is always opposite to that of the SAF, which weakens the inter-model differences induced by SAF. The contributions Fig. 8 Differences in temperature between the five experiments and multi-model mean over Eastern TP (a) and Western TP (b). For the horizontal axis, the first item is the T2m, the second item is the surface temperature (Ts), the third to eighth items are the partial temperature differences due to surface albedo feedback, the change in cloud radiative forcing, the clear-sky short wave radiation, the clear-sky longwave radiation, the change in ground heat storage and change in turbulent heat flux, respectively. The last item is the sum of the third to eighth items. In addition, two items including ground heat storage and turbulent heat flux are multiplied by − 1 for convenient comparison of SW clr to the simulated temperature spread are negligible for all seasons, while the effects of other contributors display seasonal dependence. For example, the contributions of CRFs can't be ignored over the western TP during winter, while the LW clr plays an important role during summer. Also, the relative contribution of CRFs is larger over the eastern TP during summer than that over the western TP, and somehow offset the disagreement induced by the LW clr , which reduces the spread over the eastern TP.
As the SAF over the TP is quite sensitive to the existence of snow, the discrepancy in the simulated snow cover fraction (SCF) among the experiments is shown in Fig. 9. Similar spatial distributions between the spread in SCF and temperature are observed, with larger inter-model variability locating over the high altitudes in summer and autumn, and occurring along the southern edge of TP and southeast parts in spring and winter (Fig. 9). The SCF exhibits a high spatial heterogeneity over complex terrains and the variability of snow accumulation is strongly affected by the interaction of mountain ranges with the large-scale circulation causing orographic precipitation (Stoelinga et al. 2013). Due to the different choices of CP, large disagreements can be found in the simulated precipitation among the five experiments over the TP, especially during summer. It is found that the large spread in the simulated SCF always occurs over the regions with diverse precipitation simulations accompanied by T2m below 0 °C (Fig. 9). Such results support the findings of He et al. (2019) and Gao et al. (2020), who prove that substantial differences in the precipitation amount have a significant influence on the simulated snow water equivalent and SCF. Therefore, different precipitation amounts induced by CPs can affect the simulation of the SCF and SAF, which ultimately exerts a profound influence on temperature simulation. However, the disagreements in SCF among the experiments could not fully account for the simulated temperature spreads. For example, there is a larger spread in temperature simulation in spring over the western corner of TP (located over 79-85°E, 35-36°N), where little spread in SCF is present. This may be associated with the age of snow, which is an important parameter for snow albedo in the Noah LSM, since CP not only determines how much precipitation falls but also when it falls (Lynn et al. 2004). Thus, there is an urgent call for improving the SAF and precipitation simulation over the TP, in order to reduce the bias and uncertainty in temperature simulation.

Diagnosis of mechanisms responsible for EDW
Various physical mechanisms and processes, such as SAF (Giorgi et al. 1997;Palazzi et al. 2017), cloud feedback (Duan and Wu 2006), longwave water vapor feedback (Rangwala 2013;Rangwala et al. 2009), and aerosol feedback (Lau et al. 2010), have been proposed to explain EDW in previous studies, and most of them use Pearson correlation method to identify the divers of EDW. However, this method cannot fully quantify the contribution of each process to EDW. Thus, Eq. (3) is also used to investigate the relative importance of each surface energy component to EDW over the TP. Consistent with the above, we consider the first ten years as the base state and the last ten-year as the perturbation state, then decompose the temperature change into contributions from changes in six surface energy components (Fig. 10). Because of the large diversity among simulations in summer and pronounced EDW in winter, we hereafter focus on these two seasons for brevity.
During summer, the temperature change induced by SAF is close to zero below 3500 m, and then significantly increases with elevation up to around 5000 m in most of the Fig. 9 Spread in the simulated precipitation and snow cover fractions due to different cumulus parameterization schemes. Regions with T2m above zero are marked by black dots in the top row simulations, leading to the simulated EDW between 4000 and 5000 m over the TP. However, the contributions of SAF are partly canceled out by the THF, which presents opposite signs to that of SAF with weaker magnitudes. Other factors, such as CRFs and LW clr , play a relatively minor role in shaping the simulated EDW. Note that the elevationdependent structure of LW clr -induced warming exhibits a similar pattern as the Ts change, maximizing first at middle elevation and then at high elevation, and these two warming rates even present comparable magnitudes over the lowmiddle elevations, implying the important role of LW clr in determining the altitudinal variation of Ts over the TP. In addition, the contributions of SW clr and Q are negligible in all simulations.
During winter, there is a generally good correspondence between the structure of SAF-induced temperature change and EDW, supporting the hypothesis that SAF exerts the dominant control on EDW over the TP (Palazzi et al. 2017). The maximum SAF-induced warming up to 0.4-2.1 °C occurs around 4500-5000 m, where the snowline always lies . The profile of LW clr -induced warming shows a modest increase in temperature with elevation above Fig. 10 As in Fig. 7, but for temperature change due to each surface energy component during summer (a) and winter (b) 3000 m, again demonstrating the secondary role of LW clr in shaping simulated EDW over the TP. From the perspective of the surface energy budget, the SW clr has an insignificant role in generating EDW, while other factors tend to offset the EDW over the TP.
During other seasons, the vertical structure of temperature change can be largely understood by the changes in SAF ( Supplementary Fig. 5). The contributions of SAF to EDW between 2500 and 4000 m during spring and autumn are almost canceled out by THF, leading to little altitudinal variations in temperature change. Changes in SAF and CRFs result in a clear negative EDW (reduced warming at the highest elevations) over regions above 4000 m during spring, and this may be the reason that the models fail to reproduce the observed EDW over the high altitudes during spring. Those results demonstrate that further improvement of the components affecting the SAF and CRFs in the model is urgently needed to better understand the mechanism of EDW.
In general, SAF plays a primary role in controlling simulated EDW over the TP, and this may be explained by the altitudinal dependence of the loss of snow cover under climate system warming, which shows a similar profile with the SAF with the correlations between 0.88 and 0.99 during summer and between 0.98 and 0.99 during winter (Fig. 11a,  b). The altitudinal dependence of the loss of SCF leads to an elevation-dependent reduction in surface albedo and increases in absorbed solar radiation and ultimately generates EDW via SAF (Pepin et al. 2015). The effect of SAF in shaping EDW is also confirmed by sensitivity experiments over the Rocky Mountain wherein the SAF is artificially suppressed (Minder et al. 2018). They found that without an active SAF, a more uniform pattern warming is generated and the EDW is mostly eliminated.
Several studies find a statistical relationship between changes in LW clr and EDW, and they attribute the changes in LW clr to the normalized changes in water vapor amount (Rangwala et al. 2009(Rangwala et al. , 2016Palazzi et al. 2017). While in our simulations, the LW clr appears to play a secondary and relatively minor role in shaping EDW over the TP, and the variance of LW clr is mainly caused by the changes in air temperature instead of water vapor amount alone ( Fig. 11c-f). Ts anomaly will always be associated with the variance of LW clr via efficient boundary layer mixing, and it can greatly impact LW clr through strong lower atmosphere feedback, which in turn drives surface warming (Vargas Zeppetello et al. 2019). These imply that the forcing and response are difficult to distinguish due to the short equilibration timescale, and that new techniques are still required to improve such analyses.
Cloud feedback has also been considered as a potential contributor to EDW (Duan and Wu 2005;Liu et al. 2009), as changes in cloud cover and cloud properties significantly affect the surface energy budget. Our results suggest that the disagreements in EDW above 4000 m during spring between simulations and observation may be closely related to the CRFs. While in other seasons, there is no obvious evidence for the role of CRFs in contributing to EDW. The CRFs in our study consists of two parts: the cloud short-wave radiative forcing (CRF_SW) and long-wave radiative forcing Fig. 11 As in Fig. 7, but for snow cover fraction (SCF) changes, normalized changes in water vapor and lower-level air temperature changes. The normalized changes refer to the ratio of the changes between the two periods and the climatology mean for [2000][2001][2002][2003][2004][2005][2006][2007][2008][2009] (CRF_LW), and their contributions may be different. As is shown in Fig. 12, there is a general increase in warming rates due to the CRF_LW during summer, contributing to the EDW. However, this EDW signal could be canceled out by the CRF_SW, which shows a weakened warming rate with elevations below 4500 m. Therefore, the CRFs play a relatively minor role in shaping the simulated EDW during summer. Little changes in cloud cover by all experiments (Supplementary Fig. 6) result in small changes in CRF_SW during winter (Fig. 12). Although the temperature changes due to CRF_LW decrease with altitudes between 2500 and 4500 m and increase with altitudes above 4500 m, the contribution of CRFs to EDW during winter can be neglected compared to that of SAF. It should be noted that these conclusions apply specifically to the TP and the model settings considered in this study. For instance, there is no strong relationship between the patterns of LW clr changes and warming over the Rocky Mountains (Minder et al. 2018). Thus, multimodel simulations are needed to identify the character and causes of EDW across different mountain regions.

The impact of Cumulus schemes on EDW
To understand the sensitivity of EDW to the CP choice, the contributions of each surface energy component to EDW among the simulations were compared (Fig. 10). The inter-model differences caused by surface energy components commonly rise with elevation and reach their maximum around 5000-5500 m during summer and around 4000-4500 m during winter, especially for temperature change induced by SAF. As discussed above, the different precipitation changes due to the choice of CPs yield to diversity SCF variations, with the largest spread from − 4.9 to 1.9% locating 5000-5500 m during summer and from − 6.2 to 1.9% occurring 4000-4500 m during winter (Fig. 11). Such discrepancies in SCF simulations and the accompanying different land-atmosphere interactions lead to inconsistent temperature simulations. While for the CRFs, a larger range among simulations is found around 4000-4500 m during summer and above 5500 m during winter. Moreover, the spread in simulated EDW among experiments appears to be explained by SAF for the middle-high altitudes and by LW clr for the low altitudes in summer and autumn, and by SAF for all altitudes in winter and spring.

Summary, discussion and conclusions
Based on an observational dataset and modeling results with five cumulus parameterization schemes (CPs), this study analyzes the simulated temperature distribution over the Tibetan Plateau (TP), in terms of climatology and elevationdependent warming (EDW) during 1990-2009. The major contributors to EDW and sensitivity of simulated EDW to Fig. 12 As in Fig. 7, but for the temperature changes induced by the cloud short-wave radiative forcing (CRF_SW, a and b) and long-wave radiative forcing (CRF_LW, c and d) during summer (left column) and winter (right column) CPs are also investigated by applying the surface energy balance decomposition method.
Although all experiments could reasonably capture the observed temperature distributions over the TP, systematic cold biases still exist. The commonly cold biases in models over the complex terrain throughout the world may be related to the weaknesses in the model physics (Giorgi et al. 1989) and observational uncertainty due to lack of dense observational networks (Walker and Diffenbaugh 2009), both of which needs to be improved to better understand the processes as well as future climate projections. The spreads in simulated temperature due to CP is larger in summer than in other seasons, and they commonly increase with elevation with the maximum occurring around 4000-5000 m. A clear EDW signal below 5000 m is found in the observations and simulations during winter, which is consistent with the results from in-situ data (Liu and Chen 2000) and satellite data (Guo et al. 2019;Qin et al. 2009). However, all simulations fail to reproduce the observed EDW above 4500 m in spring and generate an EDW in summer around 4000-5000 m where the observation exhibits no clear EDW.
Further analysis reveals that the surface albedo feedback (SAF) is a primary contributor to the simulated temperature spread, and its maximum contribution occurs in summer over the western TP and in spring over the eastern TP. It also acts as the essential driver to the EDW over the TP during winter, evidenced by the good correspondence between the change in SAF and EDW in the simulations. The clearsky downward longwave radiation (LW clr ) appears to play a secondary role in shaping simulated EDW. Note that the models' inability to reproduce EDW above 4500 m during spring is closely related to the SAF and the surface cloud radiative forcing (CRFs), implying that further improvement of the components affecting the SAF and CRFs in the model is urgently needed to better understand the mechanism of EDW. Furthermore, the magnitude and structure of simulated EDW are sensitive to the choice of CPs. Different CPs can affect the simulation of snow cover fraction by modulating the precipitation amount, which alters the simulated SAF and its effect on EDW. The choice of CPs could also influence the LW clr via changing the low-level air temperature. Moreover, the spreads in the simulated EDW among the experiments are mainly caused by the SAF for the middle-high altitudes and by the LW clr for the low altitudes in summer and autumn, and by the SAF for all altitudes in spring and winter.
Using a single gridded observational dataset to assess the performance of the model over the Tibetan Plateau may be considered as a limitation of this study, because gridded datasets are afflicted with uncertainties due to different interpolation method (Hofstra et al. 2008), assimilations (Mao et al. 2010), and the observations they adopted (Schmidli et al. 2001). The uncertainties in the gridded temperature dataset over the TP have been shown in supplementary Fig. 7 and previous studies (Wang et al. 2020a;Mishra 2015), and the greatest diversity is mainly located over the high altitudes above 5000 m where no in situ stations exist ( Supplementary  Fig. 7). Additionally, the sparse in-situ stations over the high altitudes in TP could lead to large uncertainty in the gridded datasets as well as their representation of EDW. Thus, other datasets, such as remotely-sensed observations and regional reanalyses, should be analyzed to better understand the character of EDW over the TP.
It needs to be emphasized that the variability in EDW over the TP could also be affected by the internal variability of large-scale circulations (Pepin et al. 2015). For instance, enhanced warming at high elevations over the European Alps is observed during positive Pacific Decadal Oscillation (PDO) winter in the late 1980s and early1990s, while this signal was reversed during years with negative PDO index (Beniston et al. 1996). Therefore, our conclusions apply specifically to the simulation periods, and future work considering decadal variability in large-scale circulations should be carried out.
The disagreements in the temperature simulation and character of EDW, especially over high altitudes where glaciers, snow, and permafrost are mainly concentrated, limit our understanding of the climate process and enhance the uncertainty in water resource management. Qu and Hall (2014) suggest that using observations to constrain modeled albedos would reduce the spread of SAF and ultimately eliminate uncertainty in predictions of regional climate change among the multi-GCMs. Including effects of snow impurities into the land surface model is another way to generate a more realistic simulation of snow albedo, surface energy budgets, and temperature (Oaida et al. 2015). Additionally, models with finer resolution would provide more reasonable snow depth distributions over the mountainous regions, because of their improved depiction of fine-scale topography and more accurate representation of climate processes (Terzago et al. 2014;Gao et al. 2020), especially these related to moisture transport (Lin et al. 2018). Under the new framework of CORDEX-FPS-CPTP (Convection-Permitting Third Pole, http:// rcg. gvc. gu. se/ cordex_ fps_ cptp/), a set of convection-permitting sensitivity experiments at the kilometer-scale over the TP will be carried out (Ou et al. 2020;Wang et al. 2020b), and their results are expected to provide a more comprehensive assessment of climate change impacts over the mountain regions and promote the understanding of the EDW.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.