1 Introduction

Mountains are regions which are most sensitive to climate change and where climate impacts represent a threat to essential ecosystem services (e.g., Beniston 2003; Viviroli et al. 2007). One of the important questions related to climate change in the mountains is whether they are warming more than adjacent lowlands or compared to the global mean, similarly to the warming amplification observed in the Arctic (e.g., Screen and Simmonds 2010; Serreze and Barry 2011).

In recent years, the number of studies that have analyzed elevation-dependent warming (EDW)—the altitudinal dependence of warming rates—has increased (see MRI 2015 for a comprehensive review on this topic). These studies differ in the type of data that have been employed—in-situ observations, satellite data, reanalyses or climate model simulations—in the considered mountain ranges, and in the methods of analysis used to identify and quantify EDW. A majority of studies on EDW is based on the analysis of in-situ observations and points toward a general amplification of the warming rates with elevation. However, there is not a general consensus and some studies exist that show no relationship between warming rates and the elevation or that highlight a more complex situation whose understanding still requires further investigations (MRI 2015).

On the observational side, a homogeneous and dense network of meteorological stations would be required to clearly document the rate and the spatial distribution of warming as a function of elevation, but this is very seldom the actual situation encountered in high-altitude regions. The number of high-elevation stations providing long-term records (longer than at least 20 years) is still not adequate to allow evaluating statistically significant temporal trends, which is the first step for the assessment and quantification of EDW. Furthermore, mountain observations are known to be biased by altitude, since most in-situ stations are installed in valleys rather than on mountain slopes and on the tops, with the risk of undersampling these regions where the possible consequences of enhanced warming are expected to be more serious than elsewhere. Regions above 5000 m above sea level (a.s.l.) are mostly unexplored (Lawrimore et al. 2011), even if stations installed at such high altitudes would be crucial to fully understand hydro-climatic processes in the mountains and to document the small-scale temperature variations driven by topography, slope, vegetation coverage and exposure. Monitoring of land surface temperatures from satellite is another possible approach for studying EDW. The clear advantage of satellite observations over in-situ station data is their homogeneous spatial coverage; the disadvantage is that their temporal coverage is still quite short for detecting climatic trends and their statistical significance. These data are also still poorly validated in high-elevation regions where cloud occurrence represents an obstacle for satellite monitoring and data interpretation. There have been very few studies which have investigated the evidence for EDW using satellite observations (as shown in MRI 2015 see, in particular, Table S1 of their supplementary information). Qin et al. (2009), for example, performed a study using land surface temperatures measured by the satellite-borne MODIS sensor over the Tibetan Plateau in the period 2000–2006 and found that the warming rate increases in the altitude band from 3000 to 4800 m above sea level, and becomes then quite stable above those altitudes, with a slight decline at the highest elevations.

Model simulations are not affected by many of the inadequacies inherent in all kinds of observations, such as sparseness and limited temporal extension of the data, and they represent useful tools to investigate the possible mechanisms responsible for EDW, both in historical simulations and future projections. Of course they are generally limited in spatial resolution and require observational data for validation, making it difficult to be sure that simulations are accurate enough to be confident on future projections. Most of the existing model studies on EDW are based on the use of global climate models (e.g., Rangwala et al. 2010; Rangwala and Miller 2012) even though the response to large-scale greenhouse gas forcing in topographically complex mountain regions could be more adequately captured using higher-resolution non-hydrostatic regional models. Few studies exist in which climate change experiments have been performed with very high-resolution models (Hawkins et al, 2012; Kinter et al. 2013) and typically meteorological timescales have been considered to simulate specific aspects of the hydrological cycle, such as rainfall extremes (Kendon et al. 2014) and the response of the cryosphere system to changes in the hydrological cycle (Rasmussen et al. 2014).

GCMs have been used to assess the mechanisms that drive the elevational dependence of surface warming trends in the Tibetan Plateau, one of the regions showing the most striking evidence of EDW in the observations (Yan and Liu 2014). Since temperature at the Earth’s surface is primarily a response to the energy balance, potential EDW drivers are the factors that affect the net flux of energy at the surface and include, but are not limited to, albedo, clouds, water vapour and the related feedbacks. In a recent model study, Rangwala et al. (2013) analyzed the output of several GCMs from the Coupled Model Intercomparison Project phase 5 (CMIP5) experiment to investigate whether the mountains in the latitude band between \(27.5^\circ \hbox {N}\) and \(40^\circ \hbox {N}\)—the Tibetan Plateau in Asia and the Rocky Mountains in North America—are projected to warm at a faster rate relative to the low-lying areas of that latitude band by the end of the twenty-first century for three different greenhouse gas emission scenarios. They found an amplification of the warming rates in the Tibetan Plateau, particularly true for the minimum temperatures in the cold season. They ascribed the high correlation of enhanced warming with elevation in that region to an increase in downward longwave radiation at higher elevations, associated with increases in water vapour content in the atmosphere. A previous study by Liu et al. (2009) used both observations and model projections for the Tibetan Plateau region and the results were in line with the aforementioned findings: an amplified positive trend in minimum temperatures at higher elevations, particularly in winter and spring, likely caused by a combination of cloud-radiation and snow-albedo feedbacks. Another recent study by Yan et al. (2016) investigates the mechanisms driving EDW in and around the Tibetan Plateau using simulations of the CCSM3 GCM and concludes that the effects of changes in the snow depth and in cloud cover in response to the increase in CO\(_2\) concentrations \((4\times \hbox {CO}_2\hbox { experiment})\) lead to greater heat storage at higher elevations and therefore to EDW in that region.

The study presented here aims to analyze elevation-dependent warming in a wider region than the Tibetan Plateau only, encompassing also the Hindu-Kush, Karakoram, Himalayan mountains and the surrounding lower-lying regions using a selected ensemble of CMIP5 GCMs. We analyze the climate warming simulated by the models for this region and its dependence on elevation in the historical simulations and in the future projections under a high-range emission scenario (RCP 8.5). The possible mechanisms leading to EDW in this region are identified in the individual models and in their multi-model mean (MMM) and a multiple linear regression model is used to evaluate their relative contribution to the simulated change in the minimum and in the maximum surface temperature. In short, we select as variables of interest for EDW those which present a significant correlation with elevation (once having verified that also warming rates vary significantly with the elevation) and whose dependence on the elevation is physically consistent with the sign of the elevational gradient of warming rates. Then we identify as actual EDW drivers the variables which, independently of elevation, still correlate with temperature changes at each grid point. This approach has the advantage that the spatial information inherent in the climate models is to some extent “preserved” and represents an alternative or complementary way of approaching the issue of identifying the EDW mechanisms with respect to other approaches found in the literature (e.g., Rangwala et al. 2016).

The paper is organized as follows. Section 2 describes the area of study and the model data and provides an overview of the methodology used in the analysis. Section 3 reports the results of the analyses performed to assess an elevation-dependent warming signal in the Tibetan Plateau-Himalayan region using the historical and projection simulations of the CMIP5 model ensemble considered in this study. Section 4 discusses the driving mechanisms associated with EDW through the use of the multiple linear regression model including a number of processes/variables as predictors, the predictands being the minimum and the maximum temperature change between the end and the beginning of the century, focusing on the twenty-first century. Section 5 finally discusses the results and concludes the paper. We include as additional online resource material a Supplementary Information (SI) file devoted to a closer inspection of the results obtained from the individual CMIP5 models. We think that this is important to provide a measure of the spread of the CMIP5 models in representing EDW in the Tibetan Plateau-Himalayas. As shown in other studies which analyzed GCMs with various purposes, given the large differences which can be found in the basic model outputs, a multi-model ensemble estimate should in fact be regarded with extreme caution (e.g., Tebaldi and Knutti 2007; Palazzi et al. 2015).

2 Study area, data and methodology

2.1 Study area

The area under study extends from \(70^\circ \hbox {E}\) to \(105^\circ \hbox {E}\) longitude and from \(25^\circ \hbox {N}\) to \(40^\circ \hbox {N}\) latitude covering predominantly the high-altitude regions of the Tibetan Plateau, with an average elevation exceeding 4500 m (a.s.l.), surrounded by the highest mountain peaks of the world such as the Himalayas, Hindu-Kush, Karakoram, Pamir and others, and also includes flatter regions of central Pakistan and northwestern India, south of the Himalayan chain (Fig. 1a). The Tibetan Plateau-Himalayan region is often referred to as the “Third Pole” of the Earth because it hosts the largest amount of snow and ice outside the polar regions providing a very large reserve of fresh water. This region has a unique topography and landscape compared to other high elevation regions of the world, and is prone to different circulation regimes (e.g., Palazzi et al. 2013; Filippi et al. 2014) which strongly affect its hydrological cycle and climate.

Fig. 1
figure 1

(Top panel) Topographic map of the study area from a high-resolution Digital Elevation Model (0.008 degrees resolution). (Bottom panels) CMIP5 model ensemble orography at \(2\times 2\) degrees spatial resolution: multi-model mean (left) and inter-model standard deviation (right)

Full size image

2.2 Data and methods

We analyze the output of twenty-seven models participating in the Coupled Model Intercomparison Project phase 5 (CMIP5), available from the Earth Science Grid Federation archive data portals (http://esgf.llnl.gov). We select only those GCMs providing the following set of variables which are relevant for the study of EDW: surface altitude (orog), daily minimum and maximum near surface air temperature (tasmin and tasmax, respectively), surface downwelling longwave radiation (rlds), surface downwelling shortwave radiation (rsds), surface upwelling shortwave radiation (rsus), near surface specific humidity (huss). We calculate the surface albedo as the ratio between the upward and downward shortwave radiation at the surface. For each GCM we analyze the output of only one ensemble member, “r1i1p1” (each CMIP5 model is run producing an ensemble including different members, such as “r1i1p1” or “r2i1p1” or so, where “r” is an abbreviation for realization, “i” for initialization method, and “p” for physics version). Each of the global climate models has associated orographic data. Table 1 lists the CMIP5 models used in this study, along with their horizontal resolution, the number of grid points in the study area, and a key reference. The models are ranked in order of decreasing horizontal resolution. Further model details and information on their configuration or features can be found in the PCMDI data portal (http://www-pcmdi.llnl.gov/) and in Chapter 9 of the latest IPCC Assessment Report (AR5, IPCC 2013).

Table 1 The CMIP5 models used in this study
Full size table

We analyze the historical simulation outputs (1870–2005) and the projection outputs (2006–2100) for the Representative Concentration Pathway emission scenario RCP 8.5 corresponding to an anthropogenic radiative forcing of 8.5 Wm\(^{-2}\) by the end of the twenty-first century (Riahi et al. 2011). We decided to analyze the model projections for only one scenario based on previous studies (e.g., Rangwala et al. 2013) showing that the spatial patterns of the seasonal temperature changes between the end and the beginning of the twenty-first century in the Tibetan Plateau region are similar among different RCPs, with differences found just in the magnitude of the changes (being strongest in the most extreme RCP 8.5 scenario).

Most analyses presented in this paper refer to the elevation-dependent warming signal found in the multi-model mean (MMM) of the considered ensemble. The MMM is calculated after regridding all individual model outputs onto a common \(2\times 2\) degrees resolution grid before computing the average. The Supplementary Information shows and discusses some of the results obtained using the outputs of the individual models at their native resolution.

To investigate possible seasonal differences in the EDW signal suggested by previous studies on this topic (e.g.,Rangwala and Miller 2012; Rangwala et al. 2013 and references therein), the analysis presented here refers to seasonal averages, using the standard definition of the seasons for the Northern Hemisphere mid-latitudes: winter (December–February, DJF), spring (March–May, MAM), summer (June–August, JJA), and autumn (September–November, SON). We also analyze the warming and EDW signal found in the minimum and maximum temperatures separately, because different mechanisms which are possible drivers of EDW can be at work during day and night.

To quantify warming signals we calculate, for each model and for the MMM, the changes between the 30-year climatology of the minimum and maximum temperatures between the periods 1971–2000 and 1871–1900 (“historical” changes) and between the periods 2071–2100 and 1971–2000 (future changes under the RCP 8.5 scenario). The relationship between temperature changes and elevation is explored by using simple linear regressions and assessing their statistical significance. The latter is determined using a Monte Carlo method based on creating a large number (typically 1000) of surrogate series in which the data have been shuffled randomly with respect to the heights. The linear fits of these surrogate datasets allow to test against the null-hypothesis of no elevational gradients in the original series (Pollard et al. 1987; Schreiber and Schmitz 2000). A significance level of \(95\,\%\) is always used.

To study the mechanisms leading to EDW, we consider the seasonal changes of rlds, rsds, huss, and albedo and we analyze their dependence on the elevation. Their relative importance in driving the temperature change in the study area is assessed using a multiple linear regression model as discussed in Sect. 4.

For the observed elevation data shown in Fig. 1a, we use the ETOPO1 Global Relief Model from the NOAA’s National Geophysical Data Center (NGDC) available at http://www.ngdc.noaa.gov/mgg/global/global.html. ETOPO1 is a digital elevation model (DEM) available for the whole globe with a grid spacing of 1 arc-minute (approximately \(0.008^\circ \)). Figure 1b shows the average elevation of the CMIP5 ensemble (MMM, 2 degrees resolution) while the inter-model standard deviation is shown in panel (c). As expected, the models agree less with each other in the regions of strongest altitudinal gradients, for example along the Himalayan chain, while their standard deviation is smaller in the more homogeneous regions (central Tibetan Plateau and flatter regions of central Pakistan and northern India). Figure 2 shows the fraction of grid cells in each 800 m bin across the latitude and longitude range of the study area (Fig. 1) shown for the CMIP5 MMM (black line), the CMIP5 range (defined by 1 standard deviation greater and less than the mean, grey shading), and observed (DEM) elevation (red line). The figure also shows the mean elevation distribution of the three highest resolution models (CCSM4, CESM1-BGC, CESM1-CAM5, dashed blue line) and of the six lowest resolution models (MIROC-ESM-CHEM, MIROC-ESM, bcc-csm1-1, BNU-ESM, CanESM2, FGOALS-g2, dashed orange line). Square and star symbols are used to indicate, for each elevation bin, the maximum and minimum fraction of grid cells across the whole CMIP5 model ensemble. Of course, the higher the model resolution, the greater the number of high-elevation grid points in the region of interest. The highest-resolution GCMs have an elevation distribution that is closer to observations than the CMIP5 average. The CMIP5 MMM overestimates the fraction of grid points in the altitudinal bins at elevations between about 1000 and 4000 m above sea level, while it underestimates the fraction of grid points at lower and higher elevations.

Fig. 2
figure 2

Elevation distribution across the region. Fraction of grid cells in each 800 m bin across the Tibetan Plateau-Himalaya region (Fig. 1) shown for observations (red line), the CMIP5 MMM (black line), the CMIP5 range (grey shading), the GCMs with the highest number of grid points in the area (dashed blue line), and the GCMs with the lowest number of grid points in the area (dashed orange line). Empty squares and star symbols indicate the maximum and minimum fraction of grid point found across the whole model ensemble

Full size image

3 The simulated elevational dependence of surface warming

3.1 Historical period: 1871–2000

Figure 3 shows, for the four seasons, the difference between the average in the period 1971–2000 and the average in the period 1871–1900 of the minimum temperature (left panels) and of the maximum temperature (right panels) as a function of elevation for the MMM of the CMIP5 ensemble. The data points have been fitted with a linear regression over the whole altitude range (red line) and from 1500 m upwards (blue line); the slopes of the two regression lines (in \(^\circ \hbox {C km}^{-1}\)) are reported in each panel (a star symbol in parentheses indicates a statistical significant trend) and summarized in Table 2. Figure 3 shows that, for both tasmin and tasmax, all grid points show positive changes (i.e., warming) and the absolute value of these changes is larger at higher elevations (positive elevational gradients of warming rates). The elevational gradients are generally enhanced when assessed over an altitude range that excludes elevations lower than 1500 m (a.s.l.). The minimum temperature shows the strongest elevational gradients of warming in spring and autumn and the weakest in summer. For the maximum temperature the seasonal differences in the elevation dependent warming are, in general, less important than for the minimum temperature. Overall, the EDW signal is stronger for tasmin than for tasmax in spring, while the opposite is found in summer. Figure 3 also highlights that, beyond considerations about the elevational gradients of warming rates, the minimum temperatures exhibit a stronger absolute change than the maximum temperatures. This finding is consistent with previous studies showing the asymmetric nature of warming rates during daytime and nighttime over the Tibetan Plateau (e.g., Liu et al. 2009). The seasons showing the strongest changes – averaged over the study area—of both tasmin and tasmax during the historical period are DJF and MAM.

Fig. 3
figure 3

Change between the period 1971–2000 and the period 1871–1900 of the minimum temperature (left panels) and of the maximum temperature (right panels) as a function of surface elevation for each season and for the CMIP5 MMM (\(2\times 2\) degrees horizontal resolution). The slope of linear regression (\(^\circ \hbox {C}\hbox { km}^{-1}\)) is indicated (see text for details); a star in parentheses indicates the statistical significance of the elevational gradients

Full size image
Table 2 Summary of the modeled elevational gradient (slope of the linear regression, in \(^\circ \hbox {C km}^{-1}\), describing the temperature changes as a function of the elevation in the twentieth and twenty-first centuries) in the Tibetan Plateau-Himalayas for the multi-model mean of the CMIP5 ensemble. Values in parentheses indicate the slope of the linear regression for the data points fitted from 1500 m upwards
Full size table

A similar analysis is performed for each model separately at its native spatial resolution: the models exhibit noticeable differences with each other (not shown here), such that some of them agree very well with the MMM while others do not and even exhibit an elevational gradient of warming rates of opposite sign (indicating decreased warming as elevation increases). The slope of the linear regression between the minimum (maximum) temperature changes and elevation for the individual models is reported in Table S1 (S2) of the Supplementary Information. We find that the inter-model spread, measured as one standard deviation greater and less than the mean, overall increases as the elevation increases, particularly from 3000 m above sea level upward (shown in Figures S1 and S2 of the SI, as well as in Fig. 4 discussed later).

Fig. 4
figure 4

Minimum and maximum temperature change between the 1971–2000 climatology and the 1871–1900 climatology as a function of surface elevation for the multi model mean data averaged in bins with width 150 m, for the four seasons. The blue (red) lines show the multi-model mean of the GCM ensemble while the line-filled (solid-filled) areas represent the range of variability of the models measured as one standard deviation greater and less than the MMM for the minimum temperature (maximum temperature)

Full size image

Figure 3 suggests that a simple linear fit describing the relationship between the temperature changes and the elevation over the whole range of altitudes is likely to be a simplistic approach. A more complex situation, in fact, arises in which the (either minimum or maximum) temperature change tends to slightly decrease as a function of the elevation or to remain almost constant, depending on the season, until about 1500 m above sea level. It then increases with elevation giving rise to a positive and statistically significant slope. Finally, above about 4500 m (a.s.l.), the elevational gradient reverses indicating a reduction of the warming rate with elevation (in agreement with results found by Qin et al. 2009). In order to make the non-linearity of the relationship between the temperature changes and the elevation clearer, Fig. 4 shows the minimum and maximum temperature change between the 1971–2000 climatology and the 1871–1900 climatology as a function of the surface elevation in the CMIP5 ensemble for the four seasons. Data have been averaged in bins with width 150 m, to compute the elevational distribution of the MMM and the range of variability in the models (inter-model spread, calculated as one standard deviation greater and less than the MMM). The figure highlights the complex regime mentioned above, which is more evident in winter, spring, and autumn particularly for the minimum temperature: the elevational gradient of the minimum temperature change is observed to vary and even be reversed at altitudes around 1500–2000 m above sea level, as also discernible from the slopes of the linear fit calculated over the two different altitude ranges shown in Fig.  3 (red and blue lines and figures in the legend). Figure  4 also reiterates, in a very clear way, that the temperature changes between the last three decades of the twentieth century and the last three decades of the nineteenth century are more pronounced for the minimum than for the maximum temperature, corroborating previous studies that identified the asymmetric nature of warming rates in the Tibetan Plateau (e.g., Liu et al. 2009). A stronger warming trend in the daily minimum than maximum temperatures was also reported for other mountain areas, for example the Alpine region, as described in (Jungo and Beniston 2001).

3.2 Projected changes: 1971–2100

The same analyses performed for the historical period are repeated to assess the simulated elevational gradients of warming by the end of the twenty-first century. Figure 5 shows the same as Fig. 3 but the minimum and maximum temperature changes are calculated as the difference between the 2071–2100 climatology and the 1971–2000 climatology using the model projections under the RCP 8.5 scenario (from 2006 onwards). A marked elevation-dependent warming signal is found for the minimum temperature in all seasons; the strongest signal for the maximum temperatures is found in autumn.

Fig. 5
figure 5

Change between the period 2071–2100 and the period 1971–2000 of the minimum temperature (left panels) and of the maximum temperature (right panels) as a function of surface elevation for each season and for the CMIP5 MMM (\(2\times 2\) degrees horizontal resolution). The slope of linear regression (\(^\circ \hbox {C}\hbox { km}^{-1}\)) is indicated (see text for details); a star in parentheses indicates the statistical significance of the elevational trend

Full size image

The elevational gradients of warming rates in the scenario simulations are higher than in the historical period: the slopes calculated over the whole range of altitudes (red lines and the corresponding figures in the legend of Fig. 5) are from about 10 (minimum temperatures in spring and maximum temperatures in winter) to about 30 (minimum temperatures in winter and maximum temperatures in summer) times higher than those evaluated for the historical period. The absolute value of warming in itself is considerably higher than in the historical model simulations.

As already pointed out for the elevational gradients of warming rates in the twentieth century, the relationship between the temperature change and the elevation emerging from Fig. 5 is far from simply linear, which is also well visible in the top panels of Fig. 6 (the same as Fig. 4, but for the future scenarios). The two-fold behaviour observed in the historical period is emphasized in the future projections and the stronger absolute value of warming likely allows to better separate the actual signals from the background noise. The inter-model spread of the minimum and maximum temperature changes at different elevations is relatively large, larger in the scenario simulations than in the historical ones (Fig. 6 compared to Fig. 4). Interestingly, the inter-model spread tends to increase as the elevation increases, particularly from \(\sim\)3000 m above sea level upward (see also Figure S2 of the Supplementary Information). This highlights that model simulations tend to be in less agreement with each other at higher altitudes, possible owing to problems in representing the complex orography of high-elevation regions. Figure 6, as well as Fig. 4, shows the results for the MMM and the inter-model standard deviation, but the same conclusions can be drawn if the multi-model median and the percentiles are used instead.

Fig. 6
figure 6

Minimum and maximum temperature change between the 2071–2100 climatology and the 1971–2000 climatology as a function of surface elevation for the multi model mean data averaged in bins of width 150 m (top panels) and as a function of the mean temperature for the multi model mean data averaged in bins of width 1\(^\circ \hbox {C}\) for the four seasons (bottom panels). The blue (red) lines show the multi model mean of the GCM ensemble while the line-filled (solid-filled) areas represent the range of variability of the models measured as one standard deviation greater and less than the MMM for the minimum temperature (maximum temperature)

Full size image

The two-fold behaviour mentioned above is clearly evident in the bottom panels of Fig. 6, where the minimum and maximum temperature changes are plotted as a function of the mean surface temperature, rather than as a function of the mean elevation. The mean temperature is calculated as the average between tasmin and tasmax also averaged over the 30-years long time periods. In particular, the wintertime signal in the minimum temperature is strongly emphasized and the switch between the two regimes occurs close to the \(0^\circ \hbox {C}\) isotherm. It is interesting to note that these considerations hold true not only for the MMM and the inter-model spread, but also for each individual model (not shown here), which gives robustness to these findings. The existence of different temperature responses in regions with temperatures above and below the freezing point observed for tasmin in winter seems to occur also in the other seasons and for tasmax (Fig. 6) even if the signal is generally less evident and the switch between the two regimes is smoother and less localized across the zero-degree isotherm. The two-fold behaviour in the relationship between temperature change and elevation is not evident during summer: the bottom panels of Fig. 6 suggest that the two regimes might not exist in this season simply because there are almost no points with mean temperatures below zero.

Interestingly, Fig. 6 shows that daily asymmetry in warming rates observed in the historical period is strongly reduced and almost disappears in the projections. The winter season represents the only notable exception: from about 1500 m upward, in fact, the minimum temperatures are projected to warm much more than the maximum temperatures.

Looking at the individual model outputs presented in the SI (Tables S3 and S4), it is interesting to note that, unlike in the historical simulations, almost all models simulate positive elevational gradients of temperature changes in the future and in the various seasons, particularly true for the maximum temperature (Table S4). Despite that, the inter-model spread remains very high and it increases at higher elevations (Figures S3 and S4). Table 2 summarizes the slopes of the linear regression between the minimum/maximum temperature changes and the elevation for the CMIP5 MMM in the future projections (as observed for the historical simulations, all values are statistically significant).

4 Mechanisms

The analyses reported in Sect. 3 indicate that the CMIP5 models and their multi-model mean depict a situation in which warming is amplified with elevation and this is found both in the historical simulations and in the projections under the RCP 8.5 scenario, the signal being much stronger in the latter case. Also, the link between temperature changes and elevation is described by a more complex relationship than a simple linear one and two main regimes can be distinguished, corresponding to temperatures above and below the zero-degree isotherm. This would suggest a leading role of mechanisms exhibiting enhanced sensitivity around \(0\,^\circ \hbox {C}\) and therefore likely involving water phase changes, the presence/absence of snow and the snow-albedo feedback.

In this section we explore further the mechanisms leading to EDW in the Tibetan Plateau-Himalayas. For the sake of conciseness, we consider only the cases (i.e. variables, seasons and periods) showing the strongest EDW signal, namely the minimum temperature during winter and the maximum temperature during autumn in the RCP 8.5 scenario (see in particular Table 2). We analyze other CMIP5 model variables whose variations may be important for determining temperature variations and their dependence on the elevation. The variables which we focus on are albedo, surface downwelling longwave radiation, surface downwelling shortwave radiation, and near-surface specific humidity, which may affect the minimum and maximum surface temperature changes both directly and indirectly by means of, e.g., feedback mechanisms involving clouds. Changes of these variables between the end and the beginning of the twenty-first century (hereafter called \(\varDelta albedo\), \(\varDelta rlds\), \(\varDelta rsds\), and \(\varDelta huss\)) are computed in the same way as the minimum and maximum temperature changes and are considered as possible drivers of EDW. In addition, we include three other possible drivers, namely the normalized changes (i.e., the fractional changes relative to the averaged climatology between the mean in the years 1971–2000 and the mean in the years 2071–2100) in surface downwelling longwave and shortwave radiation (\(\varDelta rlds/rlds_0\), \(\varDelta rsds/rsds_0\)) and in near-surface specific humidity (\(\varDelta huss/huss_0\)) since they could be more effective in determining elevation-dependent temperature signals, as highlighted in previous studies (e.g., Rangwala et al. 2013; Rangwala et al. 2016). In particular, the authors of those studies pointed out that the normalized change in specific humidity is a more appropriate metric than the absolute change, because of the sensitivity of the downward longwave radiation to water vapour changes. For the same increase in atmospheric water vapour amount, indeed, the downward longwave radiation increases more, as percentage, when the mean state is drier (as in the high-elevation sites) and less when it is wetter (as in lower-elevation sites), leading to higher warming rates at higher elevations. Summarizing, we identify the following seven possible drivers of EDW: \(\varDelta albedo\), \(\varDelta rlds\), \(\varDelta rsds\), \(\varDelta huss\), \(\varDelta rlds/rlds_0\), \(\varDelta rsds/rsds_0\), \(\varDelta huss/huss_0\). In order to understand if these changes are actually driving EDW in the study area considered here, we check if they satisfy the following three conditions:

  1. 1.

    they have to exhibit a significant dependence on elevation, as the temperature change does

  2. 2.

    the sign of this dependence has to be consistent with the observed changes in surface temperatures, i.e. with the amplification of the warming rate at higher elevations

  3. 3.

    they have to spatially correlate with temperature changes even when the dependence on elevation is removed.

Since we select as variables of interest those which present significant correlations with elevation (condition 1), they will also all be correlated with temperature changes, which also have been found to correlate with altitude. Condition 3 is useful to identify those variables which, independently of elevation, still correlate with temperature changes. Condition 2 represents a check if there is a known physical mechanism which can link their variations with altitude with the altitudinal dependence of temperature changes.

The first step is to check for the existence of a significant relationship between \(\varDelta albedo\), \(\varDelta rlds\), \(\varDelta rsds\), \(\varDelta huss\), \(\varDelta rlds/rlds_0\), \(\varDelta rsds/rsds_0\), \(\varDelta huss/huss_0\) and elevation. Figure 7 shows, for each individual GCM and for the MMM, the correlation coefficients between each of the seven possible drivers and elevation during winter; the correlation coefficient between \(\varDelta tasmin\) and the elevation is also displayed for completeness. White boxes appear when the correlations are not statistically significant. The MMM shows that all correlations are statistically significant except those involving the changes of the surface downward shortwave radiation (both absolute and normalized). This is consistent with the signal emerging from the individual GCMs, showing a general consensus in the sign of the correlation coefficients involving all variables but \(\varDelta rsds\) and \(\varDelta rsds/rsds_0\), for which correlations of either signs occur. The same result is found in autumn (Fig. 8), the season in which the strongest EDW signal is observed in tasmax. Another interesting feature emerging from Fig. 7 is that, in spite of the significant correlation found in the MMM between \(\varDelta albedo\) and the elevation in DJF, the lowest resolution models (from IPSL-CM5B-LR to FGOALS-g2, except bcc-csm1-1) show no significant correlations instead, suggesting that their poorer depiction of the topography in the Tibetan Plateau-Himalayas hampers an accurate representation of surface processes in these models. This first-level analysis suggests not to include the changes (both absolute and normalized) of rsds in the successive steps, since they do not show a clear dependence on the elevation (i.e., they do not satisfy the first condition listed above).

Fig. 7
figure 7

Correlation coefficient between the following variables - \(\varDelta rsds\), \(\varDelta rlds\), \(\varDelta huss\), \(\varDelta albedo\), \(\varDelta rsds/rsds_0\), \(\varDelta rlds/rlds_0\), \(\varDelta huss/huss_0\), \(\varDelta tasmin\) - and the elevation for all CMIP5 models and for the multi-model mean (MMM) in DJF and in the study area

Full size image
Fig. 8
figure 8

The same as Fig. 7 but in SON. In this case the correlation between the change in tasmax and elevation is shown

Full size image

To verify the fulfilment of the second condition, the sign of the correlation coefficients has to be considered. Basic physical considerations imply that, in order to be consistent with the elevational gradient of the temperature change, changes of rlds and huss have to exhibit the same elevational dependence as the temperature change does, since an increase in these variables determines an increase in temperature as well. Vice versa, albedo changes have to exhibit an elevational gradient of opposite sign, since an increase in albedo causes less absorption of heat by the surface and a decrease in surface temperature. Figures 7 and 8 confirm that the temperature changes are positively correlated with elevation. This holds true for both \(\varDelta tasmin\) and \(\varDelta tasmax\). The figures also show that \(\varDelta rlds/rlds_0\) and \(\varDelta huss/huss_0\) are positively correlated with elevation, while \(\varDelta albedo\), \(\varDelta rlds\) and \(\varDelta huss\) are negatively correlated. This excludes from the ensemble of the possible EDW predictors the absolute changes of rlds and huss since their negative correlation with the elevation is not physically consistent with the enhanced warming at higher elevations.

The considerations made so far suggest that the three following variables may drive the positive elevational gradient of \(\varDelta tasmin\) in winter and of \(\varDelta tasmax\) in autumn: \(\varDelta albedo\), \(\varDelta rlds/rlds_0\), and \(\varDelta huss/huss_0\). However, although all of them vary with elevation and their elevational gradient has the “right sign”, their correlation with the temperature change could arise simply because they exhibit a dependence on the elevation as the temperature change does, without any causal correlation with temperature changes and the associated EDW. For this reason, we check for the third condition to be satisfied: \(\varDelta albedo\), \(\varDelta rlds/rlds_0\), and \(\varDelta huss/huss_0\) have to spatially correlate with the temperature change even when we “control for” elevation (i.e., when we remove the linear dependence of these variables and of the temperature change on the elevation) in order to be considered actual drivers of EDW. To this end, a multiple regression model is used where the three possible EDW drivers are used as predictors for the seasonal change in either tasmin or tasmax (the predictand). The resulting multiple linear regression model can be written as:

$$\varDelta T' = a_1 \varDelta {albedo}' + a_2 (\varDelta {huss}/{huss_0})' + a_3 (\varDelta {rlds}/{rlds_0})' + \eta$$
(1)

where the apex is used to indicate that both the predictors and the predictand are “altitude-detrended” by removing their linear fit on elevation and that they are standardized variables (each change is divided by its standard deviation over the whole spatial domain). \(\varDelta \hbox {T}\) indicates the seasonal change in either tasmin or tasmax. This approach allows to test all the possible combinations of the three predictors including all or part of them. Each combination leads to a different regression model. Overall, the possible models are \((2^n-1)=7\), with \(n=3\) predictors and include the three regression models with one single predictor, the three models with the combinations of two predictors and the model with all three predictors. The predictive ability of the seven regression models is assessed by examining the proportion of the variance they can explain, which can be measured by the coefficient of determination \(R^2\): the closer \(R^2\) is to 1, the better the model predictive ability is. Since by construction the regression model including all predictors is associated with the highest value of \(R^2\), the Akaike information criterion corrected for finite sample sizes is used to measure the relative quality of the models (AICc, Burnham and Anderson 2003). This metric allows to rank the regression models with different number of predictors penalizing those with more and favoring those with less predictors. The lower the AICc, the better the model. The seven multiple regression models are applied to all 27 GCMs and to their MMM.

Fig. 9
figure 9

Top coefficient of determination \(R^2\) of the seven regression models for each GCM and for the MMM. The averaged value of \(R^2\) across the GCMs is also reported \(\left< R^2 \right>\). Bottom the seven regression models obtained through the combinations of the three predictors, a black box indicating which predictor is included in each model

Full size image
Fig. 10
figure 10

The same as Fig. 9 but for \(\varDelta tasmax\) in SON

Full size image

The resulting coefficients of determination are graphically shown in Fig. 9 (prediction of the minimum temperature change in DJF) and in Fig. 10 (prediction of the maximum temperature change in SON). For completeness, the average of the \(R^2\) values across the different GCMs is also reported, \(\left< R^2 \right>\). The bottom part of Figs. 9 and 10 indicates which predictors are “switched on” (black squares) or “switched off” (blank) in each regression model. To facilitate the interpretation of the results, Table 3 summarizes the results for the MMM. In this table, the seven regression models are ranked in order of increasing AICc values. As a starting point, we consider the prediction of \(\varDelta tasmin\) in DJF and focus on the MMM of the CMIP5 models for conciseness (Fig. 9 and upper part of Table 3). In this case, \(\varDelta albedo\) emerges as the driver having a leading role for EDW. In fact, the four models including \(\varDelta albedo\) as a predictor show the highest values of \(R^2\) among the seven regression models for the MMM. In particular, among the three single-predictor models, the one with \(\varDelta albedo\) shows the highest \(R^2\). The three multi-predictor models including \(\varDelta albedo\) in conjunction with any other variable are capable of accounting for most of the variance of the predictand (more than \(75\,\%\)). The other predictors are generally less relevant than \(\varDelta albedo\) but nonetheless they play a role in the performance of the regression model at hand. As mentioned, indeed, the multi-predictor models including \(\varDelta huss/huss_0\) and/or \(\varDelta rlds/rlds_0\) in addition to \(\varDelta albedo\) are the ones showing the best (i.e. the lowest) AICc values.

Table 3 For each of the seven regression models, the top (bottom) part of the table shows the values of the coefficients \(a_1\), \(a_2\), and \(a_3\) of the regression (see Eq. 1) for the prediction of \(\varDelta tasmin\) in DJF (\(\varDelta tasmax\) in SON). For each model the coefficient of determination \(R^2\) and the AICc are shown. The models are ranked according to the AICc values. The results refer to the MMM
Full size table

Similar considerations can be drawn for the prediction of \(\varDelta tasmax\) in SON (Fig. 10 and lower part of Table 3). Also in this case the four regression models including \(\varDelta albedo\) as a predictor show the highest values of \(R^2\) for the MMM. Again, the other predictors are less relevant than \(\varDelta albedo\) but their presence still increases the variance explained by the regression model. However, in this case lower \(R^2\) values are found. This holds true for both the MMM as well as for most of the individual GCMs. The only exception is for a group of GCMs including CNRM-CM5 and the GFDL and GISS “model families”, for which the \(R^2\) values of the regression models including \(\varDelta albedo\) are rather high. The low \(R^2\) values shown by the regression models during SON indicate that the three predictors are not able to explain much of the variance of \(\varDelta tasmax\) in this season, and that some other relevant drivers may be at work. In this regard, it has to be noted that the predictors considered here have been chosen also considering processes and variables suggested by previous studies which have mainly focused on minimum temperatures, the winter season, and slightly different areas (see for example Liu et al. 2009; Rangwala et al. 2016). As such, lower \(R^2\) for prediction of the maximum temperature changes in other seasons can be expected and are not so surprising.

It is worth noting that the GCMs exhibit very large differences with each other (see Figs. 9, 10). This significant inter-model spread suggests that the processes associated with the positive temperature anomalies and driving the EDW signal in different GCMs are treated differently. It is interesting to note that the GCMs belonging to the same “model family” provide more similar signals with each other than with other independent models (see Knutti et al. 2013 for an overview of the lineage connections among the climate models used in this study). The reasons for such inter-model discrepancies are not explored in the present study, but could be related, among the others, to a different description and treatment of aerosols and clouds that can affect the radiation transfer in the atmosphere.

5 Conclusions

Elevation-dependent warming (EDW) in the Tibetan Plateau-Himalayan region has been investigated using twenty-seven global climate models of the CMIP5 ensemble. Results indicate higher warming rates at higher elevations in this region, both in the historical model simulations and in the projections under the RCP 8.5 emission scenario. The EDW signal is amplified in the future simulations with respect to the historical period, and it is particularly strong for the minimum temperatures in winter and spring and for the maximum temperatures in summer and autumn. While being in many aspects consistent with previous modeling studies (e.g., Rangwala et al. 2013; Liu et al. 2009 and others) which analyzed EDW in the high-altitude regions of the Tibetan Plateau (and therefore in a smaller region than that taken into account in our study), in this paper we emphasize specific aspects not much discussed previously. First of all we show that the relationship between the maximum and minimum temperature changes and the elevation is not simply linear. At least two different regimes can be identified occurring at temperatures below and above the zero-degree isotherm: the regions with temperatures below freezing show a much stronger warming than the regions with temperatures above, suggesting that the phase of water and/or the presence of snow play a key role.

To better understand the mechanisms that drive EDW in this region, a number of variables that may be relevant in determining the stronger warming of higher elevations have been investigated. These variables are: \(\varDelta albedo\), \(\varDelta rlds\), \(\varDelta rsds\), \(\varDelta huss\), \(\varDelta rlds/rlds_0\), \(\varDelta rsds/rsds_0\), \(\varDelta huss/huss_0\). Only three of them, however, – \(\varDelta albedo\), \(\varDelta rlds/rlds_0\), \(\varDelta huss/huss_0\) – satisfy the necessary conditions to act as possible EDW drivers: (1) to exhibit a significant dependence on the elevation, as the temperature change does; (2) to be correlated with the elevation in such a way that they can support the amplification of warming rates at higher elevations; and (3) to spatially correlate with the temperature change even when the dependence on the elevation is removed. The relative importance of these variables in determining temperature changes has been assessed using a multiple regression model. \(\varDelta albedo\) emerges as the driver having a leading role for EDW. For both prediction of \(\varDelta tasmin\) in DJF and \(\varDelta tasmax\) in SON, in fact, the four models including \(\varDelta albedo\) as a predictor show the highest values of explained variance \((R^2)\) among the seven regression models for the MMM. In particular, among the three single-predictor models, the one with \(\varDelta albedo\) shows the highest \(R^2\). The three multi-predictor models including \(\varDelta albedo\) in conjunction with any other variable are capable of accounting for most of the variance of the predictand (more than \(75\,\%\) in winter). The other two predictors are generally less relevant than \(\varDelta albedo\) but nonetheless they have a non-negligible impact on the performance of the regression model. More extensive analysis is needed to include other possible drivers of EDW, such as aerosol particles or clouds. Another important aspect would include the analysis of the potential correlations between the drivers—for example between changes in downward shortwave radiation and changes in downward longwave radiation, both associated with the presence and type of clouds, or changes in downward shortwave radiation and consequent changes in surface albedo—because they might lead to compensation or amplification of the individual effect of the drivers on the enhanced warming rates at higher elevations.

More research is needed to understand the inherent complexity of mountain regions in general and of the elevation-dependent warming in particular, especially given the environmental and socio-economical impacts of a changing climate in mountains. The model studies performed so far to investigate the processes associated with EDW have employed global climate models mainly. The change in surface albedo turns out to be a main driving factor for EDW in the Tibetan Plateau-Himalayas, suggesting the urgency of further developing models of those Earth system components which affect albedo, including (but not limited to) snow cover and mountain glaciers. Proper simulation of snow cover requires appropriate simulation of snowfall and hence of clouds and cloud dynamics and processes, which are subgrid-scale phenomena in state-of-the-art GCMs. Also, the implementation of proper parameterizations for the dependence of snow albedo on snow age, depth, terrain characteristics, aerosol deposition and others is of major importance. The study by Essery (2013), for example, found that snow-covered albedo simulations depend much more on the choice of the parameters than on the specific scheme implemented in the land-surface model itself. They suggested that a correction of the albedo-related parameters to constrain the simulated albedo on the basis of satellite observations could reduce the inter-model spread for snow-covered albedo and should be applied in the next generation of CMIP models. Modelling glaciers, their dynamics, and the effects of climate change on them is quite complex, too. Glacier expansion and retreat depend on the balance between accumulation and ablation and, therefore, winter snowfall and summer temperatures are key ingredients to assess the future glacier conditions. Previous studies (e.g., Su et al. 2013) highlighted that there is a deficiency in the way the CMIP5 models represent snow and ice albedo, which is mentioned as one of the main causes for the cold bias that still affects many GCM simulations. The largest cold biases just appear in areas with varying topography and permanent ice (Mao and Robock 1998). The standard output of the CMIP5 GCM simulations includes information on the fraction of grid cell occupied by permanent ice (i.e, glaciers), as well as other snow-related variables like snow cover, snow depth or snow water equivalent, snow melt fluxes, and others. These variables are land fields and are produced by the land vegetation models that are coupled to the other model components in the state-of-the-art Earth System models. Therefore, besides improving the parameterizations for subgrid-scale processes, also improving the land-surface models would lead to a better description of the high-mountain cryosphere system. Besides that, the use of finer resolution models would be useful to depict the complex topography of mountain regions in a more realistic way and therefore improve the way the models represent the changes of snow at ground (e.g., Terzago et al. 2014). The same considerations apply for other radiation-related processes (involving clouds in particular) and other surface processes which would benefit from the use of high-resolution climate models (possibly including non-hydrostatic approximations for the atmospheric vertical motions). At the same time, however, analyzing the capabilities of the GCMs in reproducing EDW and its driving mechanisms is important, since global climate models are the large-scale drivers of regional models and are crucial for many climate change assessments and impact studies (IPCC 2013).

Our results do not still have a solid observational counterpart because of the inadequacy or even the lack of ground-based station networks in (some portions of) these regions. More in-situ observatories, complemented by satellite data and appropriate analysis methods, would be necessary to measure EDW and to evaluate and validate the CMIP5 model performances, in order to be more confident of their projections on future elevation-dependent warming and its consequences.