Indus River Basin Glacier Melt at the Subbasin Scale

Pakistan is the most glaciated country on the planet but faces increasing water scarcity due to the vulnerability of its primary water source, the Indus River, to changes in climate and demand. Glacier melt constitutes over one-third of the Indus River’s discharge, but the impacts of glacier shrinkage from anthropogenic climate change are not equal across all eleven subbasins of the Upper Indus. We present an exploration of glacier melt contribution to Indus River flow at the subbasin scale using a distributed surface energy and mass balance model run 2001–2013 and calibrated with geodetic mass balance data. We find that the northern subbasins, the three in the Karakoram Range, contribute more glacier meltwater than the other basins combined. While glacier melt discharge tends to be large where there are more glaciers, our modeling study reveals that glacier melt does not scale directly with glaciated area. The largest volume of glacier melt comes from the Gilgit/Hunza subbasin, whose glaciers are at lower elevations than the other Karakoram subbasins. Regional application of the model allows an assessment of the dominant drivers of melt and their spatial distributions. Melt energy in the Nubra/Shyok and neighboring Zaskar subbasins is dominated by radiative fluxes, while turbulent fluxes dominate the melt signal in the west and south. This study provides a theoretical exploration of the spatial patterns to glacier melt in the Upper Indus Basin, a critical foundation for understanding when glaciers melt, information that can inform projections of water supply and scarcity in Pakistan.


INTRODUCTION
Over 300 million people receive water from the Indus River (Khan and Adams, 2019), with tributaries in China, Afghanistan, India, Pakistan, and Kashmir (Sattar et al., 2017). Because many of the international boundaries in the Indus Basin are disputed ones, it is difficult to partition the Indus Basin between different countries exactly. However, it is clear that while the Indus River flows through significant stretches of China and India before reaching Pakistan, it also traverses the entire length of Pakistan, from the Karakoram Mountains in the north to the Arabian Sea in the south.
The Indus is often characterized as the "lifeblood" of Pakistan (e.g., Craig et al., 2019). The vast majority of the country (~92%) is semi-arid to arid (Food and Agriculture Organization of the United Nations (FAO), 2011a) yet boasts the world's largest contiguous irrigation network (FAO, 2011b). Irrigated agriculture contributes a quarter of Pakistan's GDP and employs roughly half of the country's labor force (Briscoe et al., 2005). Irrigation in Pakistan is mostly (94%) canal networks, through which water flows by gravity, while groundwater supplies a minute portion of the country's irrigation (Siyal et al., 2021).
With withdrawals increasingly exceeding sustainable supply, overallocation presents significant challenges to Indus supply and reliability (Sattar et al., 2017). So, too, does climate change, which affects both the demand for and supply of water. Importantly, the river is sourced primarily from snow and glacier melt in the northern, mountainous Upper Indus Basin (UIB) (e.g., Immerzeel et al., 2010;Tahir et al., 2017).
The share of Indus River flow from glacier melt is greater for the Indus than for any other single glaciated river basin worldwide (Immerzeel et al., 2020); one estimate for the UIB places the glacier melt contribution at 40.6% of runoff (Lutz et al., 2014). For a review of such estimates, see Azam et al. (2021)'s Table 2. Estimates on the combined snow and glacier melt contribution range upwards of 72% (Immerzeel et al., 2009). Cryospheric meltwater contributes substantially to Indus River flow because most precipitation falls as snow. That precipitation comprises seasonal snowpacks or becomes incorporated into glaciers and melts in the warm season. The cryosphere, thus, modulates the timing of river flows by generating melt between April and October, with the majority of that water coming from glaciers between June and October (Azam et al., 2021).
Due to the extensive and increasing (Azam et al., 2021) demands for its water-from a dense population and 144,900 km 2 of irrigated area (Immerzeel et al., 2010)-the Indus Basin has long been highlighted as the most vulnerable to water availability changes of all ten river basins in High Mountain Asia (Immerzeel and Bierkens, 2012). A growing economy and population depend on non-renewable (e.g., glacier melt) or slowly renewable (e.g., groundwater) water resources in the Indus Basin, giving it the distinction of being "Asia's hotspot" (Immerzeel and Bierkens, 2012) and "globally the most important water tower" (Immerzeel et al., 2020).
The region's climatic and orographic complexity provides a compelling reason to look at variations in the UIB's glaciers on a finer scale (Miller et al., 2012). Precipitation sourced by monsoons and westerlies meet over the UIB (Farhan et al., 2014), and the interaction between precipitation and topography leads to effects at a range of spatial scales (Lutz et al., 2016). The glacier response to climate change in the Himalayan and Karakoram regions has been spatially heterogeneous and temporally variable (Azam et al., 2021).
Considerable work has been done to investigate the large-scale response of the Indus Basin's water supply to climate change. However, the Indus's dependence on glaciers is significant, and the climate response of glaciers across a region is not uniform (e.g., Rounce et al., 2020a;Shean et al., 2020). The Upper Indus subbasins contain different glacierized areas and generate different amounts of glacier melt; therefore, subbasins contribute differentially to the overall basin's vulnerability. There is a need to distinguish between the second-level basins (i.e., subbasins) in assessing the current and future state of the Indus Basin as a whole. The few existing detailed studies leave open the need for a robust quantification of glacier contributions to Indus flow by subbasin. Mukhopadhyay and Khan (2014) explore the components of flow in the UIB via hydrograph separation at different elevation bands using data from 11 gauging stations over several decades. This and their follow-up study (Mukhopadhyay and Khan, 2015) use correlation-based, semi-quantitative approaches and elucidate the altitude dependence of different runoff sources and their timing. They assume that percent of flow from glacier melt is proportional to glacierized area. Dahri et al. (2016) offer a look at the smaller-scale orographic effects on precipitation in the UIB and demonstrate errors associated with gridded datasets. They attribute hydrograph characteristics to glaciers but do not analyze glacier contributions quantitatively.
Modeling offers a powerful tool to investigate the dynamics of a river system and parse out the various components of its supply. The complex climatology, extreme topography, and difficult field access present challenges to studies of the cryosphere in the UIB region (Lund et al., 2020) and support the benefits of modeling. Of the two types of glacier melt models, energy balance models and temperature index models (Hock, 2003), temperature index models are more commonly applied at large spatial scales (e.g., Radić et al., 2014;Maussion et al., 2019; majority of models compared in Marzeion et al., 2020), but there is often a trade-off between simplicity and reliability (Azam et al., 2021).
We lend a mesoscale perspective to the question of how important glacier melt is to Indus River flow by quantifying relative contributions from the 11 Indus subbasins. We use a physically based, distributed surface energy and mass balance model to distinguish mass balance components and calculate the melt of all clean glaciers over the entire UIB in the detail often reserved for studies of single glaciers or small regions. Our approach adds a modeling-based, spatially detailed perspective to the discussion of how Pakistan's lifeline river is dependent upon glacier melt. Given the uncertainties associated with the model assumptions and the dynamically downscaled meteorological data used to drive it, our results provide valuable insight into glacier melt in the UIB but do not comprise a definitive quantification of such melt.

MATERIALS AND METHODS
To assess the contribution of glacier melt to Indus River flow by subbasin, we first defined the extent of the UIB and its subbasins.
Then, we used a gridded atmospheric data product to drive a glacier surface energy and mass balance model calibrated by geodetic mass balance data.

Study Area
We used the Randolph Glacier Inventory (RGI) version 6.0 (RGI Consortium, 2017) to identify glaciers within the extent of the UIB as it was defined by Wijngaard et al. (2017), which we emulated via the ArcGIS Watershed tool. Collectively, there are almost twenty-five thousand km 2 of glacierized area and over seventeen thousand individual glaciers in the UIB. The UIB extent varies in the literature because many define it as limited to the area whose water drains into the lake behind Tarbela Dam, which stores water before the river enters the lowlands in the Khyber Pakhtunkhwa province (e.g., Immerzeel et al., 2009;Mukhopadhyay and Khan, 2014;Tahir et al., 2017). We take a pure definition of the UIB: the mountainous area whose runoff ultimately drains into the Indus River (the Lower Indus Basin is the part of the Indus Basin in the lower-lying regions; e.g., Siyal et al., 2021). Our definition includes the part of the UIB in the Hindu Kush (the Kabul/Swat/Alingar subbasin in Figure 1), while UIB extents limited to the Tarbela Reservoir's catchment areas exclude it. Wijngaard et al. (2017) determined UIB extent by hydrologically correcting a 5 km digital elevation model (constructed by resampling a NASA SRTM DEM) using HydroSHEDS (Hydrological data and maps based on SHuttle Elevation Derivatives at multiple Scales; Lehner et al., 2008).
Within the UIB of Wijngaard et al. (2017), we followed the subbasin delineations of the HydroSHEDS-derived product used in Immerzeel et al. (2020). We followed the FAO's basin aggregations, such as the Gilgit and Hunza into a single "Gilgit/Hunza" subbasin, and we use the naming conventions of the FAO product in this manuscript. However, we recognize that many basins have alternate names and/or spellings. In the scientific literature, counts for subbasins within the UIB range from 5 (Ahmed et al., 2014) to 13 (Dahri et al., 2016). Following the FAO product as we do, there are 11 subbasins, which vary in size, glaciated area, and number of debris-covered glaciers, as  (Kraaijenbrink et al., 2017;RGI Consortium, 2017;Scherler et al., 2018 Figure 1, and glaciers' subbasins were defined as the basin in which they terminate.

Meteorological Inputs
The High Asia Refined Analysis (HAR; Maussion et al., 2014) is a gridded atmospheric data product generated by dynamical downscaling of the Global Forecasting System data using the Weather Research and Forecasting model. It provides hourly atmospheric data at 30 km resolution for Central Asia and at 10 km for the area surrounding the data-sparse Tibetan Plateau, including the UIB. We use the 10 km resolution HAR version 1, which includes full calendar years 2001-2013. The model described in the subsequent subsection requires the following inputs from HAR: incoming radiative fluxes at the surface, 2 m air temperature, total precipitation, surface pressure, 2 m relative humidity, emissivity, and 10 m windspeed. We imposed a precipitation cutoff for drizzle at 0.25 mm per day, consistent with instrument limits on measurement (following, e.g., Riley et al., 2021). The model is driven by meteorological input from HAR and additional variables derived from them that are necessary to calculate the turbulent and outgoing longwave fluxes (e.g., surface saturation vapor pressure). We elected to use HAR for three reasons: first, it has high temporal and spatial resolution and, thus, captures topography better than most other publicly available products. Second, HAR includes all the variables required to force our glacier surface energy and mass balance model. Finally, its validation across High Mountain Asia has been the subject of an extensive set of prior studies Mölg et al., 2014;Pritchard et al., 2019;Yoon et al., 2019). Although using HAR-or, in fact, any gridded product or even automated weather station network in complex terrain-has significant limitations, these validation studies suggest that use of HAR is justified. They show that HAR precipitation seasonality is well-represented, and precipitation is generally consistent with observed values of runoff.
HAR temperatures match observations well when considered daily and hourly but do show a cold bias in winter over the Tibetan Plateau. Biases vary seasonally, with an underestimation of diurnal temperature ranges in spring, early summer, and autumn, times when humidity biases are observed: overestimation of specific humidity in winter and relative humidity in winter, spring, and autumn. The cold temperature bias reduces in summer, along with the high bias in relative humidity. The HAR summer is characterized by underestimation of longwave and overestimation of shortwave radiation, which are results of too-low cloud cover in HAR. While temperature shows seasonal biases (Pritchard et al., 2019), HAR precipitation shows seasonality and spatial patterns consistent with an observation data set of 31 rain gauges on the Tibetan Plateau complemented by the gauge-calibrated estimates of the Tropical Rainfall Measurement Mission . Finally, there have been energy balance studies specifically using HAR on glaciers. The most directly relevant is that of , in which HAR was input to a full energy and surface mass balance model for Zhadang Glacier in the central Tibetan Plateau  used HAR only-without bias correcting or validating individual energy balance inputs-and found good agreement when compared to forcing the same glacier with nearby AWS data. Although a precipitation scaling factor was used to achieve observed mass balances in the Tibetan Plateau, the work of  does illustrate the utility of our approach.
HAR's strength lies in spatial patterns and tendencies so it is well suited to our regional research question. We do not assert that the HAR product represents reality or is free from biases (summarized above) but use it in this study as the assumed climate for the purpose of our investigation.
Input variables from HAR are either assumed constant across the glacier or adjusted for glacier-wide application, following Johnson and Rupper (2020). For example, the HAR wind speed is scaled from 10 to 2 m height (following, e.g., Oke, 1987;Stull, 1988) and then assumed uniform over the entire glacier surface. Both incoming shortwave and incoming longwave radiation are distributed over the glacier based on physical variables. Shortwave is adjusted via glacier slope, glacier aspect, and sun position following Reda and Andreas (2004), and incoming longwave is distributed using an effective emissivity (Braithwaite and Olesen, 1990). These methods of distributing incoming radiation have been used and validated in other regions of the world (e.g., Arnold et al., 1996;Hock and Holmgren, 2005;Wallace and Hobbs, 2006;Ebrahimi and Marshall, 2015;Johnson and Rupper, 2020). We note that the standard assumption of logarithmic wind speed profiles cannot account for katabatic winds. However, not simulating katabatic winds has an impact on the latent heat flux that is at least compensated for by an overestimation of temperature.
Three inputs' values across a glacier surface are determined by adjusting for the glacier-HAR elevation difference. First, the air pressure, used to calculate turbulent fluxes, is recalculated every timestep from the HAR gridcell elevation to give air pressure at the glacier surface elevation. Second, air temperatures are adjusted using day-of-the-year lapse rates calculated from the 3 × 3 matrix of mean daily temperatures of HAR gridcells centered over the modeled glacier. We calculate the daily climatological lapse rates; thus, the lapse rates are the same for a given date in each year of the model run and are used to scale temperature from the elevation of a glacier's HAR gridcell to its surface. Third, we updated the Johnson and Rupper (2020) model to include a precipitation gradient, reflecting the well-known general trend of precipitation increase with elevation (e.g., Winiger et al., 2005;Hewitt, 2014;Pritchard et al., 2019). The change in precipitation with elevation is computed by subbasin and is detailed in the subsequent section.

Surface Energy and Mass Balance Model
We use the Surface Energy and Mass Balance Model (SEBM) of Johnson and Rupper (2020), here run hourly at a resolution of 90 m over glacier surfaces, to simulate and quantify melt from Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 767411 snow and ice in the glacier system. Here, we give only a brief overview and highlight adaptations made for this study. We refer the reader to Johnson and Rupper (2020) for model details. The net energy, Q net , is the sum of the surface energy fluxes (Eq. 1): SW is shortwave (solar) radiative flux, LW is longwave radiative flux, H and LE are the sensible and latent heat terms of the turbulent energy flux, Q P is energy carried by liquid precipitation, and Q C is the conductive heat flux between the glacier surface and glacier body. While Q P is commonly neglected in SEBMs (e.g., Wagnon et al., 2003;Mölg and Hardy, 2004), it is included here for completeness. We neglect refreeze, and the only ground flux accounted for is the conductive one: Qc = K s (T 2 -T s )/d s , where K s is the thermal conductivity of the surface layer computed from the layer's density using the formula of Van Dusen (1929).
An important determiner of the shortwave radiation balance is the surface albedo which, following the methods of Johnson and Rupper (2020), is calculated as a function of the composition of the surface (snow, firn, or ice; values in Appendix Table A1), snow depth, and time since the last snowfall. When the air temperature is equal to or less than the precipitation phase transition threshold (Appendix Table A1), precipitation falls as snow and resets the albedo of the surface to that of fresh snow.
The net surface energy flux controls the evolution of glacier surface temperature such that a change in surface temperature is: ΔT s Qnet cρ s ds dt, where c is the specific heat capacity of ice (= 2097 J kg −1 K −1 ), ρ s is the surface density computed as a weighted average of the ice and snow contents in the surface layer (Appendix Table A1), and d s is the thickness of the surface layer (22 cm, after Arnold et al., 2006).
The model is a three-layer model, requiring temperature at the surface, the near-surface (22 cm depth in the glacier), and the glacier body (10 m). The model domain ranges from the surface to a depth of 10 m, and the composition of that 10 m evolves over time as snow and ice accumulate, sublimate, and melt. The model prognoses the evolving temperatures at the surface (T s ) and 22 cm depth (T 2 ), but the temperature at 10 m depth (T b ) is assumed constant and equivalent to the long-term mean air temperature (Paterson, 1994). The temperature in the second layer changes only by conduction (i.e., only when there is a temperature gradient).
The model calculates conduction through ice, snow, or a combination of the two-but not debris. Thus, we limit our study to the 81% of UIB glaciers (comprising 76% of the glacierized area) that are "clean" glaciers, defined here as glaciers whose areas are less than 15% covered by debris. Debris cover information for 17,252 of the 17,270 glaciers in the UIB (Table 1) is from Scherler et al. (2018), given by the Landsat eight bands 4:6 ratio algorithm of Hall et al. (1987). For 18 glaciers with incomplete attributes in the RGI, we use glacier area from Shean et al. (2020) and debris cover from Kraaijenbrink et al. (2017).
The surface temperature, when calculated as positive, determines the energy available for melt; physically, when energy is present to raise the surface temperature, that energy is absorbed in the endothermic melting of ice such that the actual surface temperature never exceeds 0°C. Note that warming the surface temperature to the melting point requires energy; however, it is only when the calculated surface temperature exceeds 0°C that energy is available for melt (Eq. 3).
Other significant changes from the Johnson and Rupper (2020) model include the introduction of a precipitation gradient that is consistent with HAR and calculated separately for each subbasin. Precipitation increases with elevation as a function of the amount of precipitation (Eq. 4): where P HAR p 1 h gl + p 2 is the hourly average precipitation (mm h −1 ) during hours when it precipitates. p 1 and p 2 are, respectively, the slope and y-axis intercept of the robust linear bisquare regression of elevation, limited to the range occupied by glaciers, vs. average hourly precipitation when precipitation is non-zero. Eq. 4 shows that the precipitation on the glacier (P in m) is given by the closest HAR gridcell's precipitation (in m) adjusted by the difference in elevation (h in m) multiplied by the share of the P HAR that falls during the timestep and the lapse rate for the year (p 1 , in mm h −1 ). While precipitation generally increases with elevation, precipitation does decrease in some regions above a certain elevation (e.g., Hewitt, 2014). Limited observations and 10 km HAR resolution preclude identifying those transitions in this study; there are multiple approaches to and significant uncertainties associated with downscaling precipitation in complex terrain. The effects of precipitation uncertainty are discussed in Uncertainty. The revised model also includes a stipulation that the stability correction be applied only with the Richardson number (R b ) > 0 and wind speed (U) > 1 m s −1 (Anderson et al., 2010) to accommodate a range of glacier surface temperatures and to prevent unreasonable values at very low wind speeds. The model undergoes a spin-up of 1 year, which is sufficient to establish a depth for the snowpack. Finally, we introduced a correction for large, persistent surface temperature outliers that deviate unphysically due to a combination of forcing variables. These outliers are infrequent and isolated rather than pervasive: though~95% of glaciers have extreme outliers, the outliers occur on average in 0.17% of the timesteps. Accordingly, we reset the nonphysical, outlying temperature values that vacillate more than 15°C from the previous (hourly) time step's temperature. Outlying values are set equal to the average surface temperature over the previous 24 h, a more reasonable value; this ensures results within physical reality for the whole time series. We ran the model at 90 m spatial resolution, having degraded spatial resolution from the 30 m ASTER GDEM V3 to improve run time, and at 1-h time steps over the HAR calendar years 2001-2013. Glacier shapes were given by RGI 6, and glacier mass balance is accumulation less melt and sublimation. Figure 2 shows all three mass balance components in a stacked bar graph, by subbasin. Because evaporation comes from meltwater, we do not double count it in the mass balance calculation. We assume all melt goes to runoff in the model and that, therefore, modeled melt is equivalent to modeled runoff. Additionally, note that, because the Jhelum subbasin includes so few glaciers, we include Jhelum's glaciers with the Krishen Ganga subbasin, and we report model runs and results over 10 subbasins, rather than 11 ( Table 2).

Calibration
We calibrated the model with 2000-2018 geodetic mass balance data from Shean et al. (2020), which was mostly calculated from ASTER imagery. Forcing the model with 2001-2013 HAR output but calibrating to a longer time range is not ideal, as it can introduce additional uncertainty. However, this was the most complete observational dataset for the region at the time the calibration and validation were completed, and, importantly, it includes smaller glaciers. The timeframe differences between modeled and observed mass balance are unlikely to have significant impacts on the results, as discussed further in Uncertainty. We formed a calibration sample of n = 323 glaciers and extrapolated the calibration parameter to all UIB glaciers; calibrating all would have introduced a significant computational burden and, more importantly, would not have been calibrating the variable of interest, melt. While close agreement between modeled and geodetic mass balance values across the entire study area may have improved the certainty of our melt results, calibrating the mass balance as a whole would not necessarily improve certainty in its components (melt, accumulation, sublimation). Unfortunately, measurements of melt are not available, and myriad combinations of melt, accumulation, and sublimation values can sum to the same value of mass balance. Given that it is not possible to calibrate the melt component of mass balance separately, that calibrating all UIB glaciers would not necessarily improve the certainty in our melt result, and that calibrating all glaciers would have presented a significant computational burden, we opted for minimal calibration to avoid over-calibration of our results and artificial inflation of their reliability. To select the calibration sample, we started with 30 glaciers per basin that were spatially distributed and that had distributions of elevation, slope, aspect, and area closely matching the distributions of those characteristics for all glaciers in that basin. Twenty-three additional glaciers filled spatial gaps and ensured representation of all climate settings. The elevations of the 323 sample glaciers are representative of the glacier elevation range for the broader UIB population (n = 13,912) minus its outliers.
Precipitation is a convenient and reasonable parameter to adjust given its coarse resolution and low accuracy in climate data (Rounce et al., 2020a), particularly in the UIB (Dahri et al., 2016) and other parts of High Mountain Asia (e.g., Wortmann et al., 2018). Thus, while multiple calibration methods exist (Ismail et al., 2020), we chose to reduce the model-measurement bias through a precipitation correction factor (PCF), which was the constant, annual amount by which the total precipitation magnitude (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) was adjusted from the downscaled   (Shean et al., 2020). Subbasins are ordered by increasing area to match the order in Figure 5. The average MB and melt are given, as well as corresponding melt volumes and area-weighted means. The RMSE is calculated by comparing a subbasin's individual glaciers' modeled MBs with their geodetic MBs. The percent contribution to total melt volume is provided for the northern subbasins and southern subbasins to illustrate the dominance of the northern regions on melt volume. Note that Krishen Ganga and Jhelum are combined into a single basin (KG&J) for the modeling. K/S/A is Kabul/Swat/Alingar. HAR precipitation. The PCF value was determined by iteration until modelled mass balance was within 10% of the geodetic mass balance, resulting in one PCF for each of the 323 calibration sample glaciers. The PCF is not a precipitation-bias correction as it effectively incorporated all uncertainty-including SEMB model biases, HAR biases, and downscaling biases-into a single parameter. The 323 PCFs, all of whose values are within physically reasonable bounds, were applied to the 13,912-glacier population by taking the average of PCFs for all sample glaciers within a 20 km radius or, if there were fewer than 10 in that radius, the 10 closest sample glaciers (gives UIB-wide PCF μ −0.68 and σ 0.20). This combination, out of nine algorithms identified a priori and shown in Supplementary Table S1, minimized the amount that the PCFs calculated by local averaging varied from the calibrated PCFs for the 323 subset glaciers. The way that the calibration was applied minimizes the impact of model over-fit: the 323 calibrated PCFs were not input directly to the model. Rather, they were used to calculate every glacier's individual PCF such that even the 323 glaciers in the calibration subset were modeled with PCF values calculated using the algorithm above. We calibrated an independent set of glaciers (n = 95) spatially distributed across the UIB to confirm that the choice of glaciers for calibration did not affect the results (Supplementary Figure S2B).
We performed 10-fold cross-validation to determine any bias introduced by our local averaging method to estimate PCFs for all glaciers in the UIB. Results of this cross-validation showed minimal systemic bias (μ = 0.02 m w.e. yr −1 , Supplementary Figure S2A), supporting our use of this method. The variability in calculated biases for the 10-fold samples (σ = 0.1 m w.e. yr −1 ) shows a minor sensitivity to the exact choice of calibration glacier subset. This sensitivity, however, is likely a reflection of the careful selection of the full training data set (n = 323) to be regionally and geomorphically representative, a characteristic the individual cross-validation subsets lack. The k-fold cross-validation, then, confirms the lack of any significant bias introduced either by the local averaging method for PCF estimation or by the specific glaciers used for calibration.

Uncertainty Assessment
Once calibrated, the model calculates energy and mass fluxes hourly over the surface of each glacier for the duration of the HAR forcing. Totaling the mass balance components (with melt being our focus) and aggregating by subbasin demonstrates differences in each basin's glacial response to meteorological forcing given by HAR. Calibration is one of two sources of uncertainty that we quantified and applied to the model's calculation of glacier melt. While the model-measured mass balance bias provided a mass balance uncertainty, it was important not to apply this mass balance uncertainty directly to the melt results because melt is not necessarily off if the mass balance is off (and vice versa). While melt and mass balance uncertainties likely have some relation, they are rarely equivalent.
For glaciers whose mass balance bias is small, but still more than 10%, the uncertainty is assumed to result solely from model parameterization. When mass balance varies greatly from the geodetic mass balance, error stems from calibration-related uncertainty superimposed upon that from parameter choice. We assume that when the bias is due to both, melt uncertainty is a linear function of mass balance. We assume a linear relationship given the absence of a known relation combined with the likelihood that there is generally some relationship between melt uncertainty and mass balance bias. We acknowledge this approach is subjective but use it as a semiquantitative approach to approaching uncertainties.  Appendix Table A1. They are supported by Naegeli and Huss (2017) and Bintanja and Reijmer (2001) Albedo of fresh snow (α fs ) 0.85 0.69 Kayastha et al. (1999), Paterson (1994) Roughness length for wind (z 0m ) snow (0.0005 m) and ice (0.00689 m) snow (0.0015 m) and ice (0.016 m) snow (Kayastha et al., 1999;Wagnon et al., 2003) and ice (Brocket al., 2006;Azam et al., 2014) Roughness length for temperature (z 0T ) snow (0.0008 m) and ice (0.0038 m) snow (0.0012 m) and ice (0.0042 m) The values for the temperature roughness lengths are within 1 standard deviation of assumed normal distributions about the values in Appendix Table A1.
Roughness length for vapor pressure (z 0q ) snow (0.0008 m) and ice (0.0038 m) snow (0.0012 m) and ice (0.0042 m) The values for the water vapor roughness lengths are within 1 standard deviation of assumed normal distributions about the values in Appendix Table A1. To determine the thresholds for the different categories of uncertainty and to quantify melt uncertainty, we ran model tests of 12 glaciers with different parameter and PCF inputs, as described in the subsequent subsections and Supplementary Material. The 12 glaciers were selected because they spanned all subbasins and represented the large range of glacier sizes in the UIB: half were large outliers in a distribution of glacier areas (with outlier defined using the medcouple-based definition of Hubert and Vandervieren, 2008). Modeling the 12 glaciers in the test set yielded quantification of model sensitivity to both the input parameters and the calibration through the following approaches. We apply the melt uncertainty relationships derived from the 12 glaciers to all glaciers, by subbasin, to compute uncertainty in our results.

Melt Uncertainty From Model Parameter Values
Of the parameters used in the SEBM (see Appendix Table A1), any could potentially factor into model uncertainty. However, the ten in Table 3 stand out as relatively less certain and/or more likely to elicit model sensitivity (Braithwaite, 1995;Moore, 1983;Klok and Oerlemans, 2004;Immerzeel et al., 2014;Shea et al., 2015;Giese et al., 2020;Johnson and Rupper, 2020). The albedos of bare ice and fresh snow directly affect the absorption of energy by the surface, while the temperature lapse rate and temperature at which precipitation phase transitions affect whether precipitation falls as rain or snow. This, in turn, controls the surface albedo and the accumulation component of the modelled mass balance. We also tested sensitivity to the six roughness lengths-three each for snow and ice-which are highly variable and difficult to measure directly (Smith, 2014). Roughness lengths are the height above the surface at which the surface itself no longer affects the variable in question: momentum (i.e., wind speed), temperature, or water vapor.
There were three sets of parameters for the uncertainty analysis: the original ones (Appendix Table A1 or, for lapse rate, the model's calculated value) and ones specified to minimize and maximize melt ( Table 3). For the 12 glaciers, calibration to within 10% of the geodetic mass balance under three parameter sets gave a range of melt variation (i.e., uncertainty) due solely to the model parameters.

Melt Uncertainty From Calibration
After obtaining the mass balances (MBs) in the previous section that were within 10% of the geodetic values, we altered each glacier's PCF by one standard deviation of the calibration subset's PCFs (n = 323 glaciers) under each of the three parameter sets. This yielded probable ranges of melt and mass balance due to applying the 323 PCFs to all clean glaciers in the UIB. Because applying the PCF did not include any glacier-specific weighting and was instead the mean of PCFs for the nearby glaciers in the calibration sample, we could perturb all glaciers in the uncertainty analysis with the same change in PCF (here, ± σ). This perturbation revealed that the precipitation adjustment (i.e., PCF) affects the components of mass balance equally; for the 12 test glaciers, the melt response range was 49.85% of the full MB response range (see Supplementary Figure S4).

Error/Uncertainty
We divide our model results into three categories based on how well the model performs, and we use the 12-glacier set to quantify melt uncertainties in each (a more detailed version of what follows, including a more detailed version of the figure, can be found in the Supplementary Material). These were 1) model (parameter) uncertainty only; 2) model + calibration uncertainty; and 3) model + greater calibration uncertainty for glaciers whose mass balance deviation from measured values was classified as an outlier (details are provided in Figure 3 and the Supplementary Material). Calibration was the biggest source of melt uncertainty in our results. Calibration adds uncertainty on top of the relatively small amount from model parameters, but the two are not merely additive because parameters may affect model sensitivity to calibration.
The melt uncertainty analysis performed here is likely a reasonable first-order estimate and is, therefore, useful for assessing the relative difference in and sources of melt uncertainties across the UIB. In the subsequent sections, we present and focus on our results, which are robust given the uncertainties in modeled melt.

Melt
The SEBM applied to the UIB's clean glaciers elucidates regional patterns in melt energy and its components. We emphasize again here that the significant paucity in observations precludes a complete calibration and validation of melt in this region. Thus, these results should be viewed as a theoretical assessment of plausible variations in melt, rather than a definitive quantification of melt. Within this theoretical framework, the amount of melt on any glacier is a complex function of the glacier's geomorphic characteristics, surface properties, and climate. Across the UIB, glaciers melt for different proportions of the 13-years model run ( Figure 4A). Days with melt correlate with mean daily energy available for melt ( Figure 4B), which tends to increase towards the southwest. The average daily melt energy is reflected in the area-weighted mean melt of each subbasin ( Table 2).
Melt has a stronger correlation (R 2 = 0.35) with elevation than any other component of the mass balance or even the mass balance itself. The glaciers with high melt energy tend to be those located at lower elevations, but, with an R 2 = 0.35, elevation is only a weak predictor of melt. Figures 4C,D show the relative contributions of daily mean radiative and turbulent fluxes to the daily mean melt energy on days with melt. We do not show energy associated with precipitation and subsurface conduction because of their relatively small contributions to the overall net energy. Because energy fluxes are additive (Eq. 1), the percents can exceed 100% (see also Fitzpatrick et al., 2017). Radiative fluxes seem to dominate the melt energy in the northeastern subbasins Nubra/Shyok and Zaskar, whereas turbulent fluxes dominate in the west (Kabul/Swat/Alingar, Gilgit/Hunza) and south (Chenab1, Bavi, Sutlej). The three central basins seem to be a mix, with no clear dominant contributor to the melt signal. Because turbulent heat flux equations perform poorly for sublimation (Rupper and Roe, 2008), with many of the assumptions made in the calculation of sublimation violated in extreme terrain (Stigter et al., 2018), the flux composition of glaciers whose mass balance is dominated by sublimation-mostly concentrated on the northern boundary of Kabul/Swat/Alingar and Gilgit/Hunza and the eastern part of Nubra/Shyok-is less certain. The sublimation-dominated glaciers are a minority, though, and they do not change the broad basin-wide trends in melt energy components.
While dry, arid conditions dominate across the northeastern part of the UIB (Joshi et al., 2005;Gascoin, 2021), we do not see turbulent fluxes dominate in this area as Gascoin (2021) finds for snow because turbulent fluxes can dominate snow mass loss while radiative fluxes dominate the ice mass loss. Compared to snow, glacier ice extends to lower elevations throughout the summer (experiencing warmer temperatures and, thus, more longwave fluxes); additionally, it is darker than snow for longer periods of time and, thus, absorbs more shortwave radiation. Melt is 8.5 times more efficient than sublimation; therefore, radiatively-driven melt can dominate the mass balances on a glacier while sublimation dominates mass changes in higher elevation snowpack. Spatial patterns in melt and its components are important for understanding dominant controls and regional gradients. From the perspective of water resources, though, it is necessary to compare integrated, subbasin-wide glacier melt. Basins that have more glaciers to begin with are going to be able to produce more melt. However, the relationship of melt volume with glacierized area for each subbasin is not a simple linear relationship ( Figure 5).

Mass Balance
Although this study aimed to compute and explore patterns in glacier melt, glacier mass balance is the sole metric for which model output can be evaluated against measured values (e.g., in situ, geodetic) across this vast region. Unfortunately, in situ mass balance observations are extremely limited in the region. Still, existing ones provide some information on temporal variance and spatial distribution of mass balance on glaciers. Two of the 13,912 modeled glaciers are included in the World Glacier Monitoring Service (WGMS) database and cover shorter time periods than the model does (2001): Chhota Shigri (2003) and Patsio (2011 glaciers (WGMS, 2021).   found that these observations fall within the modeled mass balance range (see Johnson and Rupper (2020)'s Figure 3). For Patsio, there are fewer data (only 3 years), and it is difficult to make broad conclusions. Still, we do note that in those years of overlap there is a negative bias in the model. Such year-toyear comparisons have limited utility because we do not expect HAR  to capture the specific climate of any given year. Far more telling are comparisons of the mean and range: for Patsio glacier, the 3 years are roughly within the range of the modeled output (but, once again, the number of years is too limited to assess whether the model generally captures the long-term mean or variability). For Chhota Shigri glacier, the model captures some of the variance well but does not capture the most negative years. The observed mean (−0.60 m w.e.) and variance (0.45 m w.e.) are greater in magnitude than the modeled mean (−0.01 m w.e.) and variance (0.12 m w.e.) for the same period. The general biases seen in the model-observations comparisons are similar to those seen when comparing to geodetic mass balance observations. The geodetic mass balance data used in this study (Shean et al., 2020) do not directly provide information on temporal variance or mass balance distribution on a glacier, but they do provide a more spatially complete dataset for data-model comparisons. Figure 7 gives a histogram of the modeled and geodetic mass balance values (excluding the modeled values for calibrated glaciers). Salient in this histogram is the long, left tail of the modeled MBs, which appears also as the large MB bias in Figure 3; Supplementary Figure S3 (model-data MB discrepancies can be extremely large). Although, as stated in Calibration, it is not possible to attribute a model-data mismatch to the various MB components specifically, when the model yields a large model-data mismatch, it likely also simulates melt poorly. Uncertainty in melt is calculated as a function of MB bias (see Supplementary Material) to account for this, and the 11.9% of glaciers with large uncertainty (Category 3 in Figure 3; Supplementary Figure S3) have larger uncertainties in their melt and melt volumes than those with lesser MB biases.
The significantly different variances between the modeled and geodetic MB values are evident in the long left tail described above-and also in both the standard deviations and the medians' departure from the means. The model and geodetic mass balance data distributions also vary in skewness (with values of −2.20 and 0.88, respectively). These differences in variance and skewness are present at the subbasin level, as well. The glacier-by-glacier RMSE by subbasin ( Table 2) indicates that the model performs best for the Nubra/Shyok and Sutlej subbasins (both have RMSE = 0.77) and worst for Astor (RMSE = 1.69) and the combined Krishen Ganga and Jhelum (RMSE = 1.62). Mass balance does not show a clear correlation with area or other individual characteristics (e.g., aspect, slope, length, hypsometric index). The modeled MB tends to deviate from the geodetic values most for lower elevation glaciers, but we did not find a strong linear correlation between elevation and the geodetic MB-model MB difference.
Of the four lowest subbasin-wide mean modeled MBs, three are in the Karakoram. Mass balances tend to be small in the Karakoram (hence, the "Karakoram Anomaly"), but, as discussed above, a small mass balance does not imply a small melt because mass balance can still be small with a large mass turnover. The melt volumes from the three Karakoram subbasins are the greatest of all UIB subbasins. In these subbasins, glaciers are large, and there are large areas over which melt can occur. All three Karakoram subbasins show skews in the distribution of their glacier areas, with many small glaciers. A relatively small number of large glaciers contribute the majority of melt. In Table 2, it is clear that mean melt (in m w.e. yr −1 ) is not notably large in the Karakoram subbasins. It is greater than the median melt, reinforcing that the melt for the largest glaciers is higher than the mean melt for all glaciers. The large, glaciated areas more than compensate for the low subbasin-wide average melt values since mass balance in m w.e. is multiplied by area to get volumetric mass balance.

Variance
Modeled mass balance shows much more variance than the geodetic MBs for the UIB as a whole (Figure 7) and for each subbasin. Notably, there is a left skew in the model MBs that is absent in the measurements; this is reflected graphically and in the model MB mean and median values, which are lower than the corresponding values for the geodetic MB. The mode of the modeled values is greater, again reflecting the long left tail that is comprised of relatively fewer glaciers (more glaciers have greater modeled MBs). The likely culprit for the larger variance and skewness of the modeled values is the way that the calibrated PCFs are applied to all UIB glaciers (described in Calibration; see also Uncertainty below). Three factors likely impact variance when it comes to the PCF. First, we averaged local PCFs without weighting by distance or geomorphic similarity to the modeled glacier. Second, we calibrated total annual precipitation, although adjusting precipitation when it is rain vs. snow has different impacts on glacier mass balance. Finally, calibration coefficients are unlikely to be constant across time. This is especially true because the process of calibration encapsulates all biases that make the model differ from measurements-not just precipitation-into the "precipitation" correction factor. Many factors are wrapped up in the calibration and, therefore, Further, geodetic mass balance and model output mass balance may not be truly identical quantities. The geodetic mass balance data are calculated by applying a single density of ρ = 850 kg m −3 (Shean et al., 2020) to a remote-sensing-based glacier volume change, whereas the model computes mass balance through explicitly representing melt, sublimation, and accumulation and accounting for the differences in ice and snow density. The constant density assumption in geodetic MB may underrepresent variability and contribute to the lower variance in the MB observations. Another potential contributor to the variation difference is that the model does not account for SW shading or topographic effects on emitted LW. The modeled MBs in the left tail of the rightskewed distribution shown in Figure 7 are for glaciers at low elevations, and the variance of the model is greater than that of the measurements, even when excluding these outliers.

The UIB's Glaciers
Modeled melt for the glaciers in Nubra/Shyok is proportionally less than that for glaciers in Gilgit/Hunza, even though they collectively have 1.8 times the glacierized area. The Kabul/Swat/Alingar has essentially the same glacierized area as Sutlej, and yet the glaciers in the Sutlej subbasin produce 12% more modeled melt volume. Conversely, Chenab1 and Zaskar produce the same amount of modeled melt, but Zaskar's glaciated area is 15% greater. The area weighted mean melt ( Table 2), generally equal to the total melt volume in a basin normalized by the total glacier area, corroborates these statements and highlights that the Nubra/Shyok has the lowest area-normalized melt rate of any subbasin. Even when the calculated uncertainty is considered, melt does not scale directly with glacierized area.
The basins show different responses to the HAR meteorological forcing due to both spatial variations in the meteorology and the basins' considerable variations in size and the glaciers they contain. The majority of modeled melt comes from a relatively small number of glaciers, which is a result of having a range of glacier sizes in each basin. For example, 60% of the glacier melt volume in the Gilgit/Hunza comes from 17.5% of its glaciers, and 60% of the glacier melt volume in the combined Krishen Ganga and Jhelum subbasins is from 34.3% of its glaciers. 80% of the total UIB glacier melt comes from 40% of glaciers ( Figure 8).
Hypsometric curves, showing glacier area vs. elevation by 50 m elevation band (RGI Consortium, 2017), provide some insight into why some basins produce more melt water volume than others. Of the three most glaciated basins-in order, Nubra/Shyok, Gilgit/ Hunza, and Indus1-the Gilgit/Hunza glaciers melt notably more and are situated at lower elevations, where they experience warmer temperatures. Gilgit/Hunza's peak elevation in the curve of Figure 9 is 400 m lower than Indus1's and 450 m lower than Nubra/Shyok's.

Uncertainty
The 12-glacier "uncertainty test set," which was meant to capture a spread of glacier sizes and regions, allowed us to place the results into the context of the uncertainties. While a full Monte Carlo experiment, or similar, could provide greater detail, the computational requirements of the model precluded such an investigation. Although our subset of 12 glaciers is not exhaustive, it nevertheless served to highlight key sensitivities to model parameters and calibration. Although uncertainties remain large and difficult to quantify at the basin-wide/regional scale, our uncertainty analysis revealed that the SEMB parameters impose smaller uncertainties compared to the calibration (parameter-only uncertainties were found up to a MB bias of only 0.03 m w.e. yr −1 ). The "uncertainty from calibration" implicitly includes biases in the climate variables and downscaling those variables to scales requisite for distributed glacier models.  Uncertainty from calibration could have potentially been decreased by calibrating all glaciers or making the PCF calculation for an individual glacier more robust through, for example, weighting by elevation or another geomorphic characteristic. The community has highlighted difficulties of applying calibration (e.g., Rounce et al., 2020a;Rounce et al., 2020b). It is essential to strike a balance between calibrating one glacier based on others (which is common practice but introduces errors) and over-calibrating (which would negate the modeling results). The difficulty we found in applying calibration regionally, despite our attempts to circumvent this somewhat by using a PCF calculated from those of multiple nearby glaciers, is similar to previous work on single glaciers (MacDougall and Flowers, 2010;Zolles et al., 2019).
Our approach avoids over-calibrating and is computationally reasonable-but at the expense of greater certainty. While the central tendency of the modeled mass balance is similar to the geodetic mass balance observations, model-data discrepancies can be extremely large. Of the 11.9% of glaciers having a model-geodetic MB bias classified as an outlier (greater than 1.37 m w.e. yr −1 according to the outlier definition of Hubert and Vandervieren, 2008, detailed in the Supplementary Material), 5.9% have biases exceeding 2 m w.e. yr −1 and 2.47% have biases exceeding 3 m w.e. yr −1 . Some of these biases average out at the subbasin scale but, given the skew towards negative mass balance bias (Figure 7), there is likely a bias towards overestimating melt in our approach.
Uncertainties in precipitation, though partly mitigated through our choice of precipitation for a calibration parameter, may still affect the model output and results. Achieving accuracy when downscaling precipitation in mountainous regions where it varies greatly over small distances (Lutz et al., 2014) is a challenge, particularly where strong vertical gradients exist (e.g., in the UIB: Hewitt, 2014;Winiger et al., 2005). The model's precipitation gradient and PCF affect the amount of snow vs. rain. Precipitation phase directly determines surface albedo, subsurface heat flux, and sensible heat flux from rain; therefore, it directly affects a glacier's surface energy balance, which affects the MB's ablation term. Solid precipitation comprises the accumulation term of mass balance. Thus, there are two ways for a change in precipitation to impact the MB directly.
An incorrect precipitation directly impacts the energy balance and directly (and indirectly) affects the modeled MB. Importantly, perturbations of the PCF reveal that 50% of the mass balance uncertainty is due to melt (Supplementary Figure S4). The difficulty of assessing the precise share of uncertainty driven by precipitation in general provided additional motivation to apply the calibrated correction to precipitation over other variables. Although isolating uncertainty from precipitation is difficult, we do assess overall uncertainty in our results.
Model parameters may change the model's sensitivity to the PCF. We did not address this specifically and note that the computed bounds shown in Figure 5 do not include this source of uncertainty. However, since the uncertainty bounds were calculated conservatively (see Supplementary Material for calculation details), they may encompass any uncertainty added by this parameter-calibration interaction. Other sources of uncertainty include equifinality, a problem that can be mitigated by calibrating with multiple datasets (Azam et al., 2021), and glacier debris (described below). The model itself does not include glacier dynamics and neglects the radiative effects of dust and black carbon on glacier surfaces. By neglecting dynamics, our mass balance estimates do not include feedbacks associated with elevation change over time (e.g., glaciers thin into warmer climates which amplify melt processes). In addition, surface velocities would advect firn/snow from higher elevations to lower elevations which would affect the areal distribution of snow versus ice seasonally. While these types of dynamic influences are assumed to be small effects over the 13 years of focus in this study, they do introduce additional uncertainties. Coupling a dynamics model to a surface energy balance one, such as done with GloGEMFlo by Zekollari et al. (2019), would be a valuable next step-and is certainly a critical one over longer timescales.
Debris-covered glaciers are an important and complex component of glacierized regions in High Mountain Asia (HMA, e.g., Herreid and Pellicciotti, 2020;Scherler et al., 2018). This study excludes the UIB's 3,299 debris-covered glaciers, which also contribute to glacier meltwater volume. Thus, our results are limited to clean glacier contributions to meltwater in the UIB. The model also assumes the small areas of debris cover on clean glaciers (up to 15% of total glacier area) are clean ice. In regions where debris exceeds a few cm in thickness, melt rates are likely overestimated in the model. In transition regions between clean ice and debris, where debris is thinner, modeled melt rates are likely underestimated. These two effects offset each other in the modeled melt, though the magnitude of this offset depends on the areal distribution of the debris thickness and is likely variable across the region. Neglecting debris-covered portions of these clean ice glaciers, therefore, introduces additional uncertainty in our model results.
Running the model 2001-2013 but calibrating with geodetic MB data from 2000-2018 introduces uncertainty into the calibration and, thus, results. However, regional aggregations of MB over decades-long ranges points to a wide range of geodetic MBs overlapping in uncertainty bounds for differing decades of coverage (Gardelle et al., 2013;Kääb et al., 2015;Maurer et al., 2019;Shean et al., 2020). In addition, these studies show similar spatial patterns in mass balance across HMA. Therefore, although we would expect greater mass loss in more recent years, we assume that the relative magnitudes and spatial patterns of geodetic MB would remain similar between the different time periods.
Generally, the area-melt relationships hold in the face of these uncertainties for two reasons. First, because uncertainty is calculated in meters of melt, it increases with area, its effect being compounded for very large glaciers and highly glaciated basins. The basins with lower elevation glaciers (more commonly misrepresented by the model) have smaller areas, and, thus, their larger errors don't compound as much in the melt volume calculation ( Figure 5).
Second, we expect the model to have a consistent bias in overor under-predicting melt. The modeled MBs relative to observations (i.e., geodetic values) shown in Figure 7 for the full UIB have a distribution-with a long tail in negative model MB values-that is present in all subbasins. Such similarity in distributions suggests that biases are consistent between subbasins (i.e., if the model overpredicts in some basins, it likely overpredicts in all). All subbasins have the same systematic bias, and it's unlikely that the causes of uncertainty shift the model MB in one basin significantly more than in another. Therefore, although Figure 5 shows increasingly large error bounds for more glaciated basins, the relatively greater melt in the Gilgit/Hunza is not likely an artifact of the model.

CONCLUSION
This study illustrates the utility of mass and energy balance studies for providing perspective on the potential physical drivers of glacier change at regional and global scales. However, the limited nature of observations of glacier and climate variables across the UIB remains a barrier to modeling glacier melt across the region accurately. We acknowledge there are significant uncertainties associated with our model results and, therefore, treat the results as estimates of melt within a purely theoretical framework given HAR climate as the forcing. Within this context, our main findings are: A) Spatial patterns exist in the dominant contributor to glacier melt. Radiative fluxes dominate the melt signal in the northeastern subbasins (Nubra/Shyok and Zaskar), whereas turbulent fluxes dominate in the western (Kabul/ Swat/Alingar, Gilgit/Hunza) and southern (Chenab1, Bavi, Sutlej) ones. The remaining basins have no clear dominant contributor to the melt signal; melt is driven by a mix of turbulent and radiative fluxes. B) The relationship between glacierized area and melt volume is not 1:1. While basins with more glaciated area have the potential to generate more glacier melt, the extent of glaciers is not a predictor for glacier melt in the UIB. Of the three most glaciated basins (Nubra/Shyok, Gilgit/Hunza, and Indus1), the Gilgit/Hunza glaciers melt notably more (2.64-7.29 km 3 yr −1 ), at least in part due to their lower elevations. After Gilgit/Hunza, Nubra/Shyok (1.89-5.00 km 3 yr −1 ) then Indus1 (1.61-3.75 km 3 yr −1 ) produce the greatest annual volumes of glacier meltwater. C) Calibration dominates the uncertainty in our results. The selection of model parameters contributes relatively little to overall uncertainty in results. There was only a weak relationship between model performance and glacier elevation, highlighting that not every low-elevation glacier was poorly represented by the model. However, the glaciers for which the model performed poorly-calculating MBs that differed substantially from geodetic ones-tended to be situated at lower elevations.
Glaciers are crucial to maintaining and regulating the flow of the Indus River, which is "dominated by temperature-driven glacier melt during summer" (Lutz et al., 2014). Our study helps elucidate where UIB glaciers melt, a crucial prerequisite for examining changes in when they melt. The timing of glacier melt informs future hydrological projections involving the annual hydrograph (Lutz et al., 2016), extremes (Wijngaard et al., 2017), and the balance of supply with variable demands (Wijngaard et al., 2018;Ciraci et al., 2020). This research on the relative glacier melt contributions of the UIB's subbasins points to the need for a fully calibrated and validated energy balance model and climate data that can accurately capture the timing, amount, and physical drivers of melt in the region. It also suggests that a complementary understanding of the spatial patterns in water withdrawals is needed. By quantifying water use at a similarly fine spatial resolution, a complete vulnerability assessment for the UIB's glaciers could be achieved.

AUTHOR CONTRIBUTIONS
AG revised and ran the model, analyzed model output, and wrote the manuscript. SR conceived the study and provided guidance on the analysis and manuscript. Guidance using model and HAR data was provided by EJ who consulted on issues and provided feedback throughout the project. DK designed Figure 1 and provided detailed feedback and discussion on the manuscript while in preparation. RF provided project direction, research feedback, and manuscript feedback. All authors reviewed and edited the manuscript.

FUNDING
This work was supported by the NASA grants awarded to PIs Summer Rupper (NNX16AQ61G and 80NSSC20K1594) and Richard Forster (80NSSC19K1498).