Evaluation and projection of snowfall changes in High Mountain Asia based on NASA's NEX-GDDP high-resolution daily downscaled dataset

Of the landscape and surface atmospheric conditions used to estimate the precipitation phase, surface temperature is the one most commonly utilized. Due to the lack of a gridded observational snowfall dataset and the fact that NEX-GDDP does not provide raw snowfall as an output variable, we first need to establish the relationship between surface temperature and snowfall fraction, and then to calculate future snowfall. The empirical formula of surface temperature and snowfall fraction can be obtained from the observation station data of surface temperature, precipitation and precipitation type, based on the method proposed by Dai (2008). This approach is widely used for estimating the precipitation phase, as it requires fewer parameters and has acceptable accuracy (Jennings et al 2018). Further, this method determines precipitation types through a sigmoidal hyperbolic tangent curve to fit the observations of snow conditional frequency and surface temperature (equation (1)).

As a first step, we collected 120 602 daily temperature, precipitation, and type of precipitation samples from 77 stations in HMA from 1961–1979 (note that since 1980, the weather stations have not recorded precipitation type). Through these sample data, we calculated the snowfall frequency every 0.5 °C T_s bin from −10 °C to 10 °C, and then used the sigmoidal hyperbolic tangent curve to fit it. Finally, we applied the empirical formula to the temperature and precipitation data of NEX-GDDP to calculate the historical and future snowfall, as follows:

$$\begin{equation}SF = a\left[ {{\text{tanh}}\left( {b\left( {{T_s} - c} \right)} \right) - d} \right]\end{equation}$$

where SF and T_s are snow condition frequency (%) and surface temperature (°C), respectively, and a, b, c and d are dimensionless parameters obtained by using the Levenberg-Marquardt method for nonlinear least squares fitting, based on larger sample size observation stations in HMA. In this study, a = −49.02, b = 0.4321, c = 2.47, and d = 1. Figure S1 (available online at https://stacks.iop.org/ERL/15/104040/mmedia) shows that the SF–Ts relationship for snow and rain can be well-fitted with a hyperbolic tangent function in HMA, giving the coefficient of determination R² as 0.9994 (p < 0.01, figure S1(a)). In order to further verify the reliability of this empirical formula, we also compared the annual snowfall of 77 stations from 1961–1979, based on observation and formula calculations. Generally speaking, the annual snowfall calculated based on the formula has a certain underestimation error, but the determinant coefficient R² also reaches 0.96 (p < 0.01) (figure S1(b)).

Four snowfall indices are evaluated in the present work: (a) mean snowfall amount (Smean), (b) snowfall days (Sd), (c) snowfall fraction (Sf), and (d) maximum 1-day snowfall (Smax). These indices, along with their units and abbreviations, are listed in table S2. Prior to performing residual analysis, daily snowfall ≤1 mm is set to zero for both observations and simulations. This is because climate models often suffer from a high incidence of drizzle, so it is difficult to measure small amounts of precipitation (Frei et al 2018).

To quantitatively assess the effects of temperature and precipitation and their interactions with snowfall variations, a more theoretical framework proposed by Krasting et al (2013) is adopted in this study. Snowfall ($S$) is expressed as the product of snowfall fraction ( $f$ ) and precipitation ($P$):

$$\begin{equation}S = P \times f\end{equation}$$

Changes in precipitation ${{ \Delta }}P$, snowfall ${{\Delta }}S$ and snow fraction ${{\Delta }}f$ can be written as:

$$\begin{equation}{{\Delta }}P = {P_1}-{P_0}\end{equation}$$

$$\begin{equation}{{\Delta }}S = {S_1}-{S_0}\end{equation}$$

$$\begin{equation}{{\Delta }}f = {f_1}-{f_0}\end{equation}$$

Through some algebraic manipulation, ${{\Delta }}S$ change may be rewritten as:

$$\begin{equation}{{\Delta }}S = {S_1} - {S_0} = {P_0}{{\Delta }}f + {f_0}{{\Delta }}P + {{\Delta }}P{{\Delta }}f{\text{ }}\end{equation}$$

Where ${P_1}$, ${S_1}$, and ${f_1}$ represent annual precipitation, snowfall, and snowfall fraction during 2070–2099, respectively, and ${P_0}$, ${S_0}$, and ${f_0}$ are the average values during 1976–2005, respectively. ${P_0}{{\Delta }}f,{\text{ }}{f_0}{{\Delta }}P,{{ \Delta }}P{{\Delta }}f$ can be regarded as the contribution of temperature and precipitation changes and their interaction with the overall change of annual snowfall.

We also investigated changes in the number of days of different snowfall events. According to the classification standard of China Meteorological Administration, 1–2.5 mm day⁻¹ is light snowfall, 2.5–5 mm day⁻¹ is moderate snowfall, 5–10 mm day⁻¹ is large snowfall, and more than 10 mm day⁻¹ is heavy snowfall (Zhou et al 2018).

Consistent with Elsner et al's (2010) and Littell et al's (2018) classification criteria, and to determine the future changes of snow-dominated areas, we further defined areas with Sf ≤ 0.1 as rain-dominated, 0.1 < Sf < 0.4 as rain-snow transitional, and Sf ≥ 0.4 as snow-dominated.

To assess the ability of models to project snowfall, the pattern correlation (PCOR) and model anomaly bias (BIAS) are employed. BIAS measures the average tendency of the simulated data with observations, while PCOR describes the temporal and spatial correlation between the simulated data and observations.

For the projection, we focus on the ensemble results at the end of the 21st century. The ensemble mean is simply taken as an arithmetic average of the 21 simulations. To provide a range of future snow scenarios for HMA, we calculated absolute and relative changes in snowfall for the years 2070–2099, relative to a baseline period of 1976–2005 under RCP4.5 and RCP8.5 emissions scenarios.

To assess the ability of models to project snowfall and spatial patterns, a comparison of snowfall indices from the NEX-GDDP dataset for the period 1976–2005 is performed against gridded observations for HMA (figures 2 and S2). Various atmospheric circulation system (East Asian monsoons, Southwest monsoon and westerly circulation), together with huge terrain and horizontal gradient effects, exert a complex control effect on the spatial pattern of snowfall in HMA (Deng et al 2017; Yao, 2012). High values for observed amounts and days of snowfall in HMA were mainly distributed in western Tienshan Mountains and southwest and southeast of the Tibetan Plateau. This spatial pattern is successfully captured by the multi-model mean ensemble (MME) results from NEX-GDDP, with spatial correlation coefficients of 0.98 (p < 0.01) for Smean and 0.94 (p < 0.01) for Sd in HMA. Similarly, the southeast-northwest gradient of snowfall fraction is also well-captured by the MME, with high spatial correlations of 0.96 (p < 0.01).

Zoom In Zoom Out Reset image size

Although NEX-GDDP can accurately simulate spatial patterns for three of the four snowfall indices adopted in this study, the underestimated snowfall bias is quite obvious between NEX-GDDP and observations for some individual models. The MME results could, however, reduce this bias and provide a better performance than the individual model, where the mean bias for Smean, Sd, and Sf is −4.32%, −6.61% and −2.21%, respectively. On the one hand, this underestimation error is due to the underestimation of simulated snowfall by the method proposed by Dai (2008). On the other hand, previous studies have confirmed that the precipitation data of NEX-GDDP has a negative deviation from the observed data (Bao and Wen 2017). In terms of Smax, NEX-GDDP shows a low performance in simulating the spatial distribution of extreme snowfall, especially in Pamir, Hindu Kush, and the Himalayas, where the maximum snowfall is substantially overestimated. The mean bias for Smax across the entire HMA region is 15.05%. Based on the above analysis, only three snowfall indices (excluding Smax) are used to project future variability of rainfall events in HMA.

Snowfall is sensitive to climate change and is controlled by both temperature and precipitation. Relative to 1976–2005, the annual mean temperature for the end of the 21st century is simulated to increase by about 3.1 °C and 5.4 °C under RCP4.5 and RCP8.5, respectively (figure S3). Overall, the temperature shows a larger increase in eastern HMA than in western HMA. An increase in total precipitation is also simulated (53.9 mm for RCP4.5 and 131.4 mm for RCP8.5; figure S3). Precipitation changes are spatially heterogeneous and show a decreasing gradient from southeast to northwest.

Against this background, figure 3 shows ensemble-projected percentage changes in the amount of total snowfall in 2070–2099 relative to 1976–2005 for both emissions scenarios. Despite the increase in total precipitation, snowfall diminishes in most parts of HMA for both. The MME decreases in mean annual snowfall by the end of the century amount to −18.9% and −32.8% for RCP4.5 and RCP8.5, respectively (table S3). Losses >50% projected by the models are mainly concentrated in East Tienshan, East Kun Lun, Qilian, South and Eastern Tibet, and Hengduan. The largest reduction in snowfall is projected for Qilian due to its low topography, with values of −85.3% under RCP4.5 and −93.1% under RCP8.5. In contrast, increases in snowfall are projected to occur under RCP4.5 in some high-elevation regions such as Hissar Alay (2.5%), Pamir (17.1%), and Karakoram (6.0%) (Table S3).

Zoom In Zoom Out Reset image size

In general, higher temperatures usually result in less snowfall. However, if there is no temperature change, greater amounts of precipitation will result in more snowfall. To further test this hypothesis and quantitatively evaluate the role of precipitation and temperature in snowfall change, the change in snowfall is decomposed into three terms (${P_0}{{\Delta }}f,{\text{ }}{f_0}{{\Delta }}P,{{ \Delta }}P{{\Delta }}f$) following equations (2)–(6) in section 4.2. Figure S4 confirms that the projected changes in snow fraction contribute to less snowfall by the end of this century in almost all areas. Specifically, our results indicate that warmer temperatures result in snowfall decreases of 29.1% and 51.4% under RCP4.5 and RCP8.5, respectively.

On the other hand, the effect of precipitation is opposite to that of temperature. Increased precipitation causes higher snowfall amounts without any changes in temperature. The models show that the average projected increases in precipitation by the end of the 21st century can result in an increase in snowfall of about 8.6% and 27.2% under RCP4.5 and RCP8.5, respectively. Moreover, the interaction term between temperature and precipitation, which typically plays a negative role in snowfall, is small and does not have a first-order effect in determining the overall sign of snowfall. In HMA, most sub-regions are dominated by ${P_0}{{\Delta }}f$, which explains the decrease in snowfall. However, in regions like Pamir and Hissar Alay, the effects of${\text{ }}{P_0}{{\Delta }}f$ on snowfall trends are offset by the effects of ${f_0}{{\Delta }}P$, indicating an increase in snow.

Decreases in snowfall frequency are very similar to changes in mean snowfall. By the end of the current century, the frequency of snowfall is projected to drop by 29.6% under the RCP4.5 scenarios and 47.3% under the RCP8.5 scenarios (figure S5). Furthermore, there is an obvious difference in the decreasing range of snowfall frequency between eastern HMA and western HMA. In most of eastern HMA, snowfall frequency experiences a 76.8% to 94.4% reduction. Qilian experiences the greatest change, with a decrease of 86.7% under RCP4.5 and 94.4% under RCP8.5. Most regions in eastern HMA show at least moderate change and some regions like Pamir and Karakoram even exhibit an overall increase in the snowfall frequency under RCP4.5.

To examine the regional differences of changes in snowfall days, the projected changes in light snowfall, moderate snowfall, large snowfall, and heavy snowfall events across two emissions scenarios were investigated (figure S6). In HMA, the frequency of light and moderate snowfall accounts for 73.2% of the total snowfall days. Thus, the changes in these two types of snowfall events can well explain the changes in snowfall frequency across the whole HMA. This is confirmed by the coincidence of the spatial structure of the changes in light snow and moderate snowfall with the spatial structure of total snowfall days. In other words, the projected increase in light and moderate events leads to an increase in snowfall days in western HMA, while declines in frequency in light snowfall and moderate events yield a negative contribution to the total snowfall days in eastern HMA. Interestingly, an increase in both large and heavy snowfall events was observed in eastern HMA but an overall decrease was observed in western HMA, except for Hissar Alay and Pamir. It should be pointed out that although the number of days for large and heavy snowfall accounts for only a small proportion of total snowfall days, it nonetheless has a significant impact on the economy and human life (Changnon and Changnon 2006).

A warming climate will cause more precipitation to fall as rain, resulting in a decline in the snowfall fraction. In HMA, relative to 1976–2005, MME decreases in the mean snowfall fraction for 2070–2099 are predicted to be −26.7% and −42.3% under RCP4.5 and RCP8.5 scenarios, respectively. Corresponding to the decrease in snowfall in East Tienshan, East Kun Lun, Qilian, Inner Tibet, South and Eastern Tibet and Hengduan, there is also a robust decline in snowfall fraction in those areas. Taking Qilian as an example, the snowfall fraction dropped by 83.0% and 91.8% under RCP4.5 and RCP8.5, respectively. Benefiting from the increase in snowfall, the snowfall fraction showed a slight decrease in Hissar Alay, Pamir, West Tienshan, Hindu kush, Karakoram, and West Kun Lun, ranging from −22.0% to 3.0% under RCP4.5 (figure S7).

To better understand the impact of future warming on snowfall rates, we further defined areas as rain-dominated, transitional, and snow-dominated (figure 4), based on classification criteria in Littell et al (2018). Under historical records, regions characterized as either rain-dominated, transitional or snow-dominated account for about 19.0%, 43.4% and 37.6% of the entire HMA, respectively. Snow-dominated areas occupied most of northwest HMA, while rain-dominated areas are mainly located in the southeastern margin of HMA and the Qaidam Basin.

Zoom In Zoom Out Reset image size

By the end of the century, snow-dominated areas are projected to drop to about 23.9% and 17.9%, whereas rain-dominated regions are expected expand to 42.0% and 53.9% under RCP4.5 and RCP8.5, respectively. Meanwhile, the transitional zones between rain and snow are shrinking and disappearing, indicating that large areas which are currently snow-dominated will become rain-dominated. In the future, the distribution of snow-dominated and precipitation-dominated areas will have distinct spatial distribution characteristics. Snow-dominant areas can only be observed in Pamir, Western Kunlun Mountains and the Karakoram, while rain-dominant areas occupy nearly all of eastern HMA.

By the end of the present century, warmer temperatures in HMA resulting from anthropogenic climate change will reduce the fraction of precipitation that falls as snow, shorten snowfall duration, and decrease snowfall amounts. These projected changes are consistent with observed and simulated changes in snowfall in most parts of the world, including the western United States of America (USA), the Alps, Switzerland, and the Arctic (Feng and Hu 2007, Serquet et al 2011, Krasting et al 2013, Bintanja and Andry 2017, Frei et al 2018). For example, using daily outputs of temperature and precipitation downscaled for 20 CMIP5 models, Klos et al (2014) further pointed out that by mid-21st century, western USA is projected to see a 30% reduction in wintertime snow-dominated area. Within a similar time period, the snow-rain transitional zone in central and eastern USA will shift to the north at magnitudes of 2 °C–6 °C under future warming (Ning and Bradley 2015). Driven by the projected warming, snowfall in most parts of the Alps will also decrease by the end of the present century. Furthermore, with the expansion of areas dominated by rain, those dominated by snow are shrinking in HMA, as has also been reported in the Arctic. Bintanja and Andry (2017) found that precipitation will be predominantly rainfall by 2091–2100, and that only Greenland will at that time experience snowfall fractions over 80%.

Our results are comparable to those in the same region. For example, Zhou et al (2018) forecasted snow cover change in the Tibetan plateau, based on a regional climate model. Zhou et al (2018) indicated that, by the end of the 21st century, snow days will decrease and mean intensity will increase by 17.7% and 3.3%, respectively, in the Eastern Tibetan plateau. Further, the ensemble-projected decrease in light snowfall amount is 117.9%, while the decrease in the number of snowfall days is 18.4%. The authors also pointed out that the main reason for the decrease in snowfall days is the decrease in light snowfall events, and that the decline in snow days and snow water equivalent in the Tibetan Plateau is more pronounced than in other regions in China (Zhou et al 2018). Regional research on HMA also indicates that snowfall under a strong anthropogenic forcing scenario in the Himalayas would likely decrease by 30%–50% in the Indus River Basin, 50%–60% in the Ganges River Basin, and 50%–70% in the Brahmaputra Basin by 2071–2100 (Viste and Sorteberg 2015). However, it should be pointed out that snowfall is less sensitive to warming in the Karakoram than in the Himalayas, due to a unique seasonal cycle (Kapnick et al 2014), and that this phenomenon was also found in our study.

At the same time, some regional studies in HMA have begun to pay attention to the impact of these changes on snow metrics, which show a reduction in snow water equivalent, springtime snow cover, and snow depth (Chen et al 2016, Smith and Bookhagen 2018, Li et al 2019, Tang et al 2019). As the temperature continues to rise, the number of days <0 °C is expected to decrease, leading to an increase in rain frequency. Subsequently, total snowfall is expected to decrease as a result of a lower probability for precipitation to fall as snowfall. Moreover, simulated HMA warming exhibits a large seasonal cycle, with warming peaking in late autumn and winter, which is typically the snowfall season. Unlike monthly increases in rainfall, snow decreases in most months, especially in spring and early autumn, which causes the snowfall fraction to substantially diminish in the April, May, November and December (figure S8).

The projected snowfall showed a more robust decline in eastern HMA than in western HMA. This difference can be explained by climatology, orographic influence, and unique seasonal recycling.

Climatology. The annual mean temperature for our baseline period (1976–2005) in HMA indicated that annual mean temperature was much higher in eastern HMA compared to western HMA. For example, in eastern HMA (e.g. Qilian, Hengduan, and South and East Tibet), baseline annual mean temperatures exceeded 6 °C, while in Kun Lun, which is located in western HMA, the reported baseline annual mean temperatures were below −6 °C. A detailed analysis comparing baseline annual mean temperature and snowfall indicated that the projected declines in annual snowfall were most acute in climatologically warmer regions (figure S9). So, for example, from −6 °C to 6 °C, there was a snowfall decline of −41.2% to −54.9% for RCP4.5 and −55.0% to −74.1% for RCP8.5. Another unfavorable factor restraining snowfall increases in eastern HMA is that the projected temperature rise in this area is much faster than the rise in western HMA, which can be clearly seen in the temperature projection changes in figure S3. For example, under RCP4.5, the temperature in western Tienshan increased by 2.2 °C, while that in the Qilian Mountains rose by 3.4 °C. This means that in the future, the annual average temperature in most eastern HMA regions will be higher than the freezing point, preventing the increased precipitation from falling as snow.

Elevation Influence. Changes in snowfall are also projected to be strongly elevation-dependent. Overall, the altitude of western HMA is higher than that of eastern HMA. Numerous studies have shown that large decreases have been projected at lower elevations, with relatively little change or even increases at the highest elevations (Brown and Mote 2009, Wi et al 2012, Lute et al 2015, Deng et al 2017, Frei et al 2018). The sensitivity analysis of snowfall to altitude shows that although snowfall in low-lying areas in HMA could be reduced by more than −80% for RCP8.5, high-elevation regions could experience slight increases (figure S9). The main reason for the increased snowfall in high altitude areas is that the temperature is still below freezing point in the future, and also higher temperatures suggest greater atmospheric moisture entrainment.

Unique seasonal recycle. The unique seasonal circulation system results in a difference in snowfall change between eastern and western HMA. Eastern HMA is mainly affected by a monsoon climate, while western HMA is mainly affected by westerly circulation. We can take the Karakoram Mountains in western HMA as an example. Under the effect of low temperatures in winter, the westerly circulation provides ample water vapor for this region, which protects the Karakoram range from snowfall losses. However, this situation differs significantly from that in eastern HMA, where precipitation is affected by a monsoon climate that is mainly concentrated in the warm season and thus not conducive to snowfall increase. Previous studies confirm that glaciers in the Karakorum region have shown mass stability and even expansion (Bolch et al 2012), while the nearby Himalayas and Tibetan Plateau have observed glacial mass loss (Kang et al 2010). Simultaneously, a recent study conducted by Kapnick et al (2014) also found a much weaker sensitivity of snowfall to warming in the Karakoram than in the Himalayas due to a unique seasonal cycle (dominated by non-monsoonal winter precipitation).

The hydrological impacts of these changes on snowfall indicators have so far been decreased spring snow and snowpack, flooding, and early snowmelt (Barnett et al 2005, Berghuijs et al 2014, Chen et al 2018, Li et al 2019). In HMA, more than one billion people depend directly or indirectly on water sourced from melting snow (Smith et al 2017). As a buffer, snow releases meltwater constantly over a longer period of time. Therefore, if there is more rainfall than snowfall in winter and early snowmelt in spring, the risks of winter and spring flooding and summer drought are elevated, thus further aggravating the contradiction between supply and demand of water in these areas.

In addition, if the current water availability level is to be maintained, any reduction in the natural storage of snow will increase the demand for reservoir capacity (Barnett et al 2005). Unless there is additional storage infrastructure, many parts of the region will face difficult choices between conflicting objectives such as hydropower, irrigation, navigation, and ecosystem support. This situation may amplify as snowfall continues to decrease, snowmelt occurs earlier in the year, and water demand increases.

Moreover, decreasing snowfall and snow cover in HMA could also influence the Asian summer monsoon and rainfall. Wu and Qian (2003) pointed out that, in Eastern and Southeastern Tibet, heavy snow years lead to a weakened Asian summer monsoon, with summer precipitation in South and Southeast Asia being less than normal. Therefore, understanding the changes in snowfall in HMA has important implications for water resources management far beyond the HMA region.

This research was supported by the National Natural Science Foundation of China (41971141, U1903208, 41630859). The authors gratefully acknowledge the Youth Innovation Promotion Association of the Chinese Academy of Sciences (No. 2018480). Data from the observation stations were obtained from the China Meteorological Data Sharing Service (http://www.cma.gov.cn/). NEX-GDDP can be downloaded from its official website (https://dataserver.nccs.nasa.gov/thredds/catalog/bypass/NEX-GDDP/bcsd/catalog.html). The gridded temperature and precipitation data came from a land surface hydrological model simulation (http://hydrology.princeton.edu/data/pgf/v1/0.25deg/daily/).