The Bordu, Kara-Batkak and Sary-Tor glaciers are located in the Inner Tien Shan, south of lake Issyk-Kul, in Kyrgyzstan (Figure 1). While Bordu and Sary-Tor glaciers are situated on the north-western slopes of the Ak-Shyirak massif, Kara-Batkak glacier is located in the Terskey Ala-Too range (Figure 1). More specifically, Bordu and Sary-Tor glaciers are located close to the Kumtor Gold Mine in the valleys east of Glacier No. 354 (Kronenberg et al., 2016). Both glaciers have an elevation range between 3,900 and 4,800 m above sea level (a.s.l.) and consist of one large glacier tongue with smaller hanging tributaries. For Sary-Tor glacier, the hanging glacier on the eastern side is not considered in this study (Figure 1).

Ice thickness measurements on Sary-Tor glacier were performed revealing a maximum ice thickness of 159 m and an ice volume of 0.135 km³ (Petrakov et al., 2014). Based on the reflectivities from the radar measurements, the Sary-Tor glacier is assumed to be a polythermal glacier containing temperate ice over the bed (Petrakov et al., 2014). Regarding Bordu glacier, ice thickness measurements in 2017 showed ice thicknesses of up to 148 m with an estimated ice volume of 0.291 km³ (Van Tricht et al., 2021). The Kara-Batkak glacier (spanning 3,300–4,400 m a.s.l.) consists of a terminus part and several flat areas separated by steep icefalls, while the accumulation area is divided into three separate catchments. Ice thickness measurements on this glacier revealed a maximum ice thickness of 114 m and an ice volume of 0.096 km³ in 2017 (Van Tricht et al., 2021).

The three selected glaciers have been subject to different glaciological measurements. Between 1984 and 1989, mass balance (MB) measurements were performed on Sary-Tor glacier and used to estimate the total balance in the entire Ak-Shyirak (Ushnurtsev, 1991; Dyurgerov et al., 1995). However, in 1990, all measurements were interrupted. In 2014, detailed MB measurements were reinitiated by the Tien Shan High Mountain Scientific Research Center as part of the CHARIS Project (Satylkanov, 2018). MB measurements on Bordu glacier were started in 2015 and are continued today. The Kara-Batkak glacier is one of the few glaciers in the area with a relatively large amount of data with a mean specific mass balance series dating back to 1957 (Dyurgerov et al., 1995; WGMS, 2020). However, between 1998 and 2013, glaciological measurements were interrupted. For the 1957–1998 period, MB calculations on Kara-Batkak glacier were made based on a vertical precipitation gradient for accumulation and a combination of stake measurements on the glacier snout as well as measurements of liquid runoff at a gauge at the meltwater lake near the glacier terminus. Hence, the mass balance data of 1957–1998 are associated with a larger uncertainty compared to direct measurements. In 2013, glaciological measurements were reinitiated (Satylkanov, 2018).

For all three glaciers, the present-day annual measurements consist of density, ablation stake and snow pit measurements at the end of the summer season, as well as measurements of snow depth and density from snow pits for the winter balance at the end of spring (typically in May). The latter generally shows the amount of snow between the end of summer and the end of spring as little or no melt occurs on the glaciers during winter months. Although in theory only the surface mass balance is measured in this way, we consider it as mass balance values, since basal and internal mass balance are assumed to be generally much smaller compared to the surface mass balance. All the point measurements are then extrapolated to mean specific mass balance values for 100 m elevation bands taking into account the area of the elevation bands (Figure 2). Eventually, the mean specific mass balance for the entire glacier is calculated as the weighted mean of all these banded values, using their areas as weights (Figure 3). One of the main error sources in these data is due to the calculation of the mass balance in the upper bands (> 4,200 m a.s.l. for Kara-Batkak glacier and >4,400 m a.s.l. for Bordu and Sary-Tor glaciers) as these elevations are often unreachable due to crevasses and the risk for avalanches. Hence, accumulation, ablation and mass balance in these areas are typically assessed, based on the available data from other bands, sparse measurements, from transient snowline dynamics, and from temporal trends of the glaciers as a whole.

The climate regime in the Inner Tien Shan is characterised by a strong seasonal variation in precipitation and low air temperatures. Precipitation increases during spring and reaches a maximum during summer. Between May and September, about 75% of the annual precipitation occurs (Aizen, 1995; Kutuzov and Shahgedanova, 2009). During winter, the Siberian high over the area brings cold and dry air masses with very limited precipitation (Aizen, 1995). Therefore, accumulation peaks in spring and prolongs to early summer, when the precipitation increases, and temperatures remain low. In summer, air temperatures rise above zero creating unfavorable conditions for snow accumulation. However, the frequent (snow) showers observed in summer temporarily create a snow layer reducing the melt (in particular for the higher areas) substantially. Hence, the selected glaciers are considered examples of spring-accumulation type glaciers (Yang et al., 2013). This type of glaciers has a significantly larger sensitivity to warming compared to winter-type glaciers. Warming on spring- and summer-accumulation type glaciers prolongs the melting season, causes a significant reduction in snow accumulation due to an increased fraction of liquid precipitation, and decreases the surface albedo (Fujita and Ageta, 2000).

The surface energy mass balance model applied in this study is driven by temperature and precipitation data at hourly resolution. To reconstruct a continuous historical dataset covering 1750–2020 for the location of the glaciers, several temperature and precipitation datasets are combined by establishing correlations between overlapping periods. This reconstruction is described in the sections below. The locations of the consulted meteorological stations are shown in Figure 1. Two separate datasets are prepared. One for the Kumtor-Tien Shan meteorological station (used for Bordu and Sary-Tor glaciers) and one for Chon-Kyzyl-Suu meteorological station (used for the Kara-Batkak glacier) (see Figure 4). The used data is summarised in Table 1. Data from the local glacier meteorological stations which were installed recently, are used to find the lapse rate and atmospheric transmissivity which is discussed further below.

The causes and exact timing of the Little Ice Age (LIA) and the maximum extent of glaciers in Kyrgyzstan remain uncertain. In the Tien Shan, the LIA moraines are believed to be typically located a few hundred metres ahead of the glacier terminus and characterised by little to no vegetation cover, massive piling of loose tills and sharp-crested ridges (Solomina et al., 2004; Li et al., 2017). Often, several moraine ridges belong to the LIA, indicating different stages of glacial stagnation within this period.

Obtained lichenometric dates of the LIA moraines range between the 17^th to mid-19^th century and can be roughly divided into three advancing phases, around 1400–1500, 1700–1750 and 1800–1850 (Savoskul and Solomina, 1996; Solomina et al., 2004). Later, some smaller advances occurred again around 1870–1890, 1910–1930 and 1940–1950. These periods correspond nicely with the meteorological data showing a dip in temperature (negative anomaly) for these periods (Figures 5B, 6B).

In the Terskey Ala-Too, the mountain range in which Kara-Batkak glacier is located, the majority of terminal and lateral moraines formed during the (mid-)19^th century and the latest significant advance was dated to the 1910s (Solomina et al., 2004; Kutuzov and Shahgedanova, 2009). Dating of moraines suggests that concerning Kara-Batkak, the maximum LIA extent likely occurred in the first half of the 19^th century (Solomina et al., 2004; Li et al., 2017). This seems to be in line with the reconstructed temperature dataset which shows a last major dip around 1850 and two smaller dips around 1910 and 1940 (last periods of advance) (Figures 5B, 6B).

A two-dimensional time-dependent surface energy mass balance model is used to model the mass balance in a typical fixed balance year (1 October–30 September). The main equation that is solved is the mass balance equation which corresponds to the sum of mass gain and mass loss for each grid cell (25 m × 25 m) at every time step (in this study every hour): B t + 1 =  B t +  P s , t  + ( M t − R w ) (1)

As shown by Eq. 1, the mass balance at every new time step B_t+1 (m w.e.) is the sum of the mass balance at the previous timestep (B_t), the accumulation or solid precipitation at that timestep (P_s,t) (m w.e.) and the ablation or melt (M_t) (m w.e.), which is determined by the net energy flux per unit area at the glacier surface Ψ (W m⁻²) (see Energy Balance Equation). Regarding ablation, the retention of meltwater in the snowpack, e.g., by capillary forces or by refreezing, is an important process since it delays and reduces runoff and bare ice exposure (Janssens and Huybrechts, 2000). Hence, ablation in the model consists of runoff rather than purely melt. When no snowpack is present, the runoff over the glacier surface equals the amount of melt (Eq. 1), which instantly runs off the surface such that ablation equals melt. When a snowpack is present, however, not all of this melt instantly generates runoff since part of it is retained in the snowpack. In this case the amount of ablation is reduced by the amount of water which is retained in the snowpack (R_w). In the model, it is assumed that all snowmelt remains in the snowpack until it saturates a maximum fraction of the snow depth (P_max) (Table 2). We take the latter to be equal to 0.6 (Reeh, 1991; Janssens and Huybrechts, 2000). Hence, mass is only lost when meltwater is not retained in the firn/snow. At every time step, the snow depth and the amount of retained water are updated.

Other contributing accumulation processes, such as snow drift or avalanches (Turchaninova et al., 2019), or mass removing processes, such as sublimation (Che et al., 2019), may have a substantial influence on the local mass balance, but are not included in the model due to a lack of data. E.g., a study by Turchaninova et al. (2019) showed that for Kara-Batkak glacier, avalanches can have a contribution of 10% to the accumulation (in 2015–2016), but especially along the margins of the terminal part.

The amount of ablation or melt at every time step is calculated using an energy balance model. In this model, the energy fluxes between the glacier and the atmosphere are calculated at every time step. The total energy flux from the atmosphere to the glacier is represented by: Ψ  = ( 1 − α ) ∗ τ ∗ Q in + C 0 + C 1 ∗ T air (2)

Ψ is determined by the flux of the incoming shortwave solar radiation Q_in (W m⁻²), the albedo or reflectivity of the surface (α) and the atmospheric transmissivity (τ). The term [C₀ + C₁ × T_air] represents the other (air temperature dependent) components of the radiation balance, i.e. the sum of the longwave radiation balance and the turbulent sensible heat exchange, linearised around the melting point. T_air is the 2 m air temperature (°C), C₀ and C₁ are tuning parameters which are calibrated for each glacier (see Calibration of the Model and Calibrated Values for C ₀ , C ₁ and Meteorological Gradients). The first term (1 − α) ∗ τ ∗ Q_in or net incoming shortwave radiation is usually the dominant term in Eq. 2. The air temperature dependent fluxes are of intermediate magnitude.

To calculate the incoming solar radiation at every grid cell, it is split up into a direct and a diffuse part, such that a distinction can be made between shaded and non-shaded grid cells. Non-shaded grid cells receive direct solar radiation as well as scattered or diffuse solar radiation, whereas shaded grid cells only receive the diffuse part of the incoming solar radiation. Partitioning the diffuse and direct radiation for each glacier is difficult since no observations of cloud cover are available. Therefore, we implement a typical value of 0.4 for the fraction of diffuse radiation (Table 2) (Voloshina, 2002). Next to that, the amount of incoming diffuse radiation depends on the solar elevation angle. The amount of incoming direct radiation, however, also depends on the angle between the solar rays and the surface. This angle is determined from the slope and aspect of the surface (at 25 m resolution), as well as the solar azimuth. The transmissivity is considered to be constant, at 0.57, which is determined by a comparison between measurements of incoming radiation on Bordu glacier and theoretical incoming radiation at the top of the atmosphere. This value indirectly takes cloud effects into account as the value for clear sky is typically much higher (Voloshina, 2002).

Surface albedo has been proven to be critical for the glacier surface energy balance in the Tien Shan (Che et al., 2019). In this study, the surface albedo is parametrised following Oerlemans and Knap (1998). For every time step, the albedo of snow is calculated using Eq. 3: α snow   =  α firn + ( α frsnow −  α firn ) ∗ e ( − t snow t * ) (3) Where α_firn is the albedo of firn (0.5), α_frsnow is the albedo of fresh snow (0.75), derived from measured values of the weather station on Bordu glacier and a recent study on glacier albedo in the area (Petrakov et al., 2019). t_snow is the time since the last snowfall event (in number of days) and t* is the e-folding time scale determining how fast the snow albedo approaches the albedo of firn after a snowfall event (21.9 days) (Oerlemans and Knap, 1998). Subsequently, the albedo is determined based on the snow depth (d) through: α  =  α snow + ( α ice −  α snow ) ∗ e ( − d d * ) (4) where α_ice is the albedo of bare ice and d* is the characteristic snow depth. Because the glaciers in the Inner Tien Shan are characterised by a high natural dust content (Fujita et al., 2011; Petrakov et al., 2019), as a result of aeolian transport of dust from the deserts and semi-deserts of Central-Asia, the albedo of ice is set at 0.25. This rather low value is confirmed by measurements on Sary-Tor glacier (Petrakov et al., 2019; Rybak et al., 2019) and Grigoriev Ice Cap, which is located nearby (Fuijita et al., 2011), as well as on other glaciers in the area (Voloshina, 2002). d* is set at 3 mm w.e. which was derived for the dry snow season (Brock et al., 2000; Chen et al., 2014). This value is rather low, which appears to be necessary to ensure a sufficiently large albedo during the dry winter months with a thin snowpack on the glaciers in the Inner Tien Shan. Changes in albedo with respect to the temperature at which the snow fell (Giesen and Oerlemans, 2010) and the presence of retained meltwater or accumulation of dirt (Constantin et al., 2020) which might have an important effect locally are not included in the model. The other variables and all constants used in the model are listed in Table 2.

The model is forced at every time step using the temperature and precipitation data from the reconstructed dataset (see Meteorological Datasets for Kumtor–Tien Shan and Meteorological Datasets for Chon-Kyzyl-Suu). For every time step, the air temperature for every grid cell covering the glacier is determined from the air temperature at the nearby meteorological station (T_ref in °C), with an elevation H_ref (m a.s.l.), by interpolating to the elevation of the grid cell (H_sur in m a.s.l.) and using a temperature lapse rate γ_T (°C m⁻¹), as shown by Eq. 5. The lapse rate is determined for every glacier individually by comparing the measured temperatures of nearby meteorological stations for overlapping periods or are it is taken from the literature (see Calibrated Values for C ₀ , C ₁ and Meteorological Gradients). T air = T ref +  γ T ∗ ( H sur  −  H ref ) (5)

The precipitation at each grid cell on the glacier (P_sur in m w.e.) is determined from the hourly precipitation at the nearest reconstruction site (see Reconstruction of Historical Climate Datasets) (P_ref,t in m w.e.), by scaling it up to the elevation of the grid cell, through multiplication with a dimensionless scaling factor P_ratio (Eq. 6). This scaling factor is determined by the vertical precipitation gradient γ_P (m w.e. m⁻¹ y⁻¹), the elevation difference between every grid cell and the station, and the yearly precipitation at the nearby meteorological station (P_ref,y m w.e. y⁻¹) (Eq. 7). P sur , t  = P ref , t ∗  P ratio (6) P ratio = P ref , y +  γ p ∗ ( H sur  −  H ref ) P ref , y (7)

Accumulation in the model is determined by the amount of solid precipitation (snow). The air temperature is used to determine whether precipitation is solid (adding mass to the glacier) or liquid. Typically, the transition from snow to rain occurs at temperatures between 1 and 2°C. To account for this uncertainty, the snow fraction is linearised around T_s (set equal to 1.5°C) with a transition range of 2° (Barandun et al., 2018; Barandun et al., 2021) (Eqs. 8–10). Hence, below 0.5°C, all precipitation falls as snow, and above 2.5°C, it is all liquid. The fraction of solid precipitation (P_s) is therefore calculated using the following expression: P s = P ∗ ( T s − T air + 1 2 )   for   T [ 0 .5-2 .5 ° C ] (8) P s = P   for   T [ <0 .5 ° C ] (9) P s = 0   for   T [ >2 .5 ° C ] (10)

In addition to temperature and precipitation of a nearby location, a Digital Elevation Model (DEM) and a glacier outline are input data for the model. Glacier outlines consist of RGIv6.0 outlines (RGI, 2017), updated using manual delineation on optical satellite data from Sentinel-2 with a spatial resolution of 10 m. The high-resolution DEMs from DigitalGlobe/maxar stereo satellite imagery (HMA DEM) are selected to represent the surface elevation (Shean et al., 2020). The limited number of gaps in these DEMs (occurring in the shaded and steep areas and locally also in the accumulation area) are filled using cubic spline interpolation. To increase the accuracy of the frontal elevation and glacier outline, the HMA DEMs are updated in the frontal area of Bordu and Sary-Tor glaciers using UAV derived Digital Surface Models (DSMs) from 2017 (Bordu) and 2019 (Sary-Tor) (see Table 3) derived from Mavic Pro (Bordu) and DJI Phantom 4 Pro (Sary-Tor) UAV images in the Pix4D software. A slight Gaussian filter is applied to smooth the transition between the HMA and the UAV DSMs.

To reconstruct the historical glacier mean specific mass balance since 1750, it is important to take into account the change in glacier geometry. Since the LIA, the glaciers lost a significant fraction of their glacier tongue. This effect needs to be considered to avoid over- or underestimation of the mean specific mass balance in the past. Therefore, the LIA glacier geometry is reconstructed for the three glaciers.

First, the 2D LIA glacier extent is identified and delineated on satellite images and DEMs, using visible end and lateral moraines attributed to the LIA (Pellitero et al., 2015; Carrivick et al., 2019). For the areas where a lateral moraine can only be observed on one side of the valley, the glacier extent is matched by elevation on the opposite side, following Carrivick et al. (2019). Next, the LIA 3D surface is reconstructed by interpolating the elevations of the end and lateral moraines on both sides of the former ablation tongues (Carrivick et al., 2012, 2019). The interpolated elevation is curved with a smooth parabolic shape with the height of the central glacier axis being approximately 10 m higher than the sides of the glacier. We have no data that support significant geometry changes in the accumulation area. Therefore, in the present-day accumulation area, above the visible lateral moraines, the geometry is kept unchanged, such that the reconstructed height of the LIA surface attaches to the current accumulation area at the equilibrium line (EL). The LIA DEM is then slightly smoothed to have a continuous surface. While for Bordu and Kara-Batkak glaciers, the moraines can easily be identified on recent satellite images, for Sary-Tor glacier, the lateral moraines are inferred from Landsat satellite images from 1999 because the Kumtor concession was extended beyond the former terminal part in recent years. Therefore, the original bedrock and moraine can no longer be recognised (on current satellite images).

The historical reconstruction of the LIA surface shows that at the end of the LIA, the glaciers were substantially larger with a glacier tongue stretching hundreds of metres into their respective valleys (Figure 7). On average, the three glaciers lost about 30% of their area and 40–50% of their volume since the end of the LIA (Table 4). These values are close to those found for glaciers in the Central Tien Shan (Li and Li, 2014). Furthermore, the decline in glacier area and volume has accelerated in recent decades as more recent estimates have shown larger relative declines (Aizen et al., 2007; Chen et al., 2016).

To obtain the glacier’s elevation at every individual year between 1850 (assumed as the end of the LIA, see Little Ice Age in the Kyrgyz Tien Shan) and 2017, the surface elevation is gradually adjusted. To do this, the total elevation difference for each glacier is calculated by subtracting the LIA DEM from the recent DEM. Then, a second order equation is derived for the total elevation difference on each point on the glacier as a function of its current elevation. This equation is then used for the entire LIA glacier outline. Subsequently, for every year between 1850 and 2017, an annual elevation difference is subtracted by dividing these total elevation differences by the number of years (167). Once the bedrock (for the areas outside the current glacier extent) or present-day surface is reached, the elevation is no longer adjusted. A limitation of this procedure is that temporary glacier advances between 1850 and 2020 are not represented (Aizen et al., 2007), nor are accelerations in the decline of glacier size. Furthermore, before the assumed end of the LIA, the glacier size might have fluctuated, but, the LIA geometry is held constant at its maximum extent before 1850 (Figure 8).

Calibration of the model is done for accumulation and ablation separately. Calibration of the accumulation in the model consists of tuning the precipitation gradient, such that the modelled and observed values of winter balance correspond optimally for the different available balance years (Bordu: 2015–2020, Sary-Tor: 2014–2020, Kara-Batkak: 2013–2020). The precipitation gradient is tuned by minimizing the mismatch between the total amount of precipitation during the winter season and the winter balance data of the different glaciers. The winter balance data of the lowest bands are mostly an underestimation of the total amount of precipitation because not all of the precipitation between September and May occurs in solid state at the lower altitudes and because there might occur some ablation as well. Hence, we use all bands except the lowest two for the minimization procedure. The highest bands are also treated with caution as the accumulation and ablation values in these areas consist mainly of approximations as these regions are often inaccessible for measurements (see Mass Balance Measurements).

Then, the constants C₀ (W m⁻²) and C₁ (W m⁻² °C⁻¹) in Eq. 2 are calibrated to establish a bulk relation simulating the effects of the air temperature dependent fluxes, which determine the ablation. An increase in C₀ shifts the net surface energy flux up or down. An increase in C₁ also leads to an increase in the surface energy flux, but in addition determines the influence of temperature on the energy flux. The higher C₁, the steeper this gradient, with highest energy fluxes at lower elevations where temperatures are higher. Hence, the calibration consists of finding the optimal combination of C₀ and C₁, such that the Root Mean Square Error (RMSE) between the observed and modelled mean specific mass balance is minimal for both the individual 100 m bands as well as for the glacier wide mean specific mass balance. Hence, RMSE_tot is minimised following: RMSE b = ∑ 1 n ( MB mod , b − MB meas , b ) 2    n (11) RMSE gl = ∑ 1 y   ( MB mod , gl − MB meas , gl ) 2   n (12) RMSE tot = 1 2 ∗ ( RMSE b + RMSE gl ) (13) with n the number of bands, y the number of years, MB_mod,gl the modelled glacier wide mean specific mass balance, MB_mod,b the modelled mean specific mass balance for each 100 m band, MB_meas,gl the measured glacier wide mean specific mass balance and MB_meas,b the measured mean specific mass balance for each 100 m band.

To have an estimate of the uncertainty of the modelled mass balances, each of the three main determining factors (model parameter values, meteorological input data and observational data uncertainty) is examined. The uncertainty of the modelled mass balances affected by the choice of the model parameters is assessed by running the model 1,000 times for the calibration period with a set of parameters that differ from the original ones (see Table 2). Within the interval of possible parameter values, random values are assigned for each run. Subsequently, the distribution of the modelled mean specific mass balance is analysed, and the standard deviation (σ_par) is used as a metric for the uncertainty. In the altered parameter sets, the values for the albedos are varied with [±0.05], the transmissivity is varied with [±0.05], corresponding to the interannual variability in incoming solar radiation caused by varying cloud conditions (Giesen et al., 2008), the temperature lapse rate is varied with [±1°C km⁻¹] and the precipitation gradient is varied with [±25%]. All the other parameters are varied with ±5%.

Secondly, the uncertainty associated with the meteorological input data is assessed by running the model over the calibration period and perturbating the temperature data with [±0.5°C] and the precipitation data with [±25%] at every time step for 1,000 runs. Likewise, the standard deviation of all the modelled mass balances is used as a metric for the uncertainty (σ_clim).

The third contributing factor is the uncertainty related to the observational data against which the tuning parameters are calibrated. It comprises limitations related to the (individual) in-situ mass balance measurements, as well as to subsequent inter- and extrapolation procedures. The main limitation associated with the in-situ glaciological measurements is the data scarcity in the accumulation area. Regarding the subsequent calculation, the partitioning into elevation bands and the use of a stratigraphic system for the mass balance measurements is the main source of error. The latter is due to the fact that the model runs from 1 October–30 September while the mass balance measurements always refer to two consecutive summer surfaces. The estimated accuracy of the measured mass balances is overall ca ± 0.15 m w.e. yr⁻¹, but this increases to ±0.30 m w.e. yr⁻¹ in the areas where a reduced number of measurements can be performed (accumulation areas). To estimate the uncertainty related to the calibration data, the model is calibrated 100 times using randomly perturbed input mass balance data based on assigned uncertainty intervals. Subsequently, the model is run for the found range of optimal C₀ values, while keeping C₁ fixed, and the standard deviation of the modelled mass balances is used as metric for the calibration uncertainty (σ_cal).

It is clear that the largest uncertainty arises from the chosen value of the model parameters and the least from the perturbed temperature and precipitation data (see Table 5). This means that the evolution of the mass balance over time can be regarded as robust. The uncertainty with regard to the calibration data lies in between the two uncertainties and is rather small. Considering all the above components together using the error law of addition, the combined total uncertainty is estimated at ± 0.28 m w.e. yr⁻¹ for Bordu glacier, ± 0.27 m w.e. yr⁻¹ for Sary-Tor glacier and ±0.31 m w.e. yr⁻¹ for Kara-Batkak glacier for each year (see Table 5). These values are in line with other studies performed on glaciers in the area (Kronenberg et al., 2016; Kenzhebaev et al., 2017).

The historical reconstruction of the cumulative mass balance has an additional uncertainty component related to the changing 3D LIA surface between 1850 and 2017. At the start and end year, we assume that the uncertainty related to the reconstruction is limited and we only consider the uncertainty related to the reconstructed DEM over time. For all years (y) between 1850 and 2017, the uncertainty is set equal to (σ_rec): σ rec = min ( abs ( y  –   2017 ) ,  abs ( y  –   1850 ) ) ( 2017 − 1850 ) ∗ ( MB LIA − MB DEM − 2017 )   (14) such that the uncertainty increases the further away from both end years. Using this method, we find a maximum additional uncertainty of ±0.25 m w.e. yr⁻¹ in 1933–1934.

The calculated average lapse rate for Kara-Batkak glacier based on the local meteorological stations is −5.2°C km⁻¹ (Table 6), which is close to the value found in Satylkanov (2018). Regarding Bordu and Sary-Tor glaciers, the period under observation of the meteorological stations is too short for an accurate estimation of γ_T. Therefore, we use a value of −6.0°C km⁻¹ for these glaciers, based on Aizen et al. (1995). Besides, this value was also used by Rybak et al. (2019) for modeling Sary-Tor glacier and lies between the values used for Batysh-Sook glacier (Kenzhebaev et al., 2017) and Glacier No. 354 (Kronenberg et al., 2016), which are located nearby.

The found precipitation gradients (Table 6) correspond to the rates found in Aizen et al. (1995). In this research, γ_P appeared to be 0.6 mm w.e. m⁻¹ for Bordu glacier, 0.5 mm w.e. m⁻¹ for Sary-Tor glacier and 0.3 mm w.e. m⁻¹ for Kara-Batkak glacier. The smaller value for Kara-Batkak glacier may be due to the more closed nature of the valley in which it is located compared to the more exposed position of the weather station (Chon-Kyzyl-Suu).

Above a critical elevation (H_crit) (Table 6), a decreasing amount of precipitation is implemented to correctly reconstruct the observed MB values with the model. Therefore, precipitation in these areas is calculated using Eq. 15: P sur , t  = P ref , t ∗ P Hcrit , y +  γ p 2 ∗ ( H sur  −  H crit ) P ref , y (15)

Such a negative precipitation gradient in the upper area of a glacier has already been reported for some inner alpine regions (Grünewald and Lehning, 2011) and Central Asian glaciers in particular (Giesen and Oerlemans, 2012). The phenomenon of upper parts of high mountains being drier than the lower slopes can generally be attributed to depletion of precipitable water from the uplifted air mass along the slope (Houze, 2012; Grünewald et al., 2014). For Batysh-Sook glacier, which is located nearby, a linear decrease in accumulation above 4,400 m was also adopted, attributed to wind processes relocating snow and hampering accumulation (Kenzhebaev et al., 2017). Hence, these processes are indirectly included in the calibrated precipitation gradients.

The very large negative gradient for Kara-Batkak Glacier (−4 mm w.e. m⁻¹) (Table 6) is likely also caused by the steepness of the terrain in the highest elevation bands, limiting accumulation. However, the lack of measurements in this area remains a crucial point of uncertainty. Nevertheless, the small size of this highest area of this glacier ensures that the impact of these precipitation amounts is small. Furthermore, concerning Kara-Batkak glacier, between 3,750 and 4,000 m a.s.l, precipitation is not subject to “growth with altitude” due to a combination of steepness and local circulation patterns. As such, a parabolic precipitation reduction factor based on elevation with a minimum of 0.6 is used to match the modelled with the observed accumulation for these elevations.

The calibrated value of 24 W m⁻² °C⁻¹ for C₁ (Table 6) lies within the range of [8.4–40.5 W m⁻² °C⁻¹] found by Giesen and Oerlemans (2012) in a study about the applicability of a similar surface energy mass balance model for glaciers in different climatic regions driven by automatic weather station data. However, the C₀ values found in this study are very negative (Table 6). Nevertheless, (substantial) parameter corrections (e.g., a very low value of τ) were also needed to capture the mass balance of Central Asian glaciers and glaciers in a dry climate in the study of Giesen and Oerlemans (2012). One of the possible mechanisms they proposed to explain this was a too low albedo change to sufficiently reduce the melt rate in dry climates. As such, in this study, the parameters for the temperature-dependent fluxes need to be adjusted to oppose a higher and more realistic value for τ, and correctly reproduce the mass balance. Hence, the very low values for C₀ found in this study are not entirely unexpected. In addition, another way in which the issues of the low albedo change have been resolved in this study is the use of a low value for d* (see Energy Balance Equation).

Furthermore, it is noticeable that the values of C₀ and C₁ for the three glaciers (and especially for Bordu and Sary-Tor glaciers) are rather similar. This suggests that these values are relatively constant in the region. In contrast, the meteorological gradients (especially the precipitation gradient) are more variable and constraining these is more complex. Such climatological variability was also observed in a recent study of Barandun et al. (2021) which revealed a complex climate-glacier interplay in the study area. This likely hampers the model to be applied for glaciers for which there are no or insufficient quality observations to tune the gradients.

We find a close correspondence between modelled and measured mass balances (Table 7). It should be noted, however, that the mass balance model is always applied to the period 1 October–30 September, while mass balance measurements were sometimes carried out earlier. This might explain larger inter-annual differences between the model and the measurements for some individual years (e.g., Sary-Tor in 2015/16, Kara-Batkak in 2017/18) (Table 7). E.g., for Kara-Batkak glacier in 2018/19, the mass balance was already measured in August while in September a significant amount of ablation was found with the model (see Figure 12B). The largest difference between the modelled and measured glacier wide mean specific mass balance is +0.44 m w.e. for Sary-Tor glacier in 2015/16. In this balance year, the overestimation is caused by the less negative MB between 4,200 and 4,400 m a.s.l. (Figure 9).

The graphs (Figure 9) show the measured and modelled mass balance profile for every individual balance year for the three glaciers. Average values for each 100 m band are interpolated to obtain a profile. In general, there is a clear correspondence between the modelled and observed profile. However, for some bands in certain years there are larger deviations (e.g., band 3,900–4,000 m a.s.l. for Kara-Batkak glacier in 2015/16), but this is mainly due to the (inter-annual) variability of the measured profiles. For other years there is an almost complete agreement between measured and modelled mass balance values (e.g., Bordu glacier in 2019/20).

By adding up the local specific mass balance multiplied by grid cell area for each grid cell that covers the glacier and dividing this sum by the surface area of the glacier, the glacier wide mean specific mass balance is determined (Figure 10). Also, for the mean specific mass balance, there is a close agreement between the measured and modelled values. The small deviations that are present in certain years are eliminated when the cumulative sum is determined over the period of observation, as shown by the modelled cumulative mass balance, which is very close to the observations for all three glaciers. Only for Kara-Batkak glacier, there is a relatively small underestimation at the end of the period. However, this is mainly due to the 2017/18 and 2018/19 balance years in which the modelled mean specific mass balance is lower than the measured value. It is striking that Kara-Batkak was the only glacier in Kyrgyzstan with a less negative mean specific mass balance in 2018/19 compared to 2017/18 (WGMS, 2020), which is likely caused by the early measurement of the summer balance in August 2019.

The maps with the average mass balance clearly show that the three glaciers have a small accumulation area (MB > 0) and a large ablation area (MB < 0) (Figure 11). For the years under consideration, the average accumulation area ratios (AAR) were 0.36 for Bordu glacier, 0.41 for Sary-Tor glacier and about 0.51 for Kara-Batkak glacier. For the latter, this value was about 0.7 in the 1970s (Kutuzov and Shahgedanova, 2009), which indicates that the EL rose significantly. The small AAR values indicate that the three glaciers are currently all out of balance with the present meteorological conditions.

To reveal the path of the cumulative mean mass balance throughout a given balance year, we investigate the evolution from day to day (Figure 12). Little happens between October and March. The cumulative mean mass balance increases slightly due to some snow events but remains relatively close to zero. This is attributed to the blockage of the Siberian high, which provides cold but predominantly dry weather during winter. Between April and early July, the cumulative mean mass balance increases substantially and rapidly. This is due to precipitation episodes that largely fall as snow on the glacier surface. The ablation season especially occurs throughout July and August, resulting in a strong decrease of the cumulative mean mass balance. In September, the amount of ablation decreases with in general by the middle to end of September a minimum in the mass balance. Hence, the glaciers can be considered as examples of the spring-accumulation type with a prolonging accumulation to the early summer period (June). It is also clear that measuring the mass balance in August usually leads to an underestimation of the amount of melt during the ablation season because in general, the cumulative mean mass balance continues to decrease into September (Figure 12).

The comparison with the WGMS data from Kara-Batkak glacier between 1958 and 1998 (see Mass Balance Measurements) shows an acceptable agreement for the results with the adjusting DEM since the LIA (Mean Absolute Error (MAE) = 0.20 m w.e. yr⁻¹, Root Mean Square Error (RMSE) = 0.32 m w.e. yr⁻¹) (Figure 13). The largest differences are found for 1984/85 and 1985/86 for which the mean specific mass balance in the WGMS dataset appears to be very low (< −1.50 m w.e. yr⁻¹), and for 1996/97 for which the mean specific mass balance in the WGMS dataset is very high compared to the data of Tuyuksu glacier, the nearest glacier with MB measurements (Hoelzle et al., 2017). It must be stressed that the mean specific mass balance data in the WGMS dataset were reconstructed based on the water level of the proglacial lake of Kara-Batkak glacier and assumptions on the precipitation gradient, which are very likely accompanied by a substantial error margin. In addition, for the years with the largest deviations, the WGMS data is not consistent with the temperature and precipitation data of those years and their relationship to the other years. Nevertheless, the agreement between model results and the WGMS data is acceptable by the end of the measurement period (1998). There is a difference in cumulative mean specific mass balance of only −4.93 m w.e. which corresponds to −0.12 m w.e. yr⁻¹. Further, the decrease towards a more negative mass balance from 1976 onwards is clearly represented by the model results using the adjusting geometry. It is also clear that the use of the current glacier geometry or LIA geometry would over- or underestimate the cumulative mass balance substantially (Figure 13).

The curve obtained with the adjusting geometry lies midway between the version with the LIA geometry and the current geometry, from around 1880–1900 onwards (Figure 14). For Kara-Batkak, the separation of the LIA DEM version and the adjusting DEM version starts later (Figure 14B), which is caused by the slower retreat between 1850–1900 due to a large ice thickness of the LIA glacier front. For the curves with the adjusting geometries and the fixed LIA DEM, a slight decrease of the cumulative mass balance can already be recognised between 1750 and 1850. This might be due to an overestimation of the glacier geometry in this period. A closer inspection revealed that especially at the front of the glaciers, the modelled MB in this period was negative, leading to a mean specific MB just below 0 m w.e yr⁻¹. Probably, the glacier geometry around 1750 was somewhat smaller than the LIA maximum, such that the real curve between 1750 and 1850 is most likely located between the 2017 version and the LIA version with hence an almost zero cumulative mean specific mass balance. As a result, the cumulative value obtained in 2020 for the adapting geometry may be slightly too negative.

Nevertheless, the results for the vast majority of years reveal a negative mean specific mass balance whereby the cumulative value is (strongly) decreasing since 1850 (Figure 14). There are periods with a near-zero or even positive mean specific mass balance, especially at the beginning of the 20^th century. The latter corresponds to the last period of glacier growth, as was demonstrated by Solomina et al. (2014), resulting from a series of colder years with increased precipitation (Figures 5, 6). A trend towards an increasingly negative mass balance can be observed in the last decades, starting approximately in 1970.

For the runs with the 2017 DEM version, the cumulative mass balance remains roughly stable (and even rises) until 1900–1910. Afterwards, the cumulative mass balance decreases at an increasing rate, but overall, it remains at a high value. The curve based on the LIA geometry on the other hand, shows a much stronger decrease in the cumulative mass balance, especially over the last decades. This is due to the remaining large LIA glacier tongue resulting in very negative mass balance values over the lowest parts of the ablation area of the LIA geometry. In reality, however, glaciers constantly respond to such negative mass balance trends and adapt to changing climatic conditions through a change in their shape (retreating front), see Reconstruction of Little Ice Age Glacier Geometry.

The sensitivity of the model with respect to the different parameters is calculated by changing the parameters and analyzing how it affects the modelled mass balances, following Che et al. (2019). The sensitivity is defined as %MB change per 10% parameter change. This range is selected to determine how much variation in the model results occurs for a fixed amount of uncertainty in the parameters, as in Anslow et al. (2008). The modelled mass balances are clearly most sensitive to the albedo of snow, the transmissivity, and the temperature gradient for all three glaciers (Table 8). The latter determines both the amount of melt (by the parameterised long wave radiation and sensible heat exchange), and the distinction between snow and rain. The significantly larger sensitivity of the results to the albedo of snow compared to the albedo of ice (and firn) (Table 8) was also observed by Che et al. (2019) and can likely be attributed to the importance of summer snow which regularly interrupts the melting season. This is why glaciers with spring-accumulation and summer snow showers are more sensitive to temperature increases during the summer months. For the implementation of the model, it is therefore crucial to carefully determine these parameters in particular.

An independent validation of the mass balance reconstruction is performed by comparing it with geodetic mass balances, derived from DEM subtraction (Table 9). In general, we find a close agreement. However, for certain periods, the differences are more significant. For example, the geodetic mass balance calculated for 2000–2018 for Kara-Batkak glacier derived by Shean et al. (2020) indicates only a very limited mass loss while the observations and the modelled values are much more negative. This could indicate a substantial underestimation of height differences for Kara-Batkak glacier in their study. The best correlations are found for Sary-Tor glacier, for which all compared values correspond closely. This thus constitutes a thorough quality control of the modelled mass balance series, as it indicates that there is no significant over- or underestimation for the periods covered by all the geodetic mass balances.

As mentioned in the introduction, the mass balances of glaciers in the study area were reconstructed in several other studies. Amongst others Popovnin et al. (2021) used a simple relationship to reconstruct the mean specific mass balance of Sary-Tor glacier since 1930. In their study, the glacier wide total accumulation was reconstructed using a linear relationship with air temperature and precipitation sums measured at the Kumtor-Tien Shan meteorological station. Ablation in their approach was based only on average temperatures between May and August, as derived by Konovalov (1977). They found an average value of −0.53 m w.e. yr⁻¹ for the period 1985–2019. This ties well with the value found in this study (−0.52 ± 0.05 m w.e. yr⁻¹). It shows that their simple parameterisation was able to produce results similar to the energy balance model used in this study. For the period 1930–1985, Popovnin et al. (2021) found an average value of -0.31 m w.e. yr⁻¹, while in this study we found a value of −0.39 ± 0.05 m w.e. yr⁻¹. However, when comparing these results, it must be pointed out that they used a constant present-day glacier geometry, whereas in this study the glacier geometry was adapted since the LIA. In the 1950s, Sary-Tor glacier had a much larger glacier tongue (Aizen, 2007), so the mean specific mass balance was much more negative. Further, the modelled mass balance in this study is more in line with the observed geodetic mass balance (see Comparison With Geodetic Mass Balances), which indicates that an adjustment of the glacier geometry is crucial for historical reconstructions.

Kronenberg et al. (2016) reconstructed the seasonal mass balance of Glacier No. 354, located in the valley to the west of Bordu glacier, using a distributed accumulation and temperature-index melt model and data from the Kumtor-Tien Shan station. They found a close agreement between modelled (−0.40 m w.e. yr⁻¹) and observed (geodetic, −0.48 m w.e. yr⁻¹) mean mass balance for the period 2003–2012. For the same period, we find a mass balance of −0.41 ± 0.10 m w.e. yr⁻¹ for Sary-Tor glacier and −0.40 ± 0.10 m w.e. yr⁻¹ for Bordu glacier. These values are very similar and seem to indicate that the mass balance of these three adjacent glaciers located in adjacent parallel valleys is very consistent.

A comparison between the modelled mass balance values for Kara-Batkak glacier with the reconstructed annual mass balance values (2000–2018) of Barandun et al. (2021), obtained from transient snowline and geodetic surveys, reveals a rather large mean absolute difference of 0.50 m w.e. yr⁻¹. However, over the 18 years, the cumulative difference is only 0.69 m w.e. A comparison with the mass balance observations for the overlapping period shows that the reconstructed mass balance from Barandun et al. (2021) deviates strongly in some years, while our model results are more in line with the observations. For example, for the 2015–2016 balance year, MB measurements revealed a mean specific mass balance of −0.39 m w.e. (see Figure 3), the model results showed a mean specific mass balance of −0.37 ± 0.30 m w.e. (see Figure 10B), while the results of Barandun et al. (2021) revealed −1.69 m w.e.

For all three studied glaciers, we find an average mass balance gradient between -0.50 and −0.60 m w.e. yr⁻¹ per 100 m. This value closely matches the value of −0.45 m w.e. yr⁻¹ per 100 m for the nearby Batysh Sook glacier for the period 2004–2016 (Kenzhebaev et al., 2017) and is slightly higher than the values found for other nearby glaciers in the past (Aizen and Zakharov, 1989; Ushnurtsev, 1991). However, Kronenberg et al. (2016) found a value of −0.68 m w.e. yr⁻¹ for Glacier No. 354. Research showed that an ongoing change towards higher gradients might be expected for summer-accumulation type glaciers (Oerlemans, 2001), since for this type of glaciers, a small change in summer accumulation can significantly impact the surface type in the ablation area (snow/ice) and thus considerably change the mass balance gradient (Kronenberg et al., 2016). The same can be said for spring-accumulation type glaciers.

The lack of validation data for the historical reconstruction up to 1750 does not allow for an independent quality analysis, such as for the last decades, by comparing with geodetic mass balances. It is therefore crucial to point out that the modelled MB series depends on the climatological series that has been reconstructed. Intrinsic problems related to the temperature and precipitation reconstruction from tree rings are the loss of variability and the fact that every month exhibits the same positive or negative signal (see Temperature and Precipitation Forcing). Besides, a drought index has been used as a measure for summer precipitation. As a result, the precipitation series only encompasses the signal surpassing the bi-annual variability.

For the observational part of the climatological series, reliability and homogeneity were maximised by carefully searching for correlations between the different datasets, by correcting for possible biases and by comparing reconstructed values with observations. Here, the main limitations are possible biases or errors in the measured data from the weather stations (e.g., due to the relocation of the Kumtor-Tien Shan station) or by using the monthly anomalies created from different weather stations in the area (such as the CruTEM or CruTS datasets).

Other limitations are related to the simplicity of the model. Different processes such as refreezing and sublimation are not included as independent processes. However, they are included indirectly as a part of accumulation (calibrated through the precipitation gradient) and ablation (calibrated in the value of C₀ and C₁) in the model. Although assumed to be of minor importance, the melt rates (ablation/runoff) might be overestimated because of the absence of sublimation. Another limitation of the study is the fixed geometry of the accumulation area. There is no direct historical evidence that the accumulation areas reached higher elevations or were much broader during the LIA, but it is possible that accumulation was slightly underestimated in the historical reconstruction. To better constrain the geometry of the glaciers for the historical reconstruction, a three-dimensional glacier model with included ice dynamics could be used.