Snowmelt and Snow Sublimation in the Indus Basin

: The Indus basin is considered as the one with the highest dependence on snowmelt runoff in High Mountain Asia. The recent High Mountain Asia snow reanalysis enables us to go beyond previous studies by evaluating both snowmelt and snow sublimation at the basin scale. Over 2000– 2016, basin-average snowmelt was 101 ± 11 Gt.a − 1 (121 ± 13 mm.a − 1 ), which represents about 25–30% of basin-average annual precipitation. Snow sublimation accounts for 11% of the mean annual snow ablation, but with a large spatial variability across the basin.


Introduction
Available water resources in the Indus river basin are abstracted almost entirely, mostly for crop irrigation in Pakistan [1]. Whereas the population and crops are concentrated in the low-elevation plain of the Indus, a large part of runoff and groundwater recharge comes from the high elevation regions of the basin in the Hindu Kush, Karakoram and the Himalayas [2][3][4]. This is due to the orographic enhancement of precipitation and the reduction in evaporative demand with elevation ( Figure 1). In the Indus basin, however, the precipitation trend with elevation is not monotonous because a vast high elevation region in the north east lies in the rain shadow area of the Himalayas.
In the upper Indus basin, snowmelt is the dominant contributor to water resources, providing 74% of annual runoff above 2000 m during the period 2001-2014 (23% rainfall, 3% ice melt) [5]. Another study estimated the contribution of snowmelt to annual streamflow to 55.4 ± 2.4% [6], but this estimate excluded seasonal snowmelt from the glacier area, which may explain the lower percentage. Snowmelt contribution to runoff is also significant in the high elevation regions of the Ganges and Brahmaputra basins, but not as much as in the Indus [5][6][7]. Above estimates of snowmelt in the Indus basin were obtained using temperature index models [5,6]. These empirical models rely on air temperature to compute snowmelt while the energy balance of the snow cover is also controlled by additional meteorological variables such as shortwave and longwave radiation, wind speed and air humidity [8,9]. Another limitation of such models is that they do not account for the sublimation of the snow cover. Field studies suggest that sublimation can be a significant component of the snowpack mass balance at the highest elevations of the Himalayas [10,11].
The recent release of the High Mountain Asia UCLA Daily Snow Reanalysis (HMASR) [12,13] enables us to advance knowledge on snow-related hydrological processes in the Indus basin. This reanalysis is based on an energy balance model, which computes sublimation unlike previous studies [14]. Another advantage is that this reanalysis also corrects biases in the precipitation input data by assimilating remote sensing observations of the snow covered area [15]. Hence, this dataset provides the opportunity to reevaluate snowmelt and to evaluate snow sublimation in the Indus basin. Indus basin map and distribution of population density, annual precipitation and potential evaporation by elevation. The basin polygon was extracted from the Global Runoff Data Centre [16] and modified in the northeast following [17]. The population density was sourced from the Gridded Population of the World, Version 4 [18]. The precipitation and evapotranspiration data were sourced from TerraClimate [19]. Base map from OpenStreetMap.

Materials and Methods
The HMASR provides several snow cover variables (albedo, snow depth, snowmelt, snow water equivalent and snow sublimation) at a spatial resolution of 16 arc seconds (~500 m) and a temporal resolution of 1 day. The dataset spans the entire High Mountain Asia region including the headwaters of the Indus basin. The lower regions of the Indus basin (tile-average elevation lower than 1500 m) were not processed due to the unlikely occurrence of snowfall [13]. The dataset covers the period from 1 October 1999 to 30 September 2017.
The HMASR products are distributed by tile of 1° longitude by 1° latitude. Each tile is a netCDF file containing the daily estimates of a variable over a water year (starting on 1 October) [12]. Using netCDF Operators (NCO) [20], Geospatial Data Abstraction Library (GDAL) [21], Orfeo Toolbox (OTB) [22] and GNU parallel [23], I merged the tiles containing the daily snowmelt and sublimation estimates, cropped the resulting multiband mosaics to the Indus river basin outline and computed temporal and spatial averages (source code is provided in the Data Availability section). I used the Indus basin shapefile provided by the Global Runoff Data Centre [16], which I modified to remove a small region in the northeast following [5,17]. The modified basin has an area of 8.31  10 5 km 2 instead of 8.65 × 10 5 km 2 . To avoid latitudinal distortion in spatial average computation, I also resampled the data from their original geographic coordinate system to an equal area projection system. I used a custom equal area coordinate reference system centered over the HMA region. Since the objective of this study is to compute basin-averaged variables and analyze large scale spatial patterns, the resampling method to interpolate the data on the equal area grid is expected to have a marginal impact on the results. Thus, I used the nearest neighbor resampling method to reduce the computation cost. The snowmelt and sublimation were aggregated to daily basin-scale fluxes (units: gigatons per day, equivalent to km 3 of liquid water per day) and time cumulated specific fluxes (units: millimeters water equivalent since 1 October) using the average of all valid pixels in the basin. The Figure 1. Indus basin map and distribution of population density, annual precipitation and potential evaporation by elevation. The basin polygon was extracted from the Global Runoff Data Centre [16] and modified in the northeast following [17]. The population density was sourced from the Gridded Population of the World, Version 4 [18]. The precipitation and evapotranspiration data were sourced from TerraClimate [19]. Base map from OpenStreetMap.

Materials and Methods
The HMASR provides several snow cover variables (albedo, snow depth, snowmelt, snow water equivalent and snow sublimation) at a spatial resolution of 16 arc seconds (~500 m) and a temporal resolution of 1 day. The dataset spans the entire High Mountain Asia region including the headwaters of the Indus basin. The lower regions of the Indus basin (tile-average elevation lower than 1500 m) were not processed due to the unlikely occurrence of snowfall [13]. The dataset covers the period from 1 October 1999 to 30 September 2017.
The HMASR products are distributed by tile of 1 • longitude by 1 • latitude. Each tile is a netCDF file containing the daily estimates of a variable over a water year (starting on 1 October) [12]. Using netCDF Operators (NCO) [20], Geospatial Data Abstraction Library (GDAL) [21], Orfeo Toolbox (OTB) [22] and GNU parallel [23], I merged the tiles containing the daily snowmelt and sublimation estimates, cropped the resulting multiband mosaics to the Indus river basin outline and computed temporal and spatial averages (source code is provided in the Data Availability section). I used the Indus basin shapefile provided by the Global Runoff Data Centre [16], which I modified to remove a small region in the northeast following [5,17]. The modified basin has an area of 8.31 × 10 5 km 2 instead of 8.65 × 10 5 km 2 . To avoid latitudinal distortion in spatial average computation, I also resampled the data from their original geographic coordinate system to an equal area projection system. I used a custom equal area coordinate reference system centered over the HMA region. Since the objective of this study is to compute basin-averaged variables and analyze large scale spatial patterns, the resampling method to interpolate the data on the equal area grid is expected to have a marginal impact on the results. Thus, I used the nearest neighbor resampling method to reduce the computation cost. The snowmelt and sublimation were aggregated to daily basin-scale fluxes (units: gigatons per day, equivalent to km 3 of liquid water per day) and time cumulated specific fluxes (units: millimeters water equivalent since 1 October) using the average of all valid pixels in the basin. The resulting timeseries were plotted by day of water year (DOWY), i.e., day since 1 October. In addition, all the daily products were aggregated along the time axis to characterize the spatial variability of mean annual snowmelt and sublimation in the basin. The sublimation ratio was defined as the fraction of snow ablation due to sublimation, where ablation is the sum of snowmelt and sublimation. is 121 mm with a standard deviation of 13 mm. This corresponds to 101 ± 11 gigatons of water per year. Snowmelt is minimal in December (between DOWY~60 and~90) and starts to rise in early February (near DOWY 150). The basin-scale annual maximum daily snowmelt is reached on DOWY 225, i.e., on 14 May ± 20 days (standard deviation). However, Figure 3 shows that peak snowmelt dates vary across the basin depending on elevation (earlier dates at low elevation) and following an east-west gradient (earlier dates in the west).

Snowmelt
Water 2021, 13,2621 resulting timeseries were plotted by day of water year (DOWY), i.e., day since 1 Oc In addition, all the daily products were aggregated along the time axis to characteri spatial variability of mean annual snowmelt and sublimation in the basin. The su tion ratio was defined as the fraction of snow ablation due to sublimation, where ab is the sum of snowmelt and sublimation.  (standard deviation). ever, Figure 3 shows that peak snowmelt dates vary across the basin depending on tion (earlier dates at low elevation) and following an east-west gradient (earlier da the west).   Figure 4 shows that the total sublimation ranges from 10 to 17 mm.a −1 over the period, with an average value ± standard deviation of 15 ± 2 mm.a −1 . This correspo 12 ± 2 gigatons of water that return to the atmosphere every year in average. The a daily sublimation peak occurs on 15 March ± 46 days. This is about a month earlie the peak snowmelt but with a larger variability, since the standard deviation of the a maximum date is more than twice the one obtained for the snowmelt. In contras snowmelt, daily sublimation rates are minimal during August (~DOWY 300 to 330) resulting timeseries were plotted by day of water year (DOWY), i.e., day since 1 October. In addition, all the daily products were aggregated along the time axis to characterize the spatial variability of mean annual snowmelt and sublimation in the basin. The sublimation ratio was defined as the fraction of snow ablation due to sublimation, where ablation is the sum of snowmelt and sublimation. is 121 mm with a standard deviation of 13 mm. This corresponds to 101 ± 11 gigatons of water per year. Snowmelt is minimal in December (between DOWY ~ 60 and ~ 90) and starts to rise in early February (near DOWY 150). The basin-scale annual maximum daily snowmelt is reached on DOWY 225, i.e., on 14 May ± 20 days (standard deviation). However, Figure 3 shows that peak snowmelt dates vary across the basin depending on elevation (earlier dates at low elevation) and following an east-west gradient (earlier dates in the west).   Figure 4 shows that the total sublimation ranges from 10 to 17 mm.a −1 over the study period, with an average value ± standard deviation of 15 ± 2 mm.a −1 . This corresponds to 12 ± 2 gigatons of water that return to the atmosphere every year in average. The annual daily sublimation peak occurs on 15 March ± 46 days. This is about a month earlier than the peak snowmelt but with a larger variability, since the standard deviation of the annual maximum date is more than twice the one obtained for the snowmelt. In contrast with snowmelt, daily sublimation rates are minimal during August (~DOWY 300 to 330).  Figure 4 shows that the total sublimation ranges from 10 to 17 mm.a −1 over the study period, with an average value ± standard deviation of 15 ± 2 mm.a −1 . This corresponds to 12 ± 2 gigatons of water that return to the atmosphere every year in average. The annual daily sublimation peak occurs on 15 March ± 46 days. This is about a month earlier than the peak snowmelt but with a larger variability, since the standard deviation of the annual maximum date is more than twice the one obtained for the snowmelt. In contrast with snowmelt, daily sublimation rates are minimal during August (DOWY~300 to~330).  Figure 5 shows the relative contribution of snowmelt and sublimation to the seasonal snow cover ablation. The annual sublimation ratio is relatively stable between 8% and 14% over the study period (average 11%). However, the daily sublimation ratio is highly variable, since it can vary from 0% in August to 50% in December. The sublimation generally peaked in April, about a month before snowmelt.    Figure 5 shows the relative contribution of snowmelt and sublimation to the seasonal snow cover ablation. The annual sublimation ratio is relatively stable between 8% and 14% over the study period (average 11%). However, the daily sublimation ratio is highly variable, since it can vary from 0% in August to 50% in December. The sublimation generally peaked in April, about a month before snowmelt.   Figure 5 shows the relative contribution of snowmelt and sublimation to the seasonal snow cover ablation. The annual sublimation ratio is relatively stable between 8% and 14% over the study period (average 11%). However, the daily sublimation ratio is highly variable, since it can vary from 0% in August to 50% in December. The sublimation generally peaked in April, about a month before snowmelt.

Discussion
The mean annual snow ablation (sum of snowmelt and sublimation) is 113 Gt (standard deviation 12 Gt), which is close but lower than the annual snowfall value of 121.8 ± 7.7 Gt estimated by [6] over the period 1979-2019. This discrepancy can be due to (i) different periods of study (2000-2016 versus 1979-2019), (ii) differences in the basin delimitation [17] and (iii) net accumulation of snowfall in the glacier accumulation areas.
The mean annual snowmelt value of 121 ± 13 mm represents a significant but not dominant fraction of the Indus basin water balance (Figure 7). In comparison, average basin precipitation from the WorldClim database is 424 mm [24]. Although there is probably a large uncertainty on this reference precipitation value [3], it falls within the range of previous estimates (392 to 461 mm) reported by [1]. The basin-average rainfall can be roughly estimated to 300 mm by subtracting snowmelt and sublimation from this precipitation value, which means that snowmelt accounts for about 40% of rainfall in the Indus basin.
Although snowmelt is lower than rainfall in the Indus basin, it largely exceeds glacier melt over the same period (2000-2016) [25]. The imbalance ablation (i.e., the release of water due to glacier shrinkage) was estimated to 6.2 Gt.a −1 . Neglecting glacier mass losses due to ice and snow sublimation, the glacier imbalance ablation would represent 4% of the runoff generated from snowmelt at the basin scale. The total glacier ablation (17.7 Gt.a −1 ) from [25] is also reported in Figure 5 but it is less directly comparable to snowmelt estimated here because (i) both estimates include the melting of seasonal snowpack on the glacier surface and (ii) glacier ablation also includes a fraction of sublimation. Additionally, the glacier ablation study [25] did not use the same Indus basin as delimitated by [17].

Discussion
The mean annual snow ablation (sum of snowmelt and sublimation) is 113 Gt (standard deviation 12 Gt), which is close but lower than the annual snowfall value of 121.8 ± 7.7 Gt estimated by [6] over the period 1979-2019. This discrepancy can be due to (i) different periods of study (2000-2016 versus 1979-2019), (ii) differences in the basin delimitation [17] and (iii) net accumulation of snowfall in the glacier accumulation areas.
The mean annual snowmelt value of 121 ± 13 mm represents a significant but not dominant fraction of the Indus basin water balance (Figure 7). In comparison, average basin precipitation from the WorldClim database is 424 mm [24]. Although there is probably a large uncertainty on this reference precipitation value [3], it falls within the range of previous estimates (392 to 461 mm) reported by [1]. The basin-average rainfall can be roughly estimated to 300 mm by subtracting snowmelt and sublimation from this precipitation value, which means that snowmelt accounts for about 40% of rainfall in the Indus basin.
Although snowmelt is lower than rainfall in the Indus basin, it largely exceeds glacier melt over the same period (2000-2016) [25]. The imbalance ablation (i.e., the release of water due to glacier shrinkage) was estimated to 6.2 Gt.a −1 . Neglecting glacier mass losses due to ice and snow sublimation, the glacier imbalance ablation would represent 4% of the runoff generated from snowmelt at the basin scale. The total glacier ablation (17.7 Gt.a −1 ) from [25] is also reported in Figure 5 but it is less directly comparable to snowmelt estimated here because (i) both estimates include the melting of seasonal snowpack on the glacier surface and (ii) glacier ablation also includes a fraction of sublimation. Additionally, the glacier ablation study [25] did not use the same Indus basin as delimitated by [17]. The seasonal cycle of snowmelt is defined by a single peak between May and July, which is similar to previous studies [5,6], but the peak seems to occur earlier than the snowmelt hydrograph simulated by [6], and it seems to be more concentrated in time than the one obtained by [5]. In the Indus, snowmelt rates are high during the summer, which is also the monsoon season. During that period, the combination of heavy rainfalls and meltwater runoff can cause flooding [1]. However, snowmelt rates remain high after the monsoon, hence contributing to sustain river flow after the monsoon flood.
The spatial variability of the peak snowmelt date primarily reflects the effect of elevation but also contrasting climate seasonality between the humid part of the basin in the southwest and the arid region in the northeast. This spatial heterogeneity is also evident in the distribution of the sublimation ratio, which is the highest in the dry and cold regions in the northeast of the basin.
In situ measurements of the sublimation are extremely scarce in the High Mountain Asia. Snow sublimation measured on Yala glacier, Nepal at 5350 m a.s.l. (outside the study area) was 125 mm over a winter season [10]. This is on the upper end of the HMASR sublimation rates at similar elevation ( Figure 8). HMASR sublimation rates can reach 100 mm per year between 6000 and 7000 m a.s.l. in the Indus basin.  The seasonal cycle of snowmelt is defined by a single peak between May and July, which is similar to previous studies [5,6], but the peak seems to occur earlier than the snowmelt hydrograph simulated by [6], and it seems to be more concentrated in time than the one obtained by [5]. In the Indus, snowmelt rates are high during the summer, which is also the monsoon season. During that period, the combination of heavy rainfalls and meltwater runoff can cause flooding [1]. However, snowmelt rates remain high after the monsoon, hence contributing to sustain river flow after the monsoon flood.
The spatial variability of the peak snowmelt date primarily reflects the effect of elevation but also contrasting climate seasonality between the humid part of the basin in the southwest and the arid region in the northeast. This spatial heterogeneity is also evident in the distribution of the sublimation ratio, which is the highest in the dry and cold regions in the northeast of the basin.
In situ measurements of the sublimation are extremely scarce in the High Mountain Asia. Snow sublimation measured on Yala glacier, Nepal at 5350 m a.s.l. (outside the study area) was 125 mm over a winter season [10]. This is on the upper end of the HMASR sublimation rates at similar elevation ( Figure 8). HMASR sublimation rates can reach 100 mm per year between 6000 and 7000 m a.s.l. in the Indus basin. The seasonal cycle of snowmelt is defined by a single peak between May and July, which is similar to previous studies [5,6], but the peak seems to occur earlier than the snowmelt hydrograph simulated by [6], and it seems to be more concentrated in time than the one obtained by [5]. In the Indus, snowmelt rates are high during the summer, which is also the monsoon season. During that period, the combination of heavy rainfalls and meltwater runoff can cause flooding [1]. However, snowmelt rates remain high after the monsoon, hence contributing to sustain river flow after the monsoon flood.
The spatial variability of the peak snowmelt date primarily reflects the effect of elevation but also contrasting climate seasonality between the humid part of the basin in the southwest and the arid region in the northeast. This spatial heterogeneity is also evident in the distribution of the sublimation ratio, which is the highest in the dry and cold regions in the northeast of the basin.
In situ measurements of the sublimation are extremely scarce in the High Mountain Asia. Snow sublimation measured on Yala glacier, Nepal at 5350 m a.s.l. (outside the study area) was 125 mm over a winter season [10]. This is on the upper end of the HMASR sublimation rates at similar elevation ( Figure 8). HMASR sublimation rates can reach 100 mm per year between 6000 and 7000 m a.s.l. in the Indus basin.  The HMASR is based on the SSiB 3-layer snow scheme, which computes latent heat fluxes from both the ground and canopy but neglects blowing snow processes [26,27]. Hence, the sublimation rates reported above correspond to static surface sublimation. Blowing snow sublimation can be a significant component of the snowpack mass balance in high mountain regions [28][29][30][31][32]. However, the few available studies on this topic suggest that the relative contribution of static surface and blowing snow sublimation is highly variable depending on the climate, topography and forest cover. For example, activating blowing snow sublimation resulted in a threefold increase in the total sublimation over a 22 h blowing snow event in [30], but had a negligible effect on the total sublimation over a hydrological year in [29]. The latter result may be understood by the fact that sublimation is limited by the available energy and vapor-pressure deficit in the atmosphere so that the activation of blowing snow sublimation only modified the partitioning of the latent heat flux between static surface sublimation and blowing snow sublimation, but not the total mass loss due to sublimation. Therefore, further studies are needed to characterize the hydrological significance of blowing snow physics in the Indus basin.
The effect of the vegetation on the snow melt and sublimation patterns was not specifically addressed in this study due to the small fraction of the forested areas in the basin. Considering 2000 m as the minimum elevation of the snow dominated regions of the basin [5], the forest cover accounts for 9% only according to the Copernicus Global Land Cover Layers [33]. The dominant land cover are bare/sparse vegetation (40%) and herbaceous vegetation (24%), which are expected to have a low impact on snow cover processes.

Conclusions
The HMASR provides the opportunity to better characterize snow processes and the hydrological balance in the Indus river basin. Snowmelt is a significant but not dominant input to the water balance of the Indus river basin over the period 2000-2016. In addition, the melt season overlaps with the monsoon season, when the issue is more the excess of water than a lack of water. The data indicate that 11% of the seasonal snow is "lost" by sublimation over this period, and this value might be underestimated by the lack of blowing snow physics in the model. Although the effects of climate change on snow accumulation and melt in this region are now well understood [6], it remains to evaluate how the sublimation fraction could evolve under future climate scenarios.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The HMASR is freely available from the National Snow and Ice Data Center [12]. This study can be replicated and applied to other river basins of High Mountain Asia using the scripts provided in this repository: https://github.com/sgascoin/HMA-Snow-Reanalysisscripts (accessed on 21 September 2021).