Multi-decadal mass balance series of three Kyrgyz glaciers inferred from transient snowline observations

Glacier mass balance observations in the Tien Shan and Pamir mountains are sparse and often discontinuous. Nevertheless, glaciers are one of the most important components of the high-mountain cryosphere in the region; they strongly influence water availability in the arid, continental and intensely populated downstream areas. This study provides reliable and continuous mass balance series for selected glaciers located in the Tien Shan and Pamir-Alay. A combination of three independent methods was 5 used to reconstruct for the past two decades the mass balance of the three benchmark glaciers, Abramov, Golubin and No. 354. By applying different approaches, it was possible to compensate for the limitations and shortcomings of each individual method. This study proposes the use of transient snowline observations throughout the melting season obtained from satellite imagery and terrestrial automatic cameras. By combining modelling with remotely acquired information on summer snow depletion, it was possible to infer glacier mass changes for unmeasured years. Multi-annual mass changes based on high10 accuracy digital elevation models and in situ glaciological surveys were used to validate the results for the investigated glaciers. Substantial mass loss was confirmed for the three studied glaciers by all three methods, ranging from −0.30±0.19 m w.e. a−1 to −0.41±0.33 m w.e. a−1 over the 2004-2016 period. Our results indicate that integration of snowline observations into mass balance modelling significantly narrows the uncertainty ranges of the estimates, and hence highlights the potential of the methodology for application to inaccessible glaciers at larger scales for which no direct measurements are available. 15


Introduction
Glaciers are important components of the hydrological cycle in Central Asia.In this arid continental region, the intensely populated and irrigated downstream areas strongly depend on a supply of water from the cryosphere such as glaciers and snow (Kaser et al., 2010;Schaner et al., 2012;Duethmann et al., 2014;Chen et al., 2016;Huss et al., 2017).The uncertainty of water availability in the context of a changing climate creates a major potential for political tension and builds a complex set of future threats, affecting different sectors such as water management, energy production and irrigation (Varis, 2014;Munia et al., measured at an Automatic Weather Station (AWS mod ) installed in 2011 was used from 2011 to 2016 (Fig. 2a).This station is located at an elevation of 4100 m a.s.l. at a distance of about 1.5 km from the glacier terminus.ERA-interim Reanalysis data with spatial resolution of 0.78 degrees (Dee et al., 2011) were used to fill measurement gaps (Barandun et al., 2015).
The glacier mass balance was measured intensively from 1967 to 1998 (Suslov et al., 1980;Glazirin et al., 1993) and the monitoring was re-established in 2011 (Hoelzle et al., 2017).Since then, annual glaciological surveys were continuously carried out in August.For a re-analyzed and reconstructed mass balance series for Abramov from 1968 to 2014 and a detailed description of the measurement network, see Barandun et al. (2015).

Golubin Glacier
Golubin Glacier (N42 • 26.94 , E 74 • 30.10 ) is located in the Ala Archa valley in the Kyrgyz Ala-Too in the Northern Tien Shan (Fig. 1).The glacier has an area of ∼5 km 2 (as of 2016) and a north to northwestern aspect.The front is currently located 10 at an elevation of about 3400 m a.s.l. and the glacier spans to an elevation of about 4300 m a.s.l.Long-term measurements indicated an internal accumulation due to refreezing of meltwater of about +0.08 m w.e. a −1 (Aizen et al., 1997).For Golubin, mass loss was reported by Bolch (2015) and Brun et al. (2017) for recent decades.
We used meteorological data from the Alplager station located in the Ala Archa valley, situated at an elevation of 2145 m a.s.l. at a distance of about 10 km from the glacier (Fig. 1).There are several other meteorological stations in the valley, however, the Alplager station has the most complete and continuous series at high elevation covering the entire study period.
Intense glacier monitoring started in 1958 and continued until 1994 when the monitoring programme was stopped (Aizen, 1988).In summer 2010, mass balance measurements were re-initiated.Figure 2b summarizes the monitoring network at Golubin Glacier as of 2016, including a mass balance measurement network, an AWS and two terrestrial cameras installed in 2013.
A detailed description of the monitoring strategy is provided in Hoelzle et al. (2017).

Glacier No. 354
Glacier No. 354 (N 41 • 47.62 , E 78 • 9.69 ) is situated in the Akshiirak range in the Central Tien Shan (Fig. 1).The glacier covered a surface area of about 6.4 km 2 in 2016.The accumulation zone comprises three basins and the glacier tongue is oriented to the northwest.The glacier spans an elevation range of 3750-4680 m a.s.l.Mass loss since the mid-1970s was reported by different studies (Pieczonka and Bolch, 2015;Kronenberg et al., 2016;Brun et al., 2017).Summer snowfall is frequent and fresh snow can cover the entire glacier surface for several days during the melt season, significantly reducing ablation (Kronenberg et al., 2016).Evidence of superimposed ice is found and Kronenberg et al. (2016) estimated internal accumulation to be +0.04 m w.e. a −1 .An AWS (Tien Shan (Kumtor) AWS) installed at a distance of approximately 10 km to the glacier records continuous meteorological data for the study period (Fig. 1).We used daily precipitation sums and mean daily air temperature for modelling.
Since 2010, in situ mass balance has been obtained annually in late summer (Fig. 2c).Winter measurements exist for May 2014 (Kronenberg et al., 2016).A description of the meteorological input data and the mass balance measurement network, as well as a reconstruction of the mass balance series back to 2003 are provided in Kronenberg et al. (2016).

High-resolution satellite images and DEMs
To compute geodetic mass balances for Abramov Glacier, high-resolution DEMs were used based on Pléiades data acquired in 2015, and on stereo images from 2003 and 2011 from Satellite Pour l'Observation de la Terre (SPOT) 5.For Glacier No. 354, DEMs from 2003 (QuickBird) and from 2012 (GeoEye) were available from Kronenberg et al. (2016).In addition, a SPOT6 stereo-pair acquired in 2015 was used to produce an updated DEM.High-resolution images for Golubin Glacier are sparse and we had to rely on a SPOT7 tri-stereo scene from November 2014 and on Advanced Land Observing Satellite (ALOS) Prism scenes from 2006 (Table 1).For modelling purposes, the most complete and accurate DEM available of each glacier was used to represent the surface topography (Table 1).

Optical satellite and terrestrial camera images
We used freely accessible, orthorectified and georeferenced Landsat TM/ETM+ and OLIs, Terra ASTER-L1B and Sentinel-2A scenes to repeatedly observe the glacier outlines and the snowline throughout the melting season for all three glaciers.In  addition, we used the high-resolution optical satellite images as described above for snowline and glacier outline mapping.Two terrestrial cameras (Mobotix M25) overlooking Abramov Glacier were installed in August 2011.One camera was located next to the AWS mod , and the other one at approximately 500 m distance at an elevation of 4200 m a.s.l.(Fig. 2a).Due to multiple camera failures and power supply problems, pictures were lacking from the end of the ablation season in 2012 to the end of the ablation season in 2013, and again partly for the summer months in 2014.In 2015, continuous coverage was obtained from Camera 1 but only a few images could be retrieved from Camera 2 (Fig. 2a).A complete set of data was collected from both cameras for the first time in 2016.A similar setup has been installed for Glacier No. 354 in 2014, delivering continuous coverage since implementation (Fig. 2c).The camera is located at an elevation of 4145 m a.s.l.Images from the two cameras installed at Golubin were not used here due to limited image quality.Figure 3 illustrates the number of camera and satellite images with satisfying quality that were used to obtain maps of snowlines for the three glaciers.Image availability for snowline  3 Methods

Glacier outlines
Glacier extents were mapped manually based on satellite images for all three glaciers and for each year of the corresponding study period.Only cloud-and snow-free images were selected.No. 354 and Golubin were mostly debris-free.We excluded a debris-covered part, due to strongly reduced melt rates, at the western margin of Abramov Glacier (Barandun et al., 2015).

5
Annually repeated glacier front measurements using a handheld GPS for all three glaciers were included also for mapping from 2011 onwards.The same extents were used for all three methods.Errors related to glacier outlines digitized manually on remote sensing images depend on atmospheric and topographic corrections, shading, glacier surface characteristics, snow cover and local clouds, but mainly on misinterpretation of debris cover (Paul et al., 2013(Paul et al., , 2015)).Uncertainties are expected to The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-256Manuscript under review for journal The Cryosphere Discussion started: 27 November 2017 c Author(s) 2017.CC BY 4.0 License.
be in the range of ±5% of the total glacier area for low-resolution images, and smaller for high-resolution images (Paul et al., 2013).

Meteorological data
For Abramov Glacier, the air temperature data measured at the AWS was adjusted to the elevation of the former glaciological station by applying a constant lapse rate of −6 • C km −1 (Suslov et al., 1980).The ERA-interim Reanalysis dataset was adapted by applying mean monthly additive and multiplicative biases for air temperature and precipitation, respectively.The biases were calculated from long-term in situ measurements (Barandun et al., 2015).From the corrected monthly means, daily series were generated by superimposing day-to-day variability, observed at the meteorological station from 1969 to 1994.For Abramov, we generated air temperature series from 1995 to 2011 and precipitation series from 1995 to 2016.More detailed information on data preparation and their suitability is given in Barandun et al. (2015).Mean daily air temperature data measured at the Ala Archa AWS and Tien Shan (Kumtor) AWS were extrapolated to the median elevation of the corresponding glacier with monthly temperature lapse rates for the Northern and Central Tien Shan provided in Aizen et al. (1995).

Snowline delineation
A visual pre-selection of suitable camera and satellite images was taken in order to preclude problems associated with image quality such as fresh snowfall, extensive cloud cover, among others.Oblique ground-based photographs were first corrected automatically for lens distortion, then projected and orthorectified following Corripio (2004).Every pixel on the photograph was associated to the elevation of the DEM.Georeferenced products of satellite scenes were downloaded.On each camera and satellite image, the snow-covered area was digitized manually by means of visual separation of bare ice and snow (Huss et al., 2013;Barandun et al., 2015;Kronenberg et al., 2016).Manual detection allowed the observers knowledge of the snow cover depletion patterns to be integrated, and was assumed to be less error-prone than an automatic classification (Huss et al., 2013;Rabatel et al., 2013).
Errors occurred due to the pixel size of the images, slope of the terrain, the accuracy of the georeferencing and the quality of the DEM (Rabatel et al., 2012).In view of the fact that the border between ice and snow is not a clearly defined line, operator expertise is desired and beneficial.Contrast becomes rather weak, especially when the snowline rises above the firn line (Rabatel et al., 2013;Wu et al., 2014).In order to estimate the influence of ambiguous transition areas, we conducted extensive experiments on the interpretation of the surface type (see Section 4).
We assumed the spatial depletion pattern to be approximately constant in time so that camera and satellite images with minor invisible sections of the snowline due to shading, cloud cover, Landsat 7 SLC-off void-stripes or due to the terrestrial camera view angle could be included.To fill in those data gaps, we extrapolated the snowline based on information from repeated snowline observations of images with good quality over a ≈15-year period.The effect of a misinterpretation of the snowline on the calculated mass balance was investigated in detail and is described in Section 4.

Glaciological mass balance
Ablation stakes are distributed over the entire ablation zone in order to provide an optimal representation of melt patterns (Fig. 2).Each year, they are re-drilled at the initial position.An ice density of 900 kg m −3 was assigned.Snow pits were dug to the previous end-of-summer horizon to measure snow density and snow accumulation.Annual field surveys ranged from late July to late August and for logistic reasons, can vary from year to year.Winter snow measurements were carried out to retrieve a detailed snow distribution pattern, and to compute the winter balance for No. 354 and Golubin in May 2014 (Kronenberg et al., 2016).Winter surveys from 1993 and 1994 were available for Abramov (Pertziger, 1996).A model-based extrapolation of point measurements to the entire glacier surface after Huss et al. (2009) was used to retrieve glacier-wide mass balances for all years with direct measurements.The model is a combined distributed accumulation (Huss et al., 2008) and temperature index melt model (Hock, 1999) which was automatically optimized to best represent all collected point measurements from each seasonal/annual survey.

Geodetic mass balance
For Abramov Glacier, the 4-m Pléiades DEM from 1 September 2015 was used as reference.It was created using the AMES stereo-pipeline (Shean et al., 2016) and the processing parameters that were used in Marti et al. (2016).The two SPOT5 DEMs (August 2003 and November 2011) were derived from High Resolution Stereoscopic (HRS) images by the French mapping agency (Korona et al., 2009).The steps required to adjust the two SPOT5 DEMs horizontally and vertically to the Pléiades reference DEM are similar to the ones followed in a previous study on the Mont Blanc area (Berthier et al., 2014).
We created DEMs with a spatial resolution of 5 m for Golubin and Glacier No. 354 from the available two/tri-stereo pairs of high-resolution satellite imagery using standard procedures and the software PCI Geomatica (Kronenberg et al., 2016).The two/tri-stereo pairs were connected using common tie points before DEM extraction.For Glacier No. 354 a horizontal shift between the two DEMs was corrected through a DEM co-registration procedure as proposed by Nuth and Kääb (2011).For the data covering Golubin, no horizontal shift was encountered.We thus corrected only for a mean elevation difference of 2.7 m detected over stable ground.For this vertical co-registration, only terrain sections with a slope lower than approximately 30 • were selected and areas with parallax-matching problems or significant snow cover were avoided.Steep mountain walls and shading caused problems.Areas affected by these problems were manually masked out.Unmeasured areas (Abramov: 26% in 2003(Abramov: 26% in -2015(Abramov: 26% in , and 23% in 2011(Abramov: 26% in -2015;;Golubin: 30%;Gl. No. 354: 25%) were assumed to have experienced the same elevation change as the measured areas in the same altitude band and the median of the corresponding elevation bin was used for gap-filling.For elevation bins higher than 4300 m a.s.l. at Golubin (9% of total area) and for bins higher than 4500 m a.s.l.
at Gl. No. 354 (8% of total area), obvious DEM errors dominated and not enough realistic values for median elevation-change calculation were available.There, the median of the uppermost elevation band with reliable data was used to fill in the gaps.
To derive the geodetic mass balance ∆M geod , a density ρ ∆V of 850 kg m −3 was used for volume-to-mass conversion (Huss, 2013): where, A is the average glacier area and ∆t is the time in years between the corresponding image pairs.Uncertainties in the detected elevation changes and the derived geodetic mass balances are described in Section 4.

Mass balance modelling constrained by snowline observations
An accumulation and temperature index melt model closely constrained by transient snowline observations was implemented in order to infer glacier-wide mass balance.The applied methodology is a further stage in the approach presented by Huss et al. (2013).The principle of the approach is to employ the information given by the temporal change in the position of the transient snowline throughout the ablation season to constrain both the amount of winter snow accumulation and melting by iteratively calibrating a mass balance model.Thus the daily mass balance evolution through each individual year can be inferred.The approach also allows us to temporally extend mass balance estimates to the end of the hydrological year although snowline observations do not cover the entire ablation season.The methodological steps are described in more detail in the following.
A mass balance model with a spatial resolution of 20 m was driven with daily mean air temperature and precipitation sums measured at nearby meteorological stations or inferred from Reanalysis data (see Section 2.1).We used a classical temperature index melt model (e.g., Braithwaite, 1995;Hock, 2003).Melt M was calculated for each grid cell x, y and time step t based on a linear relation with positive daily mean air temperature T air (x, y, t) as Daily air temperatures are extrapolated to each grid cell using a constant temperature lapse rate based on literature values (Table 2).Different degree-day factors DDF ice/snow were chosen for snow and ice surfaces.The surface type over the glacier area was given by the snow depth updated with modelled daily snowfall and melt.The ratio between DDF ice and DDF snow , R DDF , was held constant over time.As a wide range of different ratios can be found in literature (e.g., Hock, 2003;Zhang et al., 2006;Gao et al., 2010;Huintjes et al., 2010;Liu and Liu, 2016), we decided to constrain R DDF by mass balance measurements in the Tien Shan and Pamir (Table 2).A sensitivity test shows that a variation of R DDF by ±25%, a value exceeding the maximum range found in the literature (±22%), only causes small changes in modelled mass balance (Table 3) indicating that the calibrated R DDF improves model performance to some extent but is not essential to applying the model.For more details see Section 4.
Snow accumulation C was calculated for each grid cell x, y and time step t by Table 2. Constant model parameters.The temperature lapse rate for Abramov Glacier was adopted from Barandun et al. (2015) and for Golubin Glacier and Glacier No. 354 from Aizen et al. (1995).δT /δz, δP/δz, Hcrit and RDDF are held constant throughout the modelling period.Hcrit is the elevation where precipitation is set to no longer increase linearly.The initial parameter ranges of DDFsnow and Cprec as well as the mean value obtained from annual calibration is given with its standard deviation.where P ws is the measured precipitation at the meteorological station at elevation z ws .z(x, y) is the elevation of each grid cell.The measured precipitation was extrapolated to every grid cell with a constant precipitation gradient δP/δz calculated from winter snow surveys (Table 2).Solid precipitation occurs at T air ≤ 1.5 • C with a linear transition range of ±1 • C (e.g., Hock, 1999).C prec is a scaling factor that accounts for gauge under-catch and other systematic measurement errors of precipitation (e.g., Huss et al., 2009).In order to account for smaller measurement errors during summer related to the type of precipitation (solid/liquid, wet/dry snow), C prec was reduced for the summer months to 25% of its value (Sevruk, 1981).Above a critical elevation H crit , precipitation is set to no longer increase linearly (Alpert, 1986).Our selected value of H crit approximated the elevation for which a decrease in accumulation was observed on long-term monitored glaciers situated in the Tien Shan and Caucasus (WGMS, 2013).Parameters used are summarized in Table 2.

Model calibration
We calibrated C prec and DDF snow keeping R DDF constant.C prec and DDF snow were calibrated annually and for each glacier separately, to correctly represent the winter snow accumulation and the melt rate.To calibrate C prec , we relied on the fact that,   As C prec and DDF snow depend on each other, overestimation of DDF snow could cause an overestimation of C prec , and vice versa, without affecting the modelled position of the snowline and the SCAF.To overcome this problem, the values of C prec and DDF snow were constrained to remain within realistic bounds.Through an iterative calibration procedure, we aimed at finding the best possible parameter combination without the need of any additional information.
First, we defined a plausible range for DDF snow and C prec for all three glaciers based on literature (Hock, 2003;Liu and Liu, 2016) (Table 2).For each DDF snow , an optimal parameter C prec was calibrated through iteratively narrowing a plausible range of initial values of C prec (Fig. 4).In this way, the RMSE between the directly modelled winter snow accumulation C modelled and the modelled cumulative melt from the onset of the ablation season to each observation date M modelled was minimized until no further improvement of the RMSE was observed.
Second, the performance of each DDF snow with its optimal C prec pair, was evaluated based on the RMSE SCAF between the observed SCAF obs and the modelled SCAF modelled for all available snowline observations within one year (Fig. 4).The range of DDF snow was narrowed around the best solution and the optimization process was restarted until no further significant improvement of the RMSE SCAF was observed.The calibration procedure was repeated for each year individually.
A minimum of two images was needed to enable application of our calibration approach.The influence of the image frequency and distribution was assessed in detail with sensitivity experiments described in Section 4. In a last step, the calibrated model was re-run with the ideal parameter set.This snowline-constrained mass balance model was thus applied to derive con- tinuous daily mass balance series that agreed with the snow depletion patterns observed by remote sensing imagery.In the following, we refer to the results obtained by the methodology described above as snowline-derived mass balance.

Adjustments to enable comparison of different methods
Geodetic surveys provide an estimate of the total mass change of a glacier ∆M geod , whereas snowline-derived mass balance series refer to the surface balance B a and do not account for internal and basal components of the mass balance (Cogley et al., 2011).For comparing the results, we adjusted the glaciological and snowline-derived surface mass balances with an estimate of the internal/basal mass balance of +0.07 m w.e. a −1 for Abramov (Barandun et al., 2015), +0.08 m w.e. a −1 for Golubin (Aizen et al., 1997) and +0.04 m w.e. a −1 for No. 354 (Kronenberg et al., 2016).Values are positive due to a significant amount of refreezing of meltwater in cold firn.
To compare the results of the different methods, the time periods covered by the datasets also needed to be homogenized.We thus adjusted the observation period of the snowline-derived mass balance to exactly match the respective periods of geodetic and glaciological mass balance.However, the final results B a(fix) derived from snowlines are presented for the fixed dates of the hydrological year (1 October to 30 September) and only include internal/basal mass balance.
4 Uncertainties and model sensitivity

Glaciological mass balance
Uncertainty σ glac related to the direct glaciological measurements for Abramov was adopted from Barandun et al. (2015), and for Glacier No. 354 from Kronenberg et al. (2016).Uncertainties concerning the glaciological mass balance of Golubin were calculated after Kronenberg et al. (2016).Uncertainties regarding all three glaciers are summarized in Table 3 and range between 0.24 and 0.30 m w.e. a −1 .

Geodetic mass balance
The total uncertainty of the geodetic mass balance estimate includes a random and systematic error.We followed Brun et al.
(2017) for the calculation of the random error on the geodetic mass balance estimate but did not assess the systematic error.
Uncertainties in elevation differences were quantified by computing the area-weighted mean of the absolute difference off- satellite data to elevation differences measured in situ (Berthier et al., 2014).For Golubin, a value of 1.8 m indicates a slightly lower DEM quality.The uncertainty related to the density assumption for converting volume to mass change was assumed to be ±60 kg m −3 for time intervals larger than five years (Huss, 2013).For shorter periods, we used a more conservative estimate of ±120 kg m −3 .The elevation uncertainty for unmeasured glacier zones was roughly estimated to be five times as large than 13 The Cryosphere Discuss.(Table 3).The individual components were estimated as follows: (1) The accuracy of the mapped SCAF is dependent on the positioning and the transect of the snowline, the georeferencing of the images, and the extrapolation of the snowline to invisible areas (Huss et al., 2013).The limit between snow-and ice-covered areas is often not a clear line but rather a transition zone, especially for glaciers with superimposed ice.To account for the total uncertainty related to the mapping procedure, we identified an upper-and lowermost position of the surface that could be classified as either snow or ice on each available image.Hence this zone included all ambiguous areas observed, such as cloudcovered regions, shading, superimposed ice or invisibility due to reduced image quality (e.g., Landsat 7 SLC-off void-stripes, invisible areas on photographs).We interpreted the zone to be either entirely snow-covered or entirely snow-free.The standard deviation of the minimal, maximal and optimal SCAF was used as an uncertainty.This uncertainty was calculated for each image individually.To evaluate the corresponding effects on calculated mass balance, the model was re-run with the maximal and the minimal possible SCAF.The standard deviation of the mass balance, σ map , ranged between 0.06 to 0.09 m w.e. a −1 for the three glaciers.
(2) To estimate the effect of varying image availability, we repeated the modelling using different snowline observation frequencies and temporal distributions throughout the summer for calibration.Due to limited image availability, this could only be conducted for the few years when many images were available (Fig. 3).We used the results to create a look-up table that linked the image frequency, the distribution over the ablation season and the last observation date of the season to an uncertainty estimate in the calculated annual mass balance, σ dis .Tests showed that the model reacts more sensitively to the image distribution than to reduced image frequency (Fig. 5).A minimum of two images well distributed throughout the ablation season (i.e., at the beginning/middle and at the end) is sufficient to achieve reliable mass balance estimates.Greater uncertainties were found if images were concentrated on for example, a few days at the beginning of the ablation season (Fig. 5).In this case, image frequency cannot compensate for the missing information on the snow depletion pattern.An image taken towards the end of the ablation season is more important than images from the beginning of the summer.Our assessment of image availability resulted in smaller uncertainties for Abramov (σ dis = 0.09 m w.e. a −1 ) than for Golubin Glacier (σ dis = 0.16 m w.e. a −1 ) and for Glacier No. 354 (σ dis = 0.18 m w.e. a −1 ) (Table 3).
(3) To estimate the uncertainty caused by the DEM used for the modelling, we compared our results to those obtained from model runs that used lower-resolution DEMs.For this experiment, we replaced the high-resolution DEM with the SRTM DEM.This enabled us to both investigate the sensitivity of the results to DEM quality and to assess our assumption of the for Golubin Glacier and of 0.14 m w.e. a −1 for Glacier No. 354.These results demonstrate a relatively low sensitivity of the presented model to meteorological input data and underline the potential of our methodology for regional application based on minimal input data.
(5) To test the uncertainty introduced by the constant (i.e.uncalibrated) model parameters, δT /δz, δP /δz and the R DDF were varied by ±25% for each glacier and year.A mean standard deviation, σ para , of around 0.17 m w.e. a −1 was found.We 10 additionally tested the behaviour of the model relative to the individual parameters and identified a higher sensitivity to δT /δz and R DDF , whereas sensitivities to the other parameters were minor (Table 4).
Components 1 to 5 are assumed to be independent of each other and are combined as RSS to represent the total error of the snowline-derived annual mass balance σ snl .

Long-term snowline-derived mass balances
We found that the mass balance model constrained by snowline observations is capable of representing the observed SCAFs on satellite and terrestrial camera images within ±8% for Abramov, ±12% for Golubin and ±7% for Glacier No. 354.Comparing the SCAF observed on camera and on spaceborne images for the same day reveals a RMSE of 2.5%.Tests showed that the 5 influence of the image source (terrestrial/space-borne) on the inferred mass balance is negligible.
Annual glacier-wide, snowline-derived surface mass balance calculated for the hydrological years 1998 to 2016 for the three benchmark glaciers located in the Tien Shan and Pamir-Alay are predominantly negative over the two last decades (Figure 8 and in Table 5).Study periods depend on the data availability for each glacier.  of September, and summer snowfalls during this month were rare (Aizen et al., 1995), likely leading to continued mass loss.
Data availability, however, was rather critical for Golubin in 2008, which was also reflected by the stronger uncertainties of the annual mass balance.The last snowline observation dates from as early as the end of July.

Comparison to glaciological and geodetic mass balances
The glaciological and geodetic surveys delivered two extensive and independent datasets for validation of the snowline-derived mass balance series.Joint analysis of the data sets permitted robust conclusions to be drawn about the mass change and its temporal dynamics over the past two decades.The glaciological mass balance measurements showed a good agreement with the mass balance inferred using snowline observations for the same time periods (Table 6).Annual glaciological mass balances were reproduced with an RMSE of less than ±0.26 m w.e. a −1 for all three glaciers using the snowline approach (Fig. 6).For Golubin and Glacier No. 354, the glaciological mass balances was slightly more negative than the snowline-derived balances (Table 6).For Abramov, on the other hand, the glaciological mass balance was somewhat less negative than the snowlinederived results.In general, a satisfactory agreement was obtained between the two methods for all three glaciers (Fig. 6).
Squared correlation coefficients between snowline-derived and glaciological mass balances between R 2 =0.63 (Abramov) and R 2 =0.90 (Golubin) were found.Table 7 and Figure 7 summarize the results obtained from the different geodetic surveys.For the comparison with the geodetic mass change, an estimate for internal/basal mass balance was added to the snowline-derived surface mass balance (see Subsection 3.6.2) and referred to as the snowline-derived total mass change.The geodetic method reveals a mass balance of −0.22 ± 0.42 m w.e. a −1 for Golubin Glacier from 8 September 2006 to 1 November 2014 (Fig. 7).The corresponding total snowline-derived mass balance was more negative with −0.38 ± 0.35 m w.e. a −1 for the same period.However, the The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-256Manuscript under review for journal The Cryosphere Discussion started: 27 November 2017 c Author(s) 2017.CC BY 4.0 License.Table 6.Annual surface mass balance B a(meas) for the measurement periods (i.e., exact dates of the surveys) based on direct glaciological surveys and on the snowline approach for the three glaciers in m w.e a −1 .results of the two methods agree within their error bars (Table 7 and Fig  6 Discussion

More accurate modelling through integrating snowline observations
In order to demonstrate the advantage of integrating snowline observations on repeated remote sensing data throughout the melting season for increasing the confidence in mass balance modelling, we ran the same accumulation and temperatureindex model without the use of snowline or any other direct observations for calibration for all three glaciers from 2004 5 to 2016 (See Section 3.6).The same constant parameters were used (Table 2).We chose C prec to account for a 20%measurement error of the recorded precipitation (Sevruk, 1981) DDF snow (5.5 mm day −1 • C −1 ) as recommended by Hock (2003) for the former Soviet territory, for all three glaciers.The parameters were held constant over time.

Intercomparison of methods to determine glacier mass balance
A satisfying agreement was found between all three independent methods used to compute glacier-wide mass balance for the three benchmark glaciers in the Central Tien Shan and Pamir-Alay for the past two decades.In the following we discuss shortcomings and advantages of each method and point out the limitations of the individual approaches.The snowline model represents the direct measurements well and performs satisfactory for all three glaciers (Table 7, 5 and Fig. 11).For Glacier No. 354, a larger misfit between the glaciological balance and the snowline model was found, especially for 2011.A possible reason might be the limited stake network at the initiation of the monitoring programme, large errors of the stake readings and an unknown measurement date for some ablation stakes, affecting the calculation of the glaciological mass balance for the considered year.Discrepancies between the two approaches are also found at Golubin for 2016 and at Abramov for 2015.For Abramov, summer snowfalls in mid-August 2015 stopped the ablation season early.This pattern is mirrored in the daily mass balance derived from the snowline approach.For Golubin, the last snowline observation is from the end of August and matches the field observations.No clear indication of a poor performance of the snowline approach could thus be identified for both of the glaciers for the considered years.
An important problem is the varying measurement periods of the glaciological mass balances for the selected glaciers.Due to changing period lengths, the data do not always represent a complete mass-balance year, and might thus not be representative, making comparison of the results with other methods, glaciers and regions difficult.Interpretation of the results, their contextualization and application in other study fields, such as in hydrology or climatology, are also hampered through the varying and irregular investigation periods.Based on our methodology, we are now able to derive homogenous glacier mass balances for comparable periods of the hydrological year.
An important factor limiting the applicability of mass balance modelling constrained by snowline observations is the dependence on good satellite imagery to map the snowline throughout the ablation season.However, the sensitivity analysis (Section 4) shows that a minimum of only two images that are well distributed throughout the ablation season are sufficient to retrieve reliable results.Image availability is most important close to the end of the ablation season.Taking into account the increasing number of satellite sensors that provide a range of possibilities to observe snowlines in the future will partly resolve this limitation.By comparing the snowline-derived daily mass balance for the years in which seasonal in situ measurements are available, we were able to investigate the effect on the results of the two parameters used for model calibration.Figure 11b shows that the mass balance of Golubin is slightly underestimated at the beginning of the ablation season, and hence the modelled melt is also too low.This shortcoming of the calibration procedure cannot be overcome without including additional data.
The geodetic mass balance and the snowline-derived results agree well, in particular for the recent years (Fig. 9).Overall, a slightly greater mass loss is calculated for Abramov and Glacier No. 354 using DEM differencing.Especially during the earlier part of our study period, the mass balance inferred with the snowline approach seems to be too positive.Limitations related to the geodetic approach are mainly connected to the limited stereo acquisitions in the first years of the 21st century.
In recent years, image availability strongly increased, but it is still not common to find a suitable scene from the end of the hydrological year for a selected glacier or region with sufficient quality for a sound geodetic evaluation.Fresh snowfall or low image contrast (in particular in the accumulation areas) interfere with the DEM quality but cannot be avoided and have thus to be corrected for, increasing the uncertainty of the result.We present geodetic mass balances for periods shorter than five years, a critical time interval for an accurate volume-to-mass conversion.Huss (2013) shows a high variability of the volumeto-mass conversion factor for short observation periods (≤ 3 years), especially for glaciers with close-to-zero mass balances in combination with strongly varying mass balance gradients.For the observation periods considered in this study, annual mass balances were predominantly negative but moderate variations of the mass balance gradients have been observed (Barandun et al., 2015;Kronenberg et al., 2016).We identified a rather large elevation change for the short observation periods, and are thus confident that the chosen conversion factor lies within the uncertainty range assigned here (see Section 4).
The glaciological and the snowline-derived mass balances refer to surface balance components only.The geodetic mass balance, on the other hand, takes into account the total mass change, thus including internal and basal ablation and accumulation.This is a limiting factor for direct comparison.Evidence of refreezing meltwater in cold firn is reported for all three glaciers (Suslov et al., 1980;Aizen et al., 1997;Dyurgerov and Mikhalenko, 1995) and can have a significant effect on the total mass change.The values used in this study to account for internal and basal mass balance are first-order approximations which improve the comparability between the different methods.However, the uncertainties in these estimates are considerable.

Comparison to other studies
We performed a comprehensive comparison of long-term averages of mass balance derived from the snowline approach to independent studies based on geodetic surveys using different sensors and modelling, both for the investigated glaciers and for regional mass balances (Fig. 12).Note that the study periods vary between the different studies and results might thus not be directly comparable.
For Abramov Glacier, we find mass balances in between the results derived by Gardelle et al. (2013) and Brun et al. (2017) based on the comparison of DEMs, overlapping within the respective uncertainty ranges.Mass changes reported by Gardelle et al. (2013) are most likely too positive as SRTM C-Band penetration depth into snow (Kääb et al., 2015;Berthier et al., 2016) might have been underestimated for the cold and dry snow of accumulation areas (Dehecq et al., 2016).Furthermore, we also compared our results for the investigated glaciers to region-wide assessments in order to investigate their regional representativeness (Fig. 12d-f 2017) but indicates smaller mass losses than other large-scale studies (Fig. 12e).Close agreement between the different regional studies is found for the Central/Inner Tien Shan where Glacier No. 354 is located (Fig. 12f).

Conclusions
In this study we used three independent methods to reconstruct robust mass balance series at high temporal resolution for Abramov, Golubin and No. 354 located in the Pamir-Alay and Tien Shan mountains for the past two decades -a period for which only little is known about glacier behaviour.We proposed a methodology to derive glacier surface mass balances series for unmeasured glaciers based on mass balance modelling constrained by repeated snowline observations.We recommended including snowline observations in the glacier monitoring strategy to reduce uncertainty and to increase robustness of the data.We used extensive geodetic and glaciological surveys to validate the results, and found satisfying agreement between the independent methods.Our snowline approach reproduced observed annual to decadal mass balances satisfactorily for all three glaciers, and enabled the calculation of daily mass balances for arbitrary periods, and is, hence, capable of covering the entire hydrological year based on minimal data input.Some of the shortcomings of the glaciological and geodetic surveys could thus be overcome.
The results of all three methods confirm a continuous mass loss of the three benchmark glaciers Abramov, Golubin and No.

Figure 1 .
Figure 1.Overview map of Central Asia showing the location of glaciers (blue) with available long-term mass balance measurements.Glaciers investigated in this study are marked in red.The insets show the position of the automatic weather station (AWS) for (a) Golubin, (b) Abramov and (c) Glacier No. 354.Note that AWS mod.indicates the new AWS installed in 2011 and AWS hist. the location of the old glaciological station at Abramov.

10Figure 3 .
Figure 3. Image availability and distribution for snowline mapping.Numbers indicate the total available scenes per year and glacier.Prior to 1998, image coverage is sparse for all three glaciers.
54 ±0.74 5.49 ±0.46 3.04 ±0.66 mm day −1 • C at the position of the transient snowline, icemelt had not yet started but all winter snow was melted.The modelled cumulative melt, calculated at the position of the observed snowline, is thus interpreted as the total amount of accumulated winter snow that melted from the onset of the ablation season until the snowline observation date.Using the melt model, we can infer the winter accumulation at the beginning of the ablation season along each observed snowline.This quantity needs to agree with The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-256Manuscript under review for journal The Cryosphere Discussion started: 27 November 2017 c Author(s) 2017.CC BY 4.0 License. the directly modelled snow accumulation at the end of the winter season.DDF snow was calibrated to best represent all SCAF observations of one ablation season (Fig.4).

Figure 4 .
Figure 4. Calibration procedure to obtain an ideal combination of DDFsnow and Cprec.An initial range for DDFsnow and Cprec wasnarrowed down through comparison to snowline observations until an optimal solution for both parameters was found.First, for each initial value of DDFsnow, the best value of Cprec was determined constraining the modelled cumulative melt (M modelled ) at the snowline position to agree with the modelled winter snow accumulation (C modelled ) for the same location.Second, the performance of each DDFsnow was evaluated to narrow down the range of DDFsnow by comparing modelled and observed snow-covered area fractions (SCAF).This was repeated until an optimal solution is reached.
glacier in 50 m altitude bins.The resulting values for Abramov, 1.0 m for 2003-2015 and 0.6 m for 2011-2015, and for No. 354, 0.7 m for 2012-2015, are in line with the uncertainty of 1.3 m found over the Mont Blanc area through comparison of similar , https://doi.org/10.5194/tc-2017-256Manuscript under review for journal The Cryosphere Discussion started: 27 November 2017 c Author(s) 2017.CC BY 4.0 License. the uncertainty determined for measured locations.We assumed independence between the different error components and combined them as Root-Sum-Square (RSS) to the total uncertainty for the geodetic mass balance, σ geod .4.3 Snowline-derived mass balance The uncertainty introduced by the mass balance model constrained by transient snowline observations σ snl depends on (1) the delineation accuracy of the SCAF, σ map , (2) the image frequency and distribution throughout the ablation season, σ dis , (3) the DEM quality, σ DEM , (4) the meteorological input data, σ meteo , and (5) the uncertainty in constant model parameters, σ para

The
Cryosphere Discuss., https://doi.org/10.5194/tc-2017-256Manuscript under review for journal The Cryosphere Discussion started: 27 November 2017 c Author(s) 2017.CC BY 4.0 License.unchanged topography for the study period.The effects of a reduced DEM quality for all three glaciers are small (σ DEM < 0.03 m w.e. a −1 ).(4) We investigated the uncertainty related to the meteorological input data, σ meteo , by re-running the model with the climatological average daily temperature and precipitation series for each glacier instead of the actual meteorological series.The test assessment revealed an RMSE of 0.13 m w.e. a −1 for the annual mass balance of Abramov Glacier, of 0.23 m w.e. a −1 SCAF (few observations, well distributed) modelled SCAF (few observations, poorly distributed) modelled mb (few observations, poorly distributed) modelled mb (few observations, well distributed)

Figure 5 .
Figure 5. Example of the (a) daily SCAF and (b) cumulative daily mass balance (modelled mb) obtained with different sets of image frequency and distribution for Abramov Glacier in 2016.Snowline observation dates used to calibrate the model are indicated with symbols.The modelled daily SCAF and cumulative daily mass balance and their corresponding snowline observations are shown with the same colour.Three images at the beginning of the ablation season (blue), three images well distributed throughout the ablation season (green), and all available images (purple) were used.
Abramov exhibited a mean annual mass balance of −0.30 ±0.19 m w.e. a −1 from 2004 to 2016.For Golubin and Glacier No. 354 a slightly more negative annual average balance 10 of −0.41 ±0.33 m w.e. a −1 and −0.36 ±0.32 m w.e. a −1 , respectively, was calculated for the same time period (Table 5).No clear mass-balance trend for the three glaciers was identified over the investigated periods.Two phases of close-to-zero mass The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-256Manuscript under review for journal The Cryosphere Discussion started: 27 November 2017 c Author(s) 2017.CC BY 4.0 License.
. 9).Comparison of digital elevation models indicated that Glacier No. 354 had a mass balance of −0.58 ±0.31 m w.e. a −1 from 27 July 2012 to 1 October 2015 (Fig.7) and −0.42 ±0.07 m w.e. a −1 from 1 September 2003 to 27 July 2012(Kronenberg et al., 2016).The total annual mass change for the same time intervals derived from the snowline approach was −0.53 ±0.43 m w.e. a −1 and −0.25 ±0.31 m w.e. a −1 , respectively.For the first period, the results are in good agreement, whereas for the second period the mass balance model constrained by snowline observations indicates a significantly less negative mass balance (Fig.9).For Abramov, a geodetic mass balance of −0.39 ±0.16 m w.e. a −1 from 27 August 2003 to 1 September 2015 and of −0.36 ±0.26 m w.e. a −1 from 29 November 2011 to 1 September 2015 was calculated.For the same periods the snowline model reveals a total mass change of −0.25 ±0.20 m w.e. a −1 and of −0.43 ±0.17 m w.e. a −1 , respectively.For the first period, the snowline-derived mass balance indicates a less negative mass balance.The second period is in good agreement (Table7and Fig.9).For all three glaciers and The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-256Manuscript under review for journal The Cryosphere Discussion started: 27 November 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 9 .
Figure 9. Cumulative snowline-derived mass change (red) in comparison to the geodetic mass change (circles) for (a) Abramov, (b) Golubin and (c) Glacier No. 354.The shading indicates the uncertainty range of the snowline-derived mass change.

Figure 10 .
Figure 10.Comparison of the cumulative mass balance derived from unconstrained mass balance modelling to snowline-derived mass balance series from 2004 to 2016.

TheFigure 11 .
Figure 11.Cumulative daily glacier-wide mass balance inferred from the snowline approach for (a) Abramov, (b) Golubin and (c) No. 354 for the year 2014 (blue line).The grey lines indicate the spread obtained by using different constant parameters to run the model (Section 4) and the red cross indicates the measured glaciological balance both for the winter and the annual period.
The average mass balance for Abramov of −0.38 ±0.10 m w.e. a −1 (2002-2014) derived by Brun et al. (2017) using multi-temporal ASTER DEMs indicates a stronger mass loss than our study.We note, however, that the start and end dates of their geodetic mass The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-256Manuscript under review for journal The Cryosphere Discussion started: 27 November 2017 c Author(s) 2017.CC BY 4.0 License.balance represents a mean over a mosaic of different dates, thus hampering the direct comparison to our results; in addition the differences are still within their error bounds.For Golubin Glacier, the inferred mass balance is in close agreement with the geodetic mass change reported by Bolch (2015) (Fig.12b).Brun et al. (2017) computed a mass balance of −0.04 ±0.19 m w.e. a −1 for ≈2002-2013, which is consistently less negative than our estimate but still lies within the respective error bounds.For Glacier No. 354, an excellent agreement between all available mass balance assessments was found (Fig.12c).Brun et al. (2017) reported a mass balance of −0.46 ±0.19 m w.e. a −1 for ≈2002-2014.
).Brun et al. (2017) divided the Pamir-Alay and the Pamir into two different regions, whereasGardner et al. (2013),Gardelle et al. (2013),Kääb et al. (2015) andFarinotti et al. (2015) did not make this distinction.For the Pamir, widely varying mass balance estimates were presented by the different studies which might be related to the important methodological differences and inconsistent time periods considered.Our results for Abramov are close to the average of the regional studies (Fig.12d).The interannual variability and in particular a very negative mass balance for 2008 and a positive balance in 2010 for Abramov found in the present study agrees well with modelled mass balance series reported byPohl et al. (2017) for the Pamir from 2002 to 2013.The snowline-derived mass balance for Golubin agrees with the region-wide estimates byBrun et al. ( 30 ±0.19 m w.e. a −1 located in the Pamir-Alay than for the Tien Shan glaciers Golubin of −0.41 ±0.33 m w.e. a −1 and Glacier No. 354 of −0.36 ±0.32 m w.e. a −1 from 2004 to 24 The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-256Manuscript under review for journal The Cryosphere Discussion started: 27 November 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 12 .
Figure 12.Long-term average mass balances from various studies in comparison with the snowline-derived mass balance for (a) Abramov, (b) Golubin, (c) Gl.No. 354.Results from the snowline approach for the three glaciers are compared to region-wide mass balance estimates for (d) the Pamir, for (e) Tien Shan (North/West) and for (f) the Tien Shan (Inner/Central).Lines indicate the mean annual mass balance over the respective periods for each study, and shaded squares represent the corresponding uncertainty.For the study by Brun et al. (2017) results based on ASTER DEMs and IceSat are presented separately.

Table 1 .
Available data on glacier monitoring for the three glaciers used in this study.The dates marked with an asterisk indicate the DEMs used for the modelling.

Table 3 .
Overall average uncertainties, as well as uncertainties specified for each component related to the geodetic and snowline-derived mass balances for all three investigated glaciers in m w.e. a −1 .

Table 4 .
Model sensitivity to the different constant input parameters for each glacier in m w.e. a −1 .See text for details on the experiments.

Table 5 .
Annual surface mass balances for the three glaciers for the hydrological year B a(fix) in m w.e a −1 derived from the mass balance model constrained by snowline observations.At the bottom of the table the standard deviation (STD) of the mass balance for each glacier from 2004 to 2016 is given.recognizedfor2002-2004 and for 2009-2011 for all glaciers.Glacier No. 354, situated in a more continental climate regime, shows the weakest interannual variability and has a positive balance only in 2009 (Table5).For Golubin, on the other hand, most negative values were found in 2012 and 2014.Abramov and Golubin had strongly negative balances in 2006 and 2008.For Abramov, the snowline observations indicated that the snowline rose close to the upper edge of the glacier already at the end of August.For Golubin, observations of the last image of the season showed that the SCAF decreased to

Table 7 .
Geodetic mass balance ∆M geod(meas) and the total annual mass change derived from the snowline approach ∆M snl(meas) for the three glaciers and for the periods corresponding to the geodetic survey in m w.e a −1 .