Comparison of elevation-dependent warming and its drivers in the tropical and subtropical Andes

Mountain regions have been recognized to be more sensitive to climate and environmental changes, and in particular to global warming. Several studies report on elevation-dependent warming (EDW), i.e., when warming rates are different in different altitude ranges, particularly focusing on the enhancement of warming rates with elevation. The Andean chain proved to be a relevant climate change hot-spot with positive temperature trends and a widespread glacier retreat over the recent decades. To assess and to better understand elevation dependent warming in this mountain region and to identify its possible dependence on latitude, the Andean Cordillera was split into five domains, three pertaining to the tropical zone and two pertaining to the Subtropics. Further, for each area the eastern and western faces of the mountain range were separately analyzed. An ensemble of regional climate model (RCM) simulations participating in the Coordinated Regional Climate Downscaling Experiment (CORDEX), consisting of one RCM nested into eight different global climate models from the CMIP5 ensemble was considered in this study. EDW was assessed by calculating the temperature difference between the end of the century (2071–2100) and the period 1976–2005 and relating it to the elevation. Future projections refer to the RCP 8.5 high-emission scenario. Possible differences in EDW mechanisms were identified using correlation analyses between temperature changes and all the variables identified as possible EDW drivers. For the maximum temperatures, a positive EDW signal (i.e. enhancement of warming rates with elevation) was identified in each side of both the tropical and subtropical Andes and in all seasons. For the minimum temperatures, on the contrary, while a positive EDW was identified in the Subtropics (particularly evident in the western side of the chain), the Tropics are characterized by a negative EDW throughout the year. Therefore, the tropical boundary marks a transition between discordant EDW behaviours in the minimum temperature. In the Tropics and particularly in the inner Tropics, different EDW drivers were identified for the minimum temperature, whose changes are mostly associated with changes in downward longwave radiation, and for the maximum temperature, whose changes are mainly driven by changes in downward shortwave radiation. This might explain the opposite EDW signal found in the tropical Andes during daytime and nighttime. Changes in albedo are an ubiquitous driver for positive EDW in the Subtropics, for both the minimum and the maximum temperature. Changes in longwave radiation and humidity are also EDW drivers in the Subtropics but with different relevance throughout the seasons and during daytime and nighttime. Also, the western and eastern sides of the Cordillera might be influenced by different EDW drivers.


Introduction
Covering about 25% of the world's land surface and hosting about 12% of the world population, mountains provide a multitude of goods and services to both high-altitude environments and to downstream regions. Water is one of the most precious and critical resources delivered by mountains (e.g. Kohler and Maselli 2009;Price 2015), and about 40% of the world population depends on that water 1 3 for drinking purposes, agriculture, hydroelectricity. In the arid and semiarid regions of the Tropics and Subtropics, in particular, mountain regions play a very crucial role since they contain more than 80% of the freshwater made available to the countries around, especially when precipitation is scarce or absent (Viviroli et al. 2007(Viviroli et al. , 2011. For example in the Andean countries, particularly Bolivia and Peru, water supply during the dry season depends on glaciated basins located at high elevations (Vuille et al. 2008).
Mountain regions have been recognized to exhibit enhanced sensitivity to climate and environmental changes, and in particular to global warming, overall showing larger temperature trends ( 0.3 • C/decade ) compared to the globally-averaged trend ( 0.2 • C/decade , IPCC 2019, and references therein). The widespread glacier retreat observed in most mountain areas of the globe over the last decades is one of the most striking evidences of the temperature increase, though other changes, such as those in precipitation, are also at play (Bradley et al. 2006;Tennant et al. 2014).
The tropical Andes, hosting more than 99% of all the glaciers located in the Tropics (Kaser 1999), are a climate hotspot with increasing temperature trends both at the surface (Nogués-Bravo et al. 2007) and in the atmosphere (Bradley et al. 2006). The literature reports several studies which document the widespread retreat of Andean glaciers over the recent decades (Ceballos et al. 2006;Mark and McKenzie 2007;Seimon et al. 2007; Thompson et al. 2006;Vuille et al. 2008) in conjunction with the observed regional climate change (e.g. Vuille and Bradley 2000;Vuille et al. 2003). Past studies based on observations identified the increase in near-surface air temperature in the tropical Andes as the main driver of the glacier retreat occurred in the second half of the twentieth century. Also, warming was found to be markedly larger in the western side than in the eastern side of the chain (Vuille et al. 2003). According to the analysis performed by Vuille et al. (2008) over about 280 in-situ stations located between 1 • N and 23 • S , a near-surface air temperature increase of 0.68 • C was detected from 1939 to 2008, with a rate of 0.1 • C/decade . Changes in precipitation were found to be less notable and more regionally varying along the chain than the temperature changes (Vuille et al. 2008). Regional climate model (RCM) studies, focusing on projected temperature changes over South America, overall agree in showing a generalized warming trend all over the continent and identifying the Andean region as a climate hot-spot also by the end of the twenty-first century, with projected temperature increases larger than 4 • C (Solman 2013) compared to pre-industrial levels.
One of the main issues facing the research community working on mountain climate change is whether warming rates exhibit a dependence on the elevation and in particular if they are enhanced at higher elevations compared to lowerelevation counterparts. This would have several implications for the high-altitude cryosphere, mountain ecosystems, biodiversity and ultimately for the human societies which depend on the services delivered by mountains. Several global and regional studies report on elevation-dependent warming (EDW), i.e., whether warming rates are different in different altitude ranges (e.g., Palazzi et al. 2019;Pepin et al. 2015). The majority of them, analysing either observations or model simulations, point toward an amplification of warming rates with elevation (i.e. a positive EDW), but many uncertainties remain, including the strength of EDW, its dependence on the specific mountain region, and its drivers. On the one hand, these uncertainties are related to the difficulty of measuring EDW, because of the sparseness or even the lack of in-situ stations up to the highest elevations providing long-term temperature timeseries. On the other hand, uncertainties are inherent in numerical simulations with climate models with limited spatial resolution when applied in mountain environments and whose parameterizations of sub-grid processes are imperfect (Minder et al. 2018;Pepin et al. 2015). Most studies aiming at better understanding and disentangling the mechanisms leading to EDW have been based on climate models rather than on observations (Palazzi et al. 2017(Palazzi et al. , 2019Rangwala and Miller 2012) since the output of numerical models is spatially and temporally continuous and contains all the key variables that provide, at least in principle, a complete picture of the possible EDW drivers and of their mutual relationships. In studies analysing mid-latitude mountains, EDW has been frequently attributed to elevation-dependent changes in albedo (associated with, but not exclusively to, changes in snow cover), in downward longwave radiation and in specific humidity (e.g. Minder et al. 2018;Palazzi et al. 2017Palazzi et al. , 2019Rangwala et al. 2016). However, the majority of those studies also concurs in recognizing that other mechanisms are at play which are currently not (or not sufficiently well) implemented in the models (e.g. Pepin et al. 2015).
Compared to the studies focused on the northern hemisphere mid-latitude mountains, EDW in the Andes and especially in the tropical Andes has not been extensively studied. Recently, analysing satellite Land Surface Temperature (LST) trends in the period 2000-2017 in winter and their dependence on elevation in the Andean region between 7 and 20 • S , Aguilar-Lome et al. (2019) found different behaviours in daytime and nighttime LST trends as a function of elevation. While the daytime winter LST showed a strong positive trend with elevation, the same warming pattern was not found in the nighttime LST. Among the first studies analysing future climate changes in the tropical Andes ( 10 • N-27 • S ) is the paper by Urrutia and Vuille (2009) who used a regional climate model (RCM) based on two different emission scenarios (A2 and B2, belonging to the family of scenarios that preceded the Representative Concentration Pathways, RCPs, employed in the fifth IPCC assessment report) to show a significant warming in the area, enhanced at higher elevations and further amplified in the middle and upper troposphere. For the subtropical Andes (from 30 to 37 • S ), the model study by Zazulie et al. (2018) analysed future climate projections using five GCMs under the RCP 4.5 and RCP 8.5 greenhouse gas concentration scenarios (Riahi et al. 2011;Thomson et al. 2011), identifying an expected warming of 5 • C in winter at the highest elevations by 2075-2100 compared to the period 1980-2005. The authors also found a reduction in albedo by 20-60% at high elevations, which was identified as one key EDW driver.
The present paper focuses on future projections of warming rates in the southern hemisphere tropical and subtropical Andes, and in particular on the existence and driving mechanisms of EDW, highlighting possible differences as a function of the latitude, the season, the considered variable (either the minimum or the maximum temperature) and the side of the Andean Cordillera, i.e. considering its western and eastern slopes separately. Results for the end of the twenty-first century (2071-2100) are presented and compared to a historical model climatology (evaluated over the period 1976-2005) using simulations performed with one regional climate model (the RCA4 RCM, Samuelsson et al. 2011) nested into eight global model simulations whose projections are based on the RCP 8.5 greenhouse gas concentration scenario (Moss et al. 2010).
The paper is structured as follows: Sect. 2 introduces the area of study, the employed data and the methods used for their analysis; Sect. 3 presents the results on EDW, its existence and spatial variability, and considerations about its driving mechanisms; Sect. 4 discusses and concludes the paper.

Climatic characteristics of the study areas
The Andean mountain chain is a relatively narrow mountain range extending all along the western coast of South America, encompassing tropical and extra-tropical latitudes which exposes it to diverse climatic conditions. North of 30 • S , relatively cold and arid conditions prevail along the western slopes of the Andes, while the Andean eastern slopes are comparatively moister, warmer, and prone to higher precipitation (Garreaud 2009;Rangwala and Miller 2012;Solman 2013). This pattern reverses south of 30 • S where precipitation mainly occurs on the western side (owing to the influence of westerly winds) while the eastern slopes are characterized by semiarid conditions (Garreaud 2009). The seasonal shift of the inter tropical convergence zone (ITCZ) controls climate conditions over the tropical part of the Cordillera.
In terms of large-scale upper-level circulation, the tropical Andes between north of the equator and 15 • S are prone to mild easterly winds while the subtropical and extra-tropical portions of the chain are subject to westerly winds. In DJF (austral summer), the easterly winds extend beyond 15 • S and the Amazon deep convection favors the establishment of the upper level Bolivian High (Lenters and Cook 1997), one of the most prominent features of the summertime circulation in the Andes. The summertime subtropical westerly jet displaces southwards and becomes less intense. In JJA (austral winter) the easterly winds do not cross 10 • S latitude and the subtropical westerly jet has its maximum around 30 • S . Lowlevel circulation is rather complex and extremely important since it controls the transport of humidity and, thus, the precipitation patterns. Between about 10 • S and 35 • S low-level winds mostly flow from the south along the western Pacific coast and from the north along the eastern Andes slopes.
The present study considers the tropical and subtropical parts of the Andes Cordillera, between 0 and 37 • S latitude and 64 • W and 80 • W longitude. The domain study was divided into five sub-regions depicted as coloured boxes in Fig. 1 where they are sketched above two orography maps. One (panel a) is derived from the land elevation data provided by the NASA Shuttle Radar Topographic Mission (SRTM, Farr et al. 2007), while the other (panel b) is the orography used in the CORDEX regional climate model simulations employed in this study (see Sect. 2.2).

Model data
We considered an ensemble of model simulations obtained with the Rossby Centre Atmospheric regional model (RCA4, Samuelsson et al. 2011) nested into eight different global climate models -CanESM2, CSIRO-Mk3-6-0, EC-EARTH, CM5A-MR, MIROC5, HadGEM2-ES, MPI-ESM-LR, and GFDL-ESM2M from the Climate Model Intercomparison Project phase 5 (CMIP5) ensemble, whose spatial resolution and a key reference are summarised in Table 1. The RCA4 simulations are part of the Coordinated Regional Climate Downscaling Experiment (Giorgi et al. 2009)-South American domain (CORDEX-SAM) and have a spatial resolution of 0.44 degrees latitude-longitude (approx. 50 km). To date and to our knowledge, these simulations are the best available for South America in terms of spatial resolution and availability of different ensemble members. The results discussed in the paper refer to the 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 (Chou et al. 2014;Riahi et al. 2011). It is worth pointing out that all the analyses were repeated for the less extreme emission scenario RCP 4.5, pointing to a stabilization of the CO 2 concentrations by the end of the century, showing no differences in the spatial pattern of the seasonal temperature changes but just in the magnitude of the changes (which are strongest in the most extreme RCP 8.5 scenario). The projection period taken into account is 2071-2100 and it is compared to historical simulations in the period 1976-2005. Figure 1c shows the density of grid points in each 200 m-thick elevation bin for the CORDEX RCM simulations ( ∼ 50 km resolution) compared to the density distribution for high-resolution ( ∼ 90 m ) land elevation data, for the five study areas. With the exception of the highest elevations, the shape of the modeled elevation density distribution quite well reproduces the land elevation data distribution in Tp2, Tp3 and Sb1, and to a lesser extent in Sb2. The modeled elevation distribution in Tp1 is, instead, in less agreement with the observations, both in the representation of the lowest and highest elevations and in correctly simulating the shape of the distribution overall. As expected, the highest altitudes (above ∼ 3300 m for Tp1; above ∼ 4700 m for Tp2, Tp3 and Sb1; above ∼ 3700 m for Sb2) are not reproduced in the RCM grid. It is worth stressing that the smoothing of higher altitudes due to the model grid resolution could potentially compromise the assessment of EDW and its drivers, especially those related to changes in snow cover and related feedbacks. To assess the capability of the CORDEX  Farr et al. 2007) and b from CORDEX models ( 0.44 • longitude-latitude resolution). The tropical (Tp1, Tp2, Tp3) and subtropi-cal (Sb1, Sb2) study areas are delimited by coloured boxes. c Fraction of grid cells in each 200m-thick elevational bin across the five study areas for the SRTM (red) and the CORDEX (cyan) data, only grid cells with elevation above 500 m were considered  Delworth et al. (2006) models in reproducing a reasonable snow cover distribution, which is also coherent with the topography of the area, we evaluated the snow cover area in the models, limited to the multi-model mean, against the state-of-the-art Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data. We used MODIS/Terra Snow Cover Monthly L3 Global 0.05Deg CMG Version 6 data, providing monthly averaged snow cover in 0.05 degree (approx. 5 km) resolution Climate Modeling Grid (CMG) cells. Monthly averages are computed from daily snow cover observations in the MODIS/Terra Snow Cover Daily L3 Global 0.05Deg CMG (MOD10C1) dataset. MODIS data were downloaded from https:// nsidc. org/ data/ MOD10 CM/ versi ons/6. We performed the comparison on a seasonal basis by analysing the snow cover climatology over the period 2000-2020. Data from the CORDEX models have been averaged over that period through an extension of the historical runs ending in December 2005 using the RCP 4.5 scenario data. Results of this comparison are shown in Figs. S1-S3 of the Supplementary Information. Figures S1 and S2 show the climatological seasonal maps of the snow cover fraction in the period 2000-2020 for CORDEX and MODIS, respectively. The latter are in agreement with those shown in the study by Saavedra et al. (2018), who analysed the changes in Andes snow cover in the period 2000-2016 using the binary 8-days product from MODIS (MOD10A2). The only remarkable difference with Saavedra et al. (2018) are the high snow fraction values which we found around 20 • S latitude, which we did not found to be associated with the longer time period considered in our study. The seasonal snow fraction anomaly (CORDEX-MODIS) in Fig. S3 shows that the largest differences between the model mean and the satellite data are seen in the two subtropical areas (Sb1 and Sb2) along the Chile-Argentina border in JJA and SON, where CORDEX models overestimate snow cover from observations. Noteworthy, in the tropical areas (Tp1, Tp2 and Tp3) the measured and modelled data show a particularly good agreement. However, a well defined area is identified in Tp3 where the models underestimate the MODIS observations, in southwestern Bolivia near the border with Chile, around 20 • S . As mentioned above, this is related to the high snow fraction data observed in MODIS.

Methods
In order to assess and quantify EDW in the five study areas, the difference (or change) between the average in the period 2071-2100 (projection period) and the average in the period 1976-2005 (historical period) of the minimum (tasmin) and maximum (tasmax) surface air temperature was calculated and, then, fitted against elevation using, first, a linear regression model. This calculation was performed for each GCM-RCM pair of the model ensemble as well as for the multi-model mean. When not differently specified, the results of the present work mainly refer to the multi-model ensemble mean, the single GCM-RCM pairs being used to discuss the robustness of the results across the ensemble. A seasonal based analysis using the standard definition of seasons for the Southern Hemisphere has been performed: winter (June to August, JJA), autumn (March to May, MAM), summer (December to February, DJF), and spring (September to November, SON). Minimum and maximum temperatures were separately analyzed because they can exhibit different EDW signals owing to the mechanisms at play during nighttime and daytime. Due to the different climates of the western and eastern sides of the Cordillera, the two sides are analyzed separately. With eastern and western sides we refer to grid points respectively to the East and to the West of the highest grid point at each latitude. Similarly to other studies focused on different areas (e.g. Giorgi et al. 2014;Palazzi et al. 2019), in this study to reduce some of the influence of the coastal areas, only grid cells with elevation above 500 m a.s.l. were kept in the analysis. The slope of the linear regression between warming rates and elevation along with its statistical significance provides a measure of the sign and strength of EDW. A positive slope identifies a positive EDW (i.e., the warming rates increase with elevation), while a negative slope identifies situations in which warming rates decrease as the elevation increases (negative EDW). The statistical significance of the linear regression slope is assessed following Bendat and Piersol (2000), as briefly summarized in the following. Given the random variable w = 1 2 ln 1+ 1− , where is the theoretical covariance of the two variables (the temperature change and the elevation in our case), having a Gaussian distribution with mean = 1 2 ln 1+r 1−r and variance 2 = 1 N−3 , where r is the empirical correlation of the two variables and N is the number of observations, at a confidence level of 0.05 the null hypothesis of no correlation between the two variables can be rejected if (Bendat and Piersol 2000): The linear regression model allows a description of the general behaviour of the relationship between warming rates and the elevation. However, it is also known from previous studies that linearity can oversimplify the real EDW pattern, which might be better represented by more than a unique slope (e.g. Palazzi et al. 2019;Vuille et al. 2003;You et al. 2010). Therefore, using the same approach as in Palazzi et al. (2019) we investigated possible departures from linearity with a LOcal regrESSion (LOESS) method. Then, the altitudes where the EDW signal changes its slope were identified with a piecewise linear regression. The piecewise regression model estimates the slope shift in the temperature change series fitting two altitude segments across a breakpoint, z (Toms and Lesperance 2003), as follows: where the index i varies between 1 and the total number of points in the model grid, z i represents the vector of elevations and T i the vector of the temperature change (either minimum or maximum temperature change), i are the fit residuals (considered as independent and identically distributed random errors with zero mean and finite variance). z represents the breakpoint, i.e. the altitude at which the regression line changes the slope. 0 , 1 and 2 are the regression coefficients: 0 is the intercept, 1 is the slope value before the breakpoint and 1 + 2 is the slope after the breakpoint, so that 2 can be interpreted as the difference in the two slopes, before and after the breakpoint. The model parameters 0 , 1 , and 2 were estimated choosing the breakpoint, z , that minimizes the residual sum of squares (Chiu et al. 2006). This piecewise fitting procedure assumes one single breakpoint and forces continuity at it (Toms and Lesperance 2003).
In the literature, EDW has been linked to several mechanisms (Pepin et al. 2015): snow-albedo feedback, changes in vegetation and the tree line, changes in cloud cover and properties, processes associated with moisture and radiative fluxes, changes in atmospheric loading and properties of aerosol particles. To evaluate their possible contribution to EDW, following Rangwala and Miller (2012) and Rangwala et al. (2013Rangwala et al. ( , 2016 we considered as drivers the model variables whose changes mostly influence the energy balance at the surface, leading to temperature changes, i.e. albedo, surface downwelling longwave and shortwave radiation (rlds, rsds), and near-surface specific humidity (huss). For each season, the absolute change between the average in the projection period and the average in the historical period for each driver mentioned above ( albedo , rlds , rsds , and huss ) was calculated. Besides considering the absolute changes of rlds, rsds, and huss, we also took into account their relative changes ( rlds∕rlds 0 , rsds∕rsds 0 , and huss∕huss 0 ), i.e. the change of each variable relative to an average climatology, since previous studies have shown that they can be more effective, at least in the mid-latitude mountains, in determining elevation-dependent warming signals (Palazzi et al. 2017(Palazzi et al. , 2019Rangwala et al. 2013Rangwala et al. , 2016. In this study we focused on the tropical and subtropical Andes and both absolute and normalized changes were considered to assess the relative role of the EDW drivers and to highlight possible differences with the more explored mid-latitude mountain regions.
In order for the simulated albedo , rlds , rsds , huss , rlds∕rlds 0 , rsds∕rsds 0 , and huss∕huss 0 to be considered as evidence of EDW mechanisms, some conditions need to be satisfied. First, the drivers have to exhibit a dependence on the elevation that is physically consistent with the EDW sign. In other words, all the drivers but albedo must increase with height if EDW (in either the minimum temperature or the maximum temperature) is positive, and decrease with height if EDW is negative. Since albedo has an opposite effect on the temperature change (i.e. decreases in albedo lead to temperature increases), for clarity in presenting the results, − albedo was used in the analysis. Second, the drivers and the temperature changes must exhibit a positive spatial correlation even when they are altitudedetrended, i.e., when their dependence on the elevation is removed. The reader is referred to Palazzi et al. (2017Palazzi et al. ( , 2019 for further details on this methodology.

Results
Before discussing the specific results on EDW following the methodology described in Sect. 2.3, we present a comparison between the model climatology and an observation-based reference climatology over a common historical period. To this aim, we considered the recent ERA5 global reanalysis dataset (Hersbach et al. 2020) and used the reference period 1979-2005 to evaluate the historical model data against ERA5.
We first compared the annual cycle climatology of both the minimum and the maximum temperature and we show the results in Figs. 2 and 3, respectively. ERA5 data are shown in red while the model results are in black (multimodel mean of the GCM-RCM ensemble) and in grey (individual GCM-RCM pairs). The mean annual cycle is calculated over the period 1979-2005 in the five study areas considered in this study (columns in the figures) and for the East and West sides of the mountain range (rows). It is interesting to observe that the models (both individual realizations and their mean) exhibit a cold bias compared to ERA5 when the minimum temperature is taken into account (Fig. 2). Instead, the models either overestimate or underestimate ERA5 when the maximum temperature is analysed (Fig. 3), but the difference between model data and the reanalysis is considerably lower than for the minimum temperature. The study by Urrutia and Vuille (2009) also showed a cold bias exhibited by one RCM over the southern Andes (compared to CRU observations) related to the excess precipitation simulated by the model, which reduces incoming solar radiation and lowers air temperature, particularly over the eastern Andean slope. Similar findings were discussed in other studies (e.g. Fernandez et al. 2006) using regional climate models. Table S1 of the Supplementary Information shows the numerical values of the biases (CORDEX-ERA5), at the seasonal and annual basis, for the minimum and the maximum temperature and the two against the ERA5 climatology (red). The five analysed areas are shown in the columns, while the two sides of the Cordillera in the rows Fig. 3 The same as Fig. 2 but for the maximum temperature sides of the Cordillera. Tp3 and Sb1 show the largest minimum temperature bias, larger for the western than for the eastern side of the Cordillera; for the maximum temperature, Tp1 shows the largest CORDEX-ERA5 discrepancy (see Table S1). It is interesting to note that, for the minimum temperature, the magnitude of the bias is higher than the magnitude of the temperature perturbation in the RCP8.5 scenario (i.e., the change of the minimum temperature between the average in the period 2071-2100 and the average in the period 1976-2005), while it is lower for the maximum temperature. This comparison is shown in Table S2 of the Supplementary Information. Nevertheless and regardless of the temperature bias magnitude, the models reproduce roughly well the shape of temperature annual cycle, overall. For the minimum temperature in east facing Tp1 and Tp2, however, the amplitude of the annual cycle appears overestimated. This could indicate excess model sensitivity to climatic drivers which decrease during the austral winter, maybe humidity-driven downward longwave radiation. For the maximum temperature in east facing Tp1 and Tp2 the CORDEX simulations produce an annual temperature maximum in September and relatively warm conditions in August which are not seen in ERA5, which could again be indicative of an enhanced model sensitivity to a specific climate driver. The maximum temperature annual cycle for west facing Tp1 is not correctly simulated by the models compared to ERA5, as both the annual maxima and minima occur in the wrong points of the annual cycle. For completeness, a further evaluation of the model skills is presented in Figs. S4 and S5 (and in Table S3) of the Supplementary Information showing the time series (and the trends) of the minimum and maximum temperature in the models and in ERA5 for the five study areas and the two sides of the Cordillera. As for the model bias, similar considerations as for the climatological annual cycle apply in this case. In addition, Table S4 of the Supplementary Information compares the interannual variability of the minimum and maximum temperature time series in ERA5 and in the model ensemble. The interannual variability is calculated as the standard deviation of the detrended time series shown in Figures S4 and S5 of the Supplementary Information. For most of the seasons, sub-regions and sides of the Cordillera the interannual variability in ERA5 is within the range of variability displayed by the GCM-RCM pairs.
As stated in Sect. 2.3, the EDW assessment is based on the analysis of the temperature difference between two climates (and its relationship with the elevation). For this, the bias of the models against ERA5 has a negligible impact on the results, because we analysed the temperature differences between the end of the twenty-first century and the historical climatology, so any bias is removed.

Elevation-dependent warming assessment
Elevation-dependent warming is assessed by exploring the dependence of the minimum and of the maximum temperature change ( tasmin and tasmax respectively) on elevation for each area in the tropical and in the subtropical Andes and for each season. We recall that the change is calculated as the difference between the 2071-2100 climatology (evaluated under the RCP 8.5 scenario) and the 1976-2005 climatology. For completeness, the spatial maps of tasmin and tasmax in the five study areas obtained from the multi-model mean of the GCM-RCM ensemble are shown in Figures S6 and  S7 in the Supplementary Information. Both variables are expected to undergo a pronounced warming in the future with respect to the past; however, the maps highlight different EDW patterns between the minimum and the maximum temperature, further the warming is not homogeneous throughout the latitude range. For example, in the tropical areas tasmin is higher at lower elevations, at variance with the situation encountered in the subtropical areas (see Figure  S6). This and other characteristics of a possible differential warming with elevation, i.e. of the EDW, in the five study areas are discussed in detail in the following analysis. Figures 4 and 5 show the scatter plot between the elevation and, respectively, tasmin and tasmax . Two different colors are used to highlight the western (blue) and eastern (red) sides of the mountain range. For each of the five study areas and for each season three different regression lines are displayed to assess the dependence of the temperature changes on elevation (see Sect. 2.3 for details): a linear regression (dark blue and dark red for the eastern and western sides respectively), a LOESS fitting curve (light blue and orange for the eastern and western sides respectively) and a piecewise regression (lavender and red for the eastern and western sides respectively). The LOESS regression lines in Figs. 4 and 5 suggest that, in some cases, a nonlinear fitting able to capture the possible complex pattern of EDW may be more appropriate than a linear relationship. It can be noticed that in quite a few cases the LOESS (light blue and orange) and piecewise (lavender and red) regression lines overlap quite well throughout the whole elevation range, indicating that splitting the altitude domain into two segments and fitting with two separate lines could more properly describe the relationship between warming rates and elevation than a single linear regression line (dark blue and dark red). To better assess the robustness of a piecewise EDW behaviour across the entire model ensemble, Figures S8 and S9 of the Supplementary Information depict the piecewise regression lines for all individual GCM-RCM pairs (colored lines) and for their mean (black lines) for tasmin and tasmax , respectively. In those figures the two sides of the Cordillera are depicted with different line types (continuous for the western side and dashed for the eastern side) and they will be used in the discussion of the following EDW analysis. In particular, these plots highlight how each different GCM-RCM identifies the breakpoints and the slopes of the piecewise regression in the different seasons, areas and sides of the Andean Cordillera.
The minimum temperature in the Tropics (Tp1, Tp2, Tp3 in Fig. 4) shows no evidence of increased warming with elevation, i.e. of positive EDW, while the two subtropical areas show an opposite behaviour characterized by positive elevational gradients of warming rates, overall. These considerations apply to both eastern and western sides of the Andean Cordillera. Nonetheless, the two sides can exhibit different warming rates. This is particularly evident in Tp1 and in the two subtropical domains, in all seasons: in Tp1 the warming is larger in the eastern side (blue) while in the subtropics the western part exhibits amplified warming. This further justifies splitting the analysis in the two sides of the Cordillera. In Tp2 and Tp3 a smaller spread in the data is present and a different warming in the two sides of the Cordillera is not clearly detectable. Figure 4 shows that, for both the eastern and western sides of the chain, the piecewise regression (lavender and red lines) almost overlaps with the LOESS curve (light blue and orange lines) in all cases but one (Tp3 in JJA). The different slopes of the two fitting lines (piecewise regression line) mark a breakpoint, i.e. an elevation where the EDW behaviour changes. In most of the seasons, areas, and sides of the Cordillera, the change in the EDW slope is remarkable. The position of the breakpoint depends on the season, on the latitude and on the side of the Cordillera. In Fig. S8 of the Supplementary Information we compare the position of the breakpoint in all GCM-RCM pairs. Though in some cases the models do not agree with each other, the model ensemble is mostly coherent in providing a robust estimate of the breakpoint height.
Unlike for the minimum temperature, the maximum temperature change (Fig. 5) is overall positively correlated Fig. 4 Scatterplots of the seasonal minimum temperature changes versus altitude for the multi-model mean in the five study areas and in the RCP 8.5 scenario. The dark blue and dark red lines represent the linear regression model for the eastern and western sides of the Cordillera, respectively, the light blue (eastern side) and orange (western side) lines represent the LOESS fit and the piecewise regression is shown in lavender (eastern side) and red (western side). Point data referring to the eastern (western) side of the mountain range are colored in blue (red) with the elevation. A distinction between the EDW behaviour found in the western and eastern sides of the mountain range is evident. This is particularly emphasised in the tropical areas where, at lower altitudes, larger temperature changes in the eastern side (blue points) than in the western side (red points) are observed. Overall, the eastern side is characterized by a twofold EDW behavior (light blue piecewise regression line). In Tp1 a negative EDW at lower altitudes is followed by a positive one at higher altitudes. In Tp2 and Tp3 a slightly positive EDW at lower altitudes and a more positive one above are found in DJF, MAM and JJA; in SON EDW is negative throughout the whole elevation range. In Tp1, Tp2 and Tp3 the western side behaves similarly throughout the seasons, showing a linear increase of tasmax with elevation. In DJF, MAM, and JJA, the height where the EDW signal changes its slope in the eastern side, also marks the height where the relationship between tasmax and the elevation is the same in the eastern and western sides of the mountain range.
In SON, the eastern and western sides present similar tasmax values only at the highest elevations. In the Tropics, SON in the season with the most pronounced difference between the eastern and western warming.
In the two subtropical areas, the difference between the two sides of the mountain chain is less evident than for the three tropical areas but still existent. In Sb1 the relationship between tasmax and the elevation is overall well described by a linear regression, especially in the western side and with the exception of JJA. In Sb2, the warming pattern is reversed with respect to the Tropics, being characterised by an amplified warming in the western part. Here, a twofold behavior with a steeper slope at the lower altitudes is found in DJF, JJA and SON. Figure  S9 of the Supplementary Information, depicting the piecewise regression lines for the individual GCM-RCM pairs of the model ensemble, shows that where the EDW slope before and after the breakpoint markedly changes in the multi-model mean (see Fig. 5), the agreement between the Fig. 5 The same as Fig. 4, but for the maximum temperature change individual GCM-RCM pairs is high in both the slopes and the breakpoint elevation.
Though the relationship between warming rates and elevation can deviate from linearity as discussed above, in the literature the strength of EDW is commonly quantified considering the slope of the linear regression over the whole elevation range. Besides, using a unique slope, rather than two slopes and a breakpoint altitude, enables a straightforward comparison among subdomains, seasons and sides of the Cordillera. Figure 6 summarizes the linear regression slopes for tasmin and tasmax in all the seasons, for all the five areas and for all GCM-RCM model pairs (colored symbols) as well as the multi-model mean (black symbol) in the RCP 8.5 scenario. In this figure the different sides of the mountain range are also depicted (West is indicated with the square red symbol while East with a blue triangle). The solid line connects, for each side of the Cordillera, the multi-model mean point of the five areas.
The upper panels referring to the minimum temperature confirm the distinct EDW behaviour between the tropical and subtropical regions. In the western mountain side (red points) the EDW is negative in the three tropical areas and positive in the subtropics. In DJF and MAM, in particular, Sb1 behaves like a transition region between the negative EDW in Tp1, Tp2, and Tp3 and the positive EDW in Sb2. In the eastern side, a less uniform behaviour is found throughout the seasons though, overall, EDW is negative in the tropics (with the exception of Tp3 in DJF) and positive (or neutral in Sb1) in the subtropics with the exception of SON. The two subtropical areas overall exhibit a larger intermodel spread compared to the tropical areas.
For the maximum temperature (lower panels), the slopes are overall positive in all seasons, for all areas and almost all individual model members. The slopes are larger in the western side than in the eastern side, except for Sb1 and Sb2 in MAM and JJA. The EDW slope gradually increases with latitude in all cases except in JJA for the western side, in which EDW increases (decreases) from lower to higher latitudes in the tropics (subtropics), and in SON for the eastern side, where EDW is constant and almost neutral in the three tropical areas while positive in the subtropics.
The different EDW sign found for the minimum and the maximum temperature in the Tropics is a quite robust feature across the whole model ensemble. Remarkably, to our knowledge, the contrasting behaviour between nighttime and daytime EDW-with minimum temperatures displaying negative slopes and maximum temperatures displaying positive slopes-differs from the behaviour most commonly observed in other mountain areas in the mid-latitudes such as the Tibetan Plateau-Himalayas (e.g., Palazzi et al. 2017) where the minimum and the maximum temperatures display a coherent EDW signal (same sign, positive) and the Fig. 6 Elevational gradients of the seasonal temperature change in the tropical (Tp1, Tp2, Tp3) and subtropical (Sb1, Sb2) Andes, for each side of the mountain chain (blue, east and red, west) and for the multi-model mean (black symbols) in the RCP8.5 scenario. The minimum and maximum temperature changes are shown in the top and bottom panels, respectively. The filled colored symbols represent statistically significant slopes while the transparent ones are not significant (Eq. 1). The colored lines connect the multi-model mean values in each area minimum temperatures generally show larger warming and larger EDW than the maximum temperatures (Rangwala and Miller 2012). The different behavior found here, for nighttime and daytime EDW, corroborates the need to separately analyse the minimum and the maximum temperature in the Andes. It is interesting to note that, in a study focused on the tropical Andes based on the analysis of satellite LST data in the period 2000-2017, Aguilar-Lome et al. (2019) also found a different elevational gradient of LST trends in winter between daytime and nighttime conditions: daytime LST trends are amplified at higher elevations while nighttime LST trends exhibit a steady increase with altitude.

Analysis of the possible EDW drivers
To better understand the different EDW signals identified in the previous section, an assessment of the possible driving mechanisms is performed applying the methodology described in Sect. 2.3. The first step is to assess whether the elevational gradient of each relevant variable identified as possible EDW driver may in principle (i.e. based on physical considerations) explain the EDW sign. For this, Fig. 7 shows the slopes of the linear regression between each variable (albedo , rlds , rsds , huss , rlds∕rlds 0 , rsds∕rsds 0 , and huss∕huss 0 ) and the elevation in the RCP 8.5 scenario, for each season in the west (red) and east (blue) side of the Andean chain. The figure refers to the multimodel mean of the GCM-RCM ensemble only. For example, tasmin in SON (upper right panel of Fig. 6) exhibits a positive EDW in Sb2 on the western side of the Cordillera (red squares). Since the slopes ofalbedo , rlds , rlds∕rlds 0 and huss∕huss 0 , shown in Fig. 7, are positive (same EDW sign), these variables are possible EDW drivers while those with negative slopes have to be excluded owing to physical considerations. As another example, tasmax in SON exhibits a negative EDW in the Tp2 area in the eastern side of the Cordillera (blue triangles) that may be explained by rlds , since this is the only driver presenting a negative slope (Fig. 7). It is worth noting that,albedo in the Subtropics and the absolute and relative change in huss are the drivers showing the largest differences between the eastern and western side of the Cordillera.
As a second step, we calculated the correlation coefficient between the (minimum and maximum) temperature change and each of the variables that satisfied condition 1. Following Palazzi et al. (2017Palazzi et al. ( , 2019, all these variables were altitude-detrended and standardized (by dividing each change by its standard deviation over each of the considered domains) before calculating the spatial correlation with the temperature change. The aforementioned correlation coefficients are displayed in Figs. 8 and 9 for the minimum and the maximum temperature, respectively, and limited to the multi-model mean of the GCM-RCM ensemble. Analogous figures are shown in the Supplementary Information to present the results for all individual GCM-RCM pairs (Figs. S10-S13). White areas in the figures identify the cases in which either conditions 1 (a given possible EDW driver shows a dependence on the elevation which is physically consistent with the EDW sign) or 2 (the correlation between a given possible EDW driver and the temperature change is positive) are not satisfied. It is worth specifying that the correlations have been calculated after grouping and averaging the temperature change data and the drivers data into elevation bins, which makes the correlations stronger and more significant but not remarkably different from those that would have been obtained without data binning (not shown here). In the following, a discussion on EDW drivers in the tropical and subtropical Andes based on the correlations shown in Figs. 8 and 9 is presented. In order to assess the difference between the EDW mechanisms in the Tropics and in the Subtropics, major emphasis in the discussion is placed to the drivers which are persistent and robust in all three tropical areas and in all two subtropical areas.
As for the minimum temperature (Fig. 8), our analysis shows that in the tropics (Tp1, Tp2, and Tp3) EDW is overall not driven by the changes in albedo and in the incoming shortwave radiation, in both the west and east sides of the Cordillera (we recall that for the minimum temperature in the Tropics a negative EDW, i.e., decreasing warming rates with the elevation, was found). Only two exceptions are found, in Tp1 West MAM and Tp3 East SON where, however, the correlations are weak, lower than about 0.2. Instead, the (absolute and relative) change in the downward longwave radiation is a clear driver in all three tropical areas, whose role seems to become less relevant approaching the Subtropics (both rlds and rlds∕rlds 0 are strong drivers in Tp1, rlds∕rlds 0 does not play a role in Tp3). For the minimum temperature, the most relevant difference in the EDW drivers between the Tropics and the Subtropics is that albedo becomes a driver in the Subtropics, with stronger correlation in the west side of the Cordillera (in Sb1 in all seasons but JJA albedo is a driver only in the West). The changes in rlds and in huss, particularly the relative changes, are also important drivers in the Subtropics. The individual GCM-RCM pairs (Figs. S10 and S11 of the Supplementary Information) overall provide the same picture as their mean. For the majority of models, rlds , rlds∕rlds 0 , and huss are recognised as EDW drivers in the tropics (the signal becomes weaker moving southward); in the subtropics, also albedo and huss∕huss 0 play a role, particularly in Sb2. Tp1 is the region presenting not only the largest correlation between the drivers ( rlds and rlds∕rlds 0 ) and tasmin but also the largest agreement in the model ensemble, in both West and East sides of the Cordillera. Tp2 immediately follows, the correlation being slightly lower than for Tp1. In the subtropics, Sb2 is the region showing the highest correlation and agreement among the models, mainly in the western side.
As for tasmax , we recall that the EDW signal in the Tropics is opposite (and positive) with respect to the one found for tasmin . This also reflects in the identified EDW drivers. Changes in tasmax are mostly associated with changes in downward shortwave radiation (in Tp1 and Tp2) while for tasmin the main driver was the absolute and relative change in rlds. However, it is worth underlining that the correlations between the EDW drivers and the maximum temperature change is weaker than for tasmin (with a few exceptions) and the picture is less coherent across the seasons and the sides of the Cordillera. In Tp3, the EDW signal for tasmax is mainly associated with rlds∕rlds 0 and, to a lesser extent, huss∕huss 0 . In Sb1, the identified drivers of EDW and their role for the minimum and the maximum temperature are very similar in the two sides of the Cordillera. While albedo is the main driver in the west side, rlds Fig. 7 Elevational gradient for the EDW drivers, for the RCP8.5 scenario, in the five areas areas (columns) and four seasons (rows) on each side of the Andean Cordillera  Fig. 8 Correlation coefficient between each of the seven possible EDW drivers and the minimum temperature change, for the RCP8.5 scenario, in the five areas areas (columns) and four seasons (rows) on each side of the Andean Cordillera. The drivers are displayed along the x-axis. White boxes identify the cases in which the correlation is negative and/or the elevational dependence of a given driver has not the same sign as EDW (see text for details) Fig. 9 Correlation coefficient between each of the seven possible EDW drivers and the maximum temperature change, for the RCP8.5 scenario, in the five areas areas (columns) and four seasons (rows) on each side of the Andean Cordillera. The drivers are displayed along the x-axis. White boxes identify the cases in which the correlation is negative and/or the elevational dependence of a given driver has not the same sign as EDW (see text for details) and rlds∕rlds 0 are the main drivers in the East (except that in SON). In Sb2, albedo is the main driver in both the western side and the eastern side of the chain with the only exception of JJA where rlds and rlds∕rlds 0 show a significant correlation with the EDW signal. With the exception of huss∕huss 0 in SON for the western side, all the EDW drivers associated with tasmax are also associated with tasmin . This suggests that in the Subtropics the minimum and maximum temperature changes may be driven by the same drivers on each side of the Cordillera. The individual model pairs (Figs. S12 and S13 of the Supplementary Information) do agree with the identification of the changes in the downward shortwave radiation as the main drivers in Tp1 and Tp2 especially in the eastern side of the Cordillera. As for Tp3, while the main EDW driver for most individual models in the west side is rlds∕rlds 0 , coherently with the picture provided by the multi-model mean, in the eastern part the absolute and relative changes in rsds turn also out to be important. In the Subtropics, almost all individual models identify albedo as one main EDW driver in the West, together with the changes in rlds and huss∕huss 0 but only in Sb2 in JJA and SON. This partly applies to the eastern side of the Cordillera too, though the signal is more scattered and rsds and rsds∕rsds 0 also play a role in Sb1 (SON) and Sb2 (MAM).
To summarise, it might be useful to highlight the drivers presenting the highest correlations with the temperature changes for each area and side of the Cordillera: -in Tp1, the changes in the minimum temperature (in both eastern and western sides of the Cordillera) are driven by the changes in rlds while the changes in the maximum temperature (especially in the East) are driven by the changes in rsds. This might explain the opposite EDW signals found in this region during nighttime and daytime (see e.g. Fig. 6). Also observed in all the individual GCM-RCM pairs, this result is very robust; -the albedo change exists as an ubiquitous EDW driver in the Subtropics (in both west and east sides in Sb2, mainly in the west side in Sb1) where it always shows a high correlation with both the minimum temperature change and the maximum temperature change. In general, the model ensemble reproduces coherently this driver. In the Tropics, while for the minimum temperature the absence of albedo as an EDW driver is quite clear, its role in driving the maximum temperature changes is less certain; -the (absolute and relative) change in shortwave radiation is an evident EDW driver only for the maximum temperature in the Tropics; -when the Tropics and Subtropics exhibit the same EDW signal, i.e. for the maximum temperature, in general no common drivers are identified in the two regions. On the other hand, when the Tropics and Subtropics show an opposite EDW signal, i.e. for the minimum temperature, the change in longwave radiation is a common driver

Discussion and conclusions
The Andean Cordillera extends over different latitudes in South America being exposed to different climate regimes, which offers a unique possibility to explore possible differences in EDW along the chain and how the possible various EDW drivers may act. In this paper, the Cordillera was divided into five domains, three pertaining the Tropics and two the Subtropics. The study employed a small ensemble of eight regional climate model simulations belonging to the CORDEX-South American domain experiment together with their multimodel mean. For future projections, the RCP 8.5 IPCC emission scenario was considered. The analysis was also performed on the RCP 4.5 scenario, not shown in the paper, and no substantial differences were found.
The CORDEX models were first compared to MODIS and ERA5 data to assess their ability in reproducing snow cover and the temperature, respectively, over a historical period. As for the minimum temperature, the climatological annual cycle is well reproduced by the CORDEX models overall, though its amplitude is overestimated in east facing Tp1 and Tp2. A cold bias is observed in the model simulations compared to ERA5. As for the maximum temperature, the model bias is reduced particularly in the subtropical regions. Also, the shape of the annual cycle is overall well captured except in east facing Tp1 and Tp2 and in west facing Tp1. We found that CORDEX (the multi-model mean) overestimates the snow cover on the high Cordillera areas in the subtropics along the Chile-Argentina border, and underestimate snow cover around the latitude of 20 • S. Though the temperature comparison between CORDEX and ERA5 overall led to satisfactorily results, the snow cover biases could have implications for the interpretation of projected temperature changes through the albedo mechanism.
The EDW signal was assessed by means of linear, LOESS and piecewise regressions. The comparison among the three methods highlighted a clear bi-modal behaviour for most of the seasons in all the study areas and sides of the Cordillera. Remarkably, when using the piecewise regression, the regression lines before and after the breakpoint showed, for most GCM-RCM model pairs, very similar slopes. This suggested that describing the EDW with a bi-modal linear behaviour might be related to physical reasons, which were not analysed in this study but deserve future investigations. Here, we opted for a more standard EDW analysis, based on the use of a single slope to assess the EDW sign and strength, in order to make the analysis of the results for the Andes more easily comparable to those published in the current EDW literature for this and other regions.
In all the three tropical regions an opposite EDW signal during daytime and nighttime (i.e. in the maximum and in the minimum temperature) was identified. The minimum temperatures, in fact, showed a negative EDW (warming rates are lower at higher elevations) while the maximum temperatures exhibited a positive EDW (amplification of warming rates with elevation). The minimum temperature elevational gradient in the Tropics also seemed to be similar in the eastern and the western sides of the mountain range, while for the maximum temperature the elevational gradient was larger in the western side. Contextually, our analysis on the possible EDW drivers showed that, in the Tropics, no common driver was identified for EDW in the maximum temperature (mostly driven by changes in downward shortwave radiation) and EDW in the minimum temperature (driven by changes in downward longwave radiation and in specific humidity), which may support the observation of the contrasting EDW signal found in the two variables. This is particularly true for Tp1 and Tp2, while in the west side of Tp3, tasmin and tasmax exhibit weak common drivers in DJF and SON.
In the subtropical Andes EDW is positive in every season and for both the minimum and the maximum temperature (in DJF, MAM and SON in the eastern side, EDW in tasmin is neutral). One clear and important EDW driver in the Subtropics, common to both nighttime and daytime temperatures, is the change in albedo. Longwave radiation and humidity are also significantly correlated to EDW in the Subtropics, but with different relevance for the minimum and maximum temperature throughout the seasons and the sides of the Andes. These findings reflect the results on EDW drivers reported in previous studies focused on mountain areas in the northern hemisphere mid-latitudes.
The analysis on the possible EDW drivers highlighted some differences between the tropical and subtropical sectors of the Andean chain, suggesting that the Tropic somehow marks a boundary, with Tp3 representing a transition area between different EDW behaviours.
It is important to stress that other possible drivers not included in this study, such as mechanisms related to aerosol particles or clouds, might influence EDW. Further, the potential cross-correlations between the drivers, that in principle may dampen or amplify the effect of the drivers taken individually, are here neglected.
A better understanding of climate changes in mountain regions and of their possible elevational gradients is of utmost importance, given the influence of high-altitude regions on downstream livelihood, economies and societies. This work showed the peculiarity of the elevational gradients of warming rates in the Andean chain, highlighting that different EDW patterns and driving mechanisms are at play in the tropical and subtropical Andes as well as in the eastern and western sides of the mountain range. Though being a crucial water reservoir for the Andean countries, this mountain range is still under-represented in the EDW literature. Further observational and model studies are thus required to deepen our knowledge on EDW and its mechanisms and to account for new or currently not well represented processes at play in this high-elevation region.