Snowfall increase counters glacier demise in Kunlun Shan and Karakoram

Glaciers in High Mountain Asia provide an important water resource for communities downstream and they are markedly impacted by global warming, yet there is a lack in understanding of the observed glacier mass balances and their spatial variability. In particular, the glaciers in the western Kunlun Shan and Karakoram ranges (WKSK) show neutral to 10 positive mass balances despite global warming. Using models of the regional climate and glacier mass balance, we reproduce the observed patterns of glacier mass balance in High Mountain Asia of the last decades. We show that low temperature sensitivities of glaciers and an increase in snowfall, for a large part caused by increases in evapotranspiration from irrigated agriculture, result in positive mass balances in WKSK. The pattern of mass balances in High Mountain Asia can thus be understood from the combination of changes in climatic forcing and glacier properties, with an important role for irrigated 15 agriculture.


Introduction
Glaciers in High Mountain Asia (HMA, see Fig. 1) show a very diverse response to the changing climate. Most glaciers are losing mass, as expected, but surprisingly, glaciers are stable or growing in a region northwest of the Tibetan plateau, a phenomenon dubbed the Karakoram anomaly (Hewitt, 2005). Recent studies have mapped glacier mass losses and velocity 20 changes in detail and have shown that the regions of largest glacier growth and acceleration are the Kunlun Shan and parts of the Pamir, with the glaciers in Karakoram being close to a steady state (Brun, Berthier, Wagnon, Kääb, & Treichler, 2017;Dehecq et al., 2019;J. Gardelle, Berthier, Arnaud, & Kääb, 2013;Julie Gardelle, Berthier, & Arnaud, 2012;Kääb, Treichler, Nuth, & Berthier, 2015;Lin, Li, Cuo, Hooper, & Ye, 2017). Part of the mass balance variability seems correlated to differences in the temperature sensitivity, i.e. the change of mass balance for a certain change in temperature, of the glaciers (Sakai & 25 Fujita, 2017), but this alone cannot explain why some glaciers are actually growing. https://doi.org/10.5194/tc-2019-228 Preprint. Discussion started: 4 November 2019 c Author(s) 2019. CC BY 4.0 License.

Regional climate model 45
Irrigation can influence the regional and global climate (Cook, Shukla, Puma, & Nazarenko, 2015;Lee, Sacks, Chase, & Foley, 2011;Lobell, Bonfils, & Faurès, 2008;Puma & Cook, 2010;Sacks, Cook, Buenning, Levis, & Helkowski, 2009). Since the regions surrounding HMA host the largest irrigated areas in the world, e.g. the Indo-Gangetic Plain (see Fig. 1), irrigation potentially influences the regional climate in HMA (Cai, Luo, He, Zhang, & Termonia, 2019;de Kok et al., 2018). However, available re-analysis datasets do not fully include irrigation, and generally have a relatively coarse spatial resolution. Hence, 50 we downscaled ERA-Interim (Dee et al., 2011) re-analysis data using the Weather Research and Forecasting model (WRF, Skamarock & Klemp, 2008) to obtain a high-resolution climate dataset between 1980-2010 for High Mountain Asia. We artificially applied irrigation to the surface by adding a precipitation term each time step, with a rate that is determined per month. The amount of irrigation applied per grid cell was based on the monthly water demand, which indicates the amount of irrigation needed to compensate evapotranspiration, that was calculated by PCR-GLOBWB (van Beek & Bierkens, 2008;van 55 Beek, Van der Esch et al., 2017;Y. Wada, Wisser, & Bierkens, 2014;Yoshihide Wada et al., 2011).
In this way, the irrigation amount is not easily overestimated, as could be the case when e.g. soil moisture would be kept constant. In reality, insufficient water might be available to meet the predicted demand, whereas inefficient irrigation will result in a larger water gift than predicted.

60
We used WRF, version 3.8.1, to downscale ERA-Interim data to a resolution of 20x20 km, with 50 vertical levels. Settings are nearly identical to our previous work (de Kok et al., 2018), and are shown in Table 1. Additionally, we now use grid nudging of the upper 35 vertical levels for horizontal wind, temperature, and humidity, as opposed to only forcing the model at the boundary in our previous study. The values of the nudging parameters have been found to perform well in similar studies (Collier & Immerzeel, 2015;Otte, Nolte, Otte, & Bowden, 2012), and are: 0.001, 0.001, and 0.00005 s -1 for horizontal winds, 65 temperature and water vapour, respectively.  (Kain, 2004) Planetary boundary layer YSU scheme (Hong, Noh, & Dudhia, 2006) (Beljaars, 1995;Dyer & Hicks, 1970;Paulson, 1970;Webb, 1970;Zhang & Anthes, 1982) Land surface Noah-MP (Niu et al., 2011) (Dlugokencky, Lang, Mund, Crotwell, & Thoning, 2018) and AGAGE (Prinn et al., 2000) data, as aggregated by the European Environment Agency (www.eea.europa.eu, accessed March 2018), and are kept constant throughout each year. We used a 10-day spin-up for each month and run each month separately to be able to include a different irrigation amount each month. Monthly initialisation of the snow cover, surface and soil temperature, and surface moisture was taken from GLDAS 2.0 (Rodell et al., 2004) monthly mean values. We 75 checked convergence between months for temperature and precipitation and they agreed within a few percent for all selected points.

Glacier model
To assess the response of the glaciers to the atmospheric forcing, we employ a glacier mass balance gradient model (Kraaijenbrink, Bierkens, Lutz, & Immerzeel, 2017). The model assumes a calibrated mass balance gradient along the glacier, 80 and includes parameterisations for the effects of supraglacial debris on the surface mass balance and the downslope mass flux.
We modelled all individual glaciers in HMA larger than 0.4 km 2 (n=33,587) transiently for the period 1980-2010 (Kraaijenbrink et al., 2017). For ease of comparison with published observations, we select only those larger than 2 km 2 for the final analysis, which represent 95% of the glacier volume in HMA. Initial mass balance conditions in 1980 were set to be stable, using ERA-Interim data to locally calibrate the mass balance gradient of each glacier, while all other initial and reference 85 conditions as described in the original study (Kraaijenbrink et al., 2017) were maintained. To modulate the mass balance gradient of the glacier over time, we applied annual precipitation changes derived from annual changes in WRF snowfall and temperature changes determined from annual changes in WRF melt season temperatures, i.e. when average daily temperature is above -5 °C. A threshold value of 0 °C did not significantly change the glacier mass balance results for most of HMA, but meant that temperature changes for the coldest points could not be determined, since they are always below 0 °C. We did not 90 take into account whether the WRF grid point was glacierised or not. To reduce potential biases imposed by spurious reference conditions, the reference for both deltas was taken to be the average of the first six modelling years, i.e. 1980-1985. We performed three separate glacier model runs to evaluate the attribution of snow and temperature to the glacier mass changes in and precipitation sensitivity of the glaciers, we also performed calculations using a fixed temperature or snowfall trend for all 95 of HMA, with the other variable kept constant.

Moisture tracking
Our moisture tracking model (Tuinenburg, Hutjes, & Kabat, 2012) follows the moisture associated with precipitation backwards in time to determine where the moisture first enters the atmosphere. It therefore establishes a direct causal link 100 between evapotranspiration and precipitation downwind. For the moisture tracking, we clustered locations that have similar climates in terms of seasonality, since these will likely also have similar moisture sources. For the clustering, we used the monthly mean values of precipitation, horizontal wind fields at 400 hPa, and 2m-temperatures, with means subtracted and divided by the standard deviations, to perform a k-means clustering using 13 clusters. We performed the clustering with different numbers of clusters and found 13 to give reasonably-sized, yet distinct areas, while also being close to the knee in 105 the total distance away from the cluster means. We perform the tracking on two of these clusters, indicated in Fig. 11, which are close geographically, but have contrasting snowfall trends.
We perform the moisture tracking by releasing moisture parcels from the target area at random positions within the area and advecting them backwards in time using interpolated wind fields on fixed pressure levels. The amount of evaporation A (mm) 110 that contributes to the precipitation in the target area, at a given location x,y and time step t, depends on the evapotranspiration ET (mm), the amount of water in the parcel W parcel (mm), the fraction of water in the parcel that evaporated from the source S target , and the total precipitable water in the column TPW (mm): (1) 115 The amount of water in the parcel is then updated every time step, including the precipitation P that adds to the parcel when moving back in time. (2) 120 The fraction if precipitation in the target area that originates from a certain source area is then updated as follows: We track the parcel until either less than 1% of the target precipitation is tracked to a source area, or the tracking time is more than 30 days.
Within the WRF domain, the parcels are advected and the moisture budget is calculated using the WRF wind fields and water fluxes. When a parcel gets within one degree of the edge of the WRF domain, there is a gradual linear change to a use of ERA-130 Interim reanalysis data to ensure continuous movement of the air parcels over the domain edge. Within one degree distance from the domain edge, the values used to do the moisture tracking are a combination of the WRF and ERA-Interim values: y int = d*y WRF + (1-d)*y ERA , where d is the distance to the edge. Outside of the WRF domain, the ERA-Interim values are used.
We noted that the surface moisture flux in ERA-Interim is on average 50% higher than in WRF, resulting in a higher mean 135 and standard deviation of the moisture sources outside the WRF domain. Unfortunately, this systematic offset between the two datasets cannot be easily remediated. Although this will not change the spatial patterns of the moisture source trends in a major way, the absolute values of the trends will be lower in the WRF domain than outside. The trends in the Tarim basin will hence be underestimated with respect to regions such as the Caspian Sea and the Junggar basin.

Statistics
Pearson correlation coefficients are calculated using the vectors of annual or seasonal mean values. The trends shown in Figs.
2 and 4 are the slopes from linear fits to the annual mean vectors.

Validation and comparison 145
Any attempt to understand the Karakoram anomaly is greatly hampered by the almost complete lack of in situ meteorological data in WKSK. The rare weather stations in the region are often situated at relatively low elevation, or in urban environments, and poorly represent the high mountain climate. Furthermore, different types of precipitation datasets seem to greatly underestimate the precipitation in mountainous terrain (Immerzeel, Wanders, Lutz, Shea, & Bierkens, 2015;Ménégoz, Gallée, & Jacobi, 2013;Palazzi, Von Hardenberg, & Provenzale, 2013). These complications imply that any meteorological dataset, 150 including reanalysis datasets, are associated with relatively large fundamental errors in WKSK, which prevents reliable validation of any model of WKSK, such as the one presented in here. Nevertheless, we compared our WRF output with data of the region around WKSK, to ensure that the WRF output is a reasonable representation of the regional climate between Another parameter that is important in our model is evapotranspiration. It cannot be directly measured remotely, but there are 175 several datasets that calculate it from other remotely sensed products, either directly or through data assimilation. These datasets vary greatly, as we illustrate in Fig. 4 for July 2010. We show evapotranspiration from GLEAM v 3.3a (Martens et al., 2017;Miralles, Gash, Holmes, De Jeu, & Dolman, 2010), which assimilates various soil moisture, temperature, radiation, and precipitation products. Furthermore, we show SSEBop (Senay, 2018)

Climatic trends
To get an impression how glaciers might have been affected by changes in the climate, we illustrate the trends for two relevant variables: the 2m-temperature in the melt season and the annual snowfall ( For each grid point, the melt season was defined as the months where the mean daily temperature is above -5 °C, since for 195 these months temperatures will likely be above freezing at least part of the time. A threshold value of 0 °C slightly increased the positive temperature trends at lower elevations in WKSK, but meant no trends for the highest elevations could be

220
Our trends show some similarities, but also major differences with respect to a similar WRF study that did not include irrigation and used another re-analysis dataset (Norris et al., 2018). For instance, our temperature trends do not exhibit the strong summer cooling at low altitudes (e.g. the Tarim basin), and are more in line with station data (Waqas & Athar, 2018;Xu, Liu, Fu, & Chen, 2010) in that respect.

Glacier mass balances
The resulting pattern of simulated mass balance (Fig. 8) shows a strong resemblance to the measured pattern of mass balances of recent decades. Most notably, we also obtain growing glaciers in WKSK, whereas the glaciers in other regions show large mass losses. In fact, all points where we model glacier growth in Fig. 3a also show growth or stable conditions in observations (Brun et al., 2017;Kääb et al., 2015), except one point in Kääb et al., (2015). A more detailed quantitative comparison of the 230 above results and the observed mass balances is hampered by the fact that our simulations only go out to 2010, but we compare our results for the period 2000-2008 in Fig. 9. The results generally match reasonably well, although our model seems to show too little growth for the growing glaciers. However, note that the errors on these observations (Brun et al., 2017) are large (~0.3 m w.e.). Furthermore, both the climate model and the glacier model will be associated with errors. By modulating the initial mass balance in the model, we find that on average 41% of the modelled mass balance in 2010 is determined by the 235 initial mass balance in 1980. Although the mass balances in 1980 were observed to be less extreme than in the 21st century (Bolch et al., 2012;Maurer, Schaefer, Rupper, & Corley, 2019), parts of HMA already had negative mass balances then, with the magnitude of initial mass balances generally less than 0.4 m w.e. yr -1 . This would result in an error on the mass balances in Fig. 8

250
We also ran the glacier mass balance model forced by changes in temperature or snowfall only, to disentangle the model sensitivities of the two different variables on the glacier mass balances (Figures 8b and 8c). These results show that the glaciers in the western and southern HMA mainly lose mass due to the increase in temperature, whereas the changes in precipitation have only a minor effect. On the other hand, in the regions where the glaciers are growing, the glaciers are barely affected by 255 the temperature increases in our model. The glacier growth in these regions is mainly caused by an increase in snow (Fig. 8c).
Furthermore, the increase in snow is possibly also responsible for moderating the temperature increases due to the high albedo of fresh snow, which leads to less energy being used for melt. However, the weak temperature response in WKSK is not only caused by the limited temperature trends, but is also due to the limited glacier temperature sensitivity there. We demonstrate this by forcing the glacier model with uniform temperature and precipitation trends (Fig. 10). The reduced temperature 260 sensitivity is in line with previous work (Sakai & Fujita, 2017;Wang, Liu, Shangguan, Radic, & Zhang, 2019). The decrease in snowfall in the western and southern HMA has a far smaller impact on the mass balance than the increase in temperature.
Especially the Himalaya show a low sensitivity to precipitation (Fig. 10). To be able to model thousands of glaciers, our mass balance model is relatively simple and does not solve the full energy balance. A full energy balance model at 1 km resolution has shown that the temperature increases can amplify melt in the monsoon-dominated Himalaya, whereas snowfall increases 265 https://doi.org/10.5194/tc-2019-228 Preprint. Discussion started: 4 November 2019 c Author(s) 2019. CC BY 4.0 License. in the melt season can amplify glacier growth in the Karakoram (Bonekamp, de Kok, Collier, & Immerzeel, 2019). Hence, more detailed models will likely strengthen our conclusion that the observed mass gains are caused by snow increases, whereas the observed mass losses are mainly caused by temperature increases. Unfortunately, modelling the climate and glaciers of the entire HMA at a sub-kilometre resolution for 30 years is currently beyond our capabilities.

Moisture sources 275
The trends in moisture source regions for WKSK (Fig. 11a,b) indicate that the largest increases in moisture from a given source to precipitation in WKSK occur in the mountains themselves. This increase in recycling occurs mainly in May, and is also the main cause of the increase in precipitation in September (see Fig. 12). The increase in recycling is probably a natural consequence of the increased precipitation there. The regions with the second largest increases are the irrigated areas in the Tarim basin, which contributes mainly in May-July (see Figs. 7 and 12). In July, the increase in Tarim irrigation still contributes 280 to increasing precipitation in WKSK, but it falls more in the form of rain, compared to May, where it is mainly snow (Fig. 7).
Another region that contributes to the increase in precipitation in WKSK is the Junggar basin, northeast of the Tarim basin. This is another arid region that has experienced rapid increases in irrigation. The increases per grid point are lower there, but spread out over a larger area. A final source region with an overall large positive trend is the Caspian Sea and the Caucasus.
Note again that, due to a systematic offset in surface moisture flux between WRF and ERA-Interim, the moisture source trends 285 in the Tarim and HMA are underestimated with respect to the other regions.   These results imply that evapotranspiration from irrigated areas in arid Northwest China play a large role in adding water to parts of HMA and hence to the observed positive mass balances. This is in line with recent work that shows that the recent wetting of Central Asia and the Tarim basin is associated with an increase in evapotranspiration in these regions (Dong et al., 300 2018;Peng & Zhou, 2017;Peng, Zhou, Zhang, & Wu, 2018). The increase in the total evapotranspiration is influenced by the increase in potential evapotranspiration (Fang, Chen, & Li, 2018), increase in water availability (Jian et al., 2018), and increase in irrigated land area. On the interannual timescale, precipitation in WKSK strongly correlates with the moisture source amount in the western Tarim basin (Pearson r=0.96 below 3500 m, r=0.68 for the entire WKSK, as indicated in Fig. 11b) for the entire WKSK), showing the importance of the large-scale weather patterns. For the Junggar basin, this correlation is weaker (r=0.65 below 3500 m, r=-0.14 for the entire WKSK), since this region contributes in winter (Fig. 12).
When performing the moisture tracking for the southwestern part of HMA, where snowfall has generally decreased (Fig.   11c,d), also the Caspian Sea and the Junggar basin positively contribute to the snowfall trend, whereas for these ranges the 310 Tarim basin does not contribute to the snowfall trend, with maximal trends in moisture sources of less than 0.1 mm yr -1 . These results show that the irrigated areas in the Tarim basin are especially important in influencing the moisture supply to the Western Kunlun Shan (de Kok et al., 2018).

Conclusion and discussion
Our simulations, based on ERA-Interim and GLDAS reanalysis data, indicate that an increase of snowfall and a low 315 temperature sensitivity are the main reasons why glaciers are growing or stable in western Kunlun Shan and Karakoram. This is the first time that the observed pattern of glacier mass balances in HMA is reproduced in a consistent way. Although the interannual variability of temperature and precipitation is reasonably reproduced, our models are associated with uncertainties, which are partly irreconcilable due to a lack of in situ measurements in WKSK. Furthermore, different parameterisations in the regional climate model, different irrigation schemes, and different glacier models will likely yield slightly different results. 320 Using an ensemble of such approaches could be used to assess the robustness of the results presented here in the future.
Furthermore, detailed studies at smaller scales will give more insight into individual glacier behaviour.
We show that the growing irrigated area in the arid region of Northwest China plays an important role in the increase in snowfall in WKSK. Previous studies have already shown that increases of irrigation in Northwest China can add precipitation 325 to neighbouring mountains (Cai et al., 2019;de Kok et al., 2018), but we now show this process of increasing irrigation is also important compared to other changes in the atmosphere over the last few decades. Already before 1980, irrigation has increased in Northwest China (Fang et al., 2018), possibly contributing to the stable glacier conditions then. Future evolution of snowfall in this part of HMA is partly linked to how the irrigated areas develop in the future. Changes in temperature, irrigated area, or irrigation efficiency are therefore important parameters in understanding future run-off from glaciers and snow in WKSK. The 330 increase in water availability for irrigation in Northwest China might be partly the result of the loss of glacier mass in Tien Shan (Dong et al., 2018). The mass loss will first result in an increase in glacier melt run-off into the Tarim basin, but ultimately the run-off will decrease as the glaciers shrink to a small size (Kraaijenbrink et al., 2017). On the other hand, if the primary source of irrigation water is groundwater, the amount of irrigation for the region will also have a limited sustainable or economic level. Once the groundwater is depleted, the glaciers in WKSK will also receive less snowfall from this region, 335 resulting in their retreat. The relative importance of groundwater extraction, melt from Tien Shan, and recycling from WKSK, for water availability in the Tarim is yet unknown and will require future study. Furthermore, improving the estimates of https://doi.org/10.5194/tc-2019-228 Preprint. Discussion started: 4 November 2019 c Author(s) 2019. CC BY 4.0 License. irrigation gifts, e.g. by remote sensing, could also improve the past climate reconstruction of WKSK. Greening and warming in West-Asia could provide additional snowfall to WKSK, together with an increase in westerly disturbances (Cannon et al., 2014;Kapnick et al., 2014), but if temperatures in HMA keep increasing, the increase in melt will probably counteract glacier 340 growth in most of HMA in the long term. It is clear that the coupling between glacier mass balance, runoff, and irrigation in different regions creates a complex problem of water availability, which will need to be researched further to inform decision makers on irrigation policies.

Code and data availability
The data underlying our results in Figs. 5-12, i.e. monthly mean output from WRF of temperature and precipitation, annual glacier mass balances, and annual moisture sources, will be directly accessible at dataverse.nl (https://hdl.handle.net/10411/ATONZD). Other data is available from the authors upon request. WRF and the glacier mass balance model are freely available. The moisture tracking model is available upon request from Obbe A. Tuinenburg. 355 Author contributions R.J.d.K. and W.W.I. designed the study, with input from all authors. R.J.d.K. performed the WRF modelling, P.D.A.K.
performed the glacier mass balance modelling, and O.A.T. performed the moisture tracking. All authors contributed to the writing and editing of the manuscript.

Competing interests 360
The authors declare that they have no conflict of interest.