Assessment of changes in mass balance of the Tuyuksu group of glaciers, northern Tien Shan between 1958 and 2016 using ground-based observations and Pleiades satellite imagery

Continuous measurements of glaciological mass balance have been conducted at the Central Tuyuksu glacier, Tuyuksu group of glaciers, Ile Alatau, northern Tien Shan since 1957, showing that cumulative mass balance was negative since the 1970s. Geodetic mass balance was calculated for the 1958–1998 and 1998–2016 periods using multi-temporal digital elevation models derived from the historic photogrammetric surveys from 1958 and 1998 and the high-resolution Pléiades satellite stereo imagery


INTRODUCTION
There is ample evidence for glacier wastage in Tien Shan, Central Asia (Unger-Shayesteh et al., 2013). In the last two decades, research focused mainly on documenting changes in the glacier area (Kutuzov and Shahgedanova, 2009;Narama et al., 2010;Severskiy et al., 2016) as repeated satellite imagery has become widely available. However, from the perspective of changes in global sea level (Zemp et al., 2019) and in regional water resources (Huss and Hock, 2018;Shahgedanova et al., 2018Shahgedanova et al., , 2020, the changes in glacier volume and in mass are most important. These changes have been assessed using geodetic mass balance method applied on a regional scale and measurements of the mass balance of individual glaciers using the glaciological (stake) method. Estimations of geodetic mass balance in Central Asia focused on comparisons of the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) obtained in 2000 with earlier DEM built from the declassified Corona and Hexagon imagery from the 1970s (e.g., Pieczonka and Bolch, 2015). While an advantage of this method is the wide regional coverage and a relatively long time span, the uncertainties of mass balance measurements are high. In many regions of Tien Shan, they are comparable with the estimated values of mass balance. Distributed geodetic mass balance was derived for the wider region of High Asia by Brun et al. (2017) using SRTM and DEM derived from ASTER stereo imagery for a single time step of 2000-2016. Geodetic mass balance was calculated for several glaciers in Kyrgyzstan using differencing of multiple DEMs produced from high-resolution stereo imagery (Barandun et al., 2018).
Measurements of glacier mass balance using glaciological (stake) method have been conducted on several benchmark glaciers in Tien Shan since the first International Geophysical Year of 1957-1958(Zemp et al., 2009). However, on most glaciers, the measurements were interrupted in the 1990s, resulting in more than 20 years of gap in data. On several glaciers in Uzbekistan and Kyrgyzstan, mass balance monitoring programs were renewed in the 2010s (Hoelzle et al., 2017). The gaps in the long-term mass balance records of several glaciers were filled using a combination of mass balance modeling and analysis of satellite data (Barandun et al., 2018).
Uninterrupted measurements continued on two benchmark glaciers, the Central Tuyuksu (often referred to as Tuyuksu) in the Ile Alatau, Kazakhstan, and the Urumqi Glacier No. 1 in China. The Tuyuksu mass balance monitoring program, which started in 1957, involves regular measurements of ablation during the melt season and surveys of snow depth at the end of the accumulation season using a network of over 100 stakes. Annual values of accumulation, ablation, and the resulting mass balance are reported to the World Glacier Monitoring Service (WGMS). Geodetic surveys are a part of the monitoring program and the glacier area has been mapped annually since 2006. There are two high-resolution DEM of the Tuyuksu and the surrounding glaciers from 1958 and 1998 (Severskiy et al., 2006). These DEMs were used to assess changes in glacier surface elevation, revealing strong glacier thinning that reached 50 m on the glacier tongue (Severskiy et al., 2006), and to verify the continuing glaciological measurements (Hagg et al., 2004). The latest verification using the 1958 and the 1998 DEMs and conducted for the Central Tuyuksu glacier showed a discrepancy between the geodetic mass balance of -12.6 m w.e. and the cumulative glaciological mass balance of -16.8 m w.e. and attributed it to the uncertainties of the glaciological method amplified in the cumulative record (Hagg et al., 2004).
More recently, however, there were no assessments of geodetic mass balance of the Tuyuksu glacier and the surrounding glaciers using contemporary high-resolution DEM. Such an assessment is needed in order to (i) assess temporal changes in the rates of glacier wastage and (ii) continue the validation of mass balance measured in situ under the conditions of rapid climatic warming. The stereo imagery from the Pléiades satellites was recently used to derive high-resolution DEM of glaciers and assess changes in elevation of glacier surface with an accuracy of ±1 m or better at multiple sites in the Andes, European Alps, Himalayas, and Iceland by Berthier et al. (2014), in Central Tien Shan by Barandun et al. (2018), and in Caucasus by Kutuzov et al. (2019). Pléiades is a constellation of satellites launched by the French Space Agency (CNES) in December 2011. It is envisaged that its application, coordinated by the Pléiades Glacier Observatory project 1 , will complement in situ observations of glacier mass balance conducted at the WGMS reference glaciers. This study is the first attempt of deriving and using Pléiades DEM of the Tuyuksu group of glaciers. This paper presents three high-resolution DEMs derived for the Tuyuksu group of glaciers from the photogrammetric surveys of 1958 and 1998 and from the Pléiades stereo imagery for 2016. Changes in surface elevation and geodetic mass balance are calculated for three time periods from these DEMs. Geodetic mass balance is compared with cumulative mass balance derived from stake measurements. The observed glacier change is discussed in the context of continuing climatic warming and its impact on the glaciers of the region.

STUDY AREA
The Tuyuksu group of glaciers is located in the headwaters of the Kishi Almaty River on the northern slope of the Ile (Zailiiskiy) Alatau (Figure 1) approximately at 43 • N and 77 • 06' E and between 3,500 and 4,300 m a.s.l. The region is characterized by strong seasonal variations in temperature and precipitation. The westerly flow dominates in autumn and in spring, resulting in the precipitation maxima in April-May on the plains shifting toward June-July in the middle and high mountains, where snow accumulation peaks in spring and in early summer. The snow and glacier melt period is limited to June-July-August, extending to September in individual years.
At the beginning of the monitoring program, the Tuyuksu group included nine glaciers with a combined area of 7.6 km 2 in 1958. Central Tuyuksu is the largest glacier of the mountainvalley type, with an area of 2.99 km 2 in 1958 and declining to 2.28 km 2 in 2016. Between 1958 and 2016, the accumulation-area ratio (AAR) averaged 41%. The equilibrium line altitude (ELA) was positioned on average at 3,830 m, with a standard deviation of ±110 m. Two other glaciers had areas in excess of 1 km 2 .
Between 1990 and 2008, the glaciers of the Ile Alatau were losing their area at an average rate of 0.9% a −1 , although the rates vary between the catchments depending on their elevation and the size of the glaciers (Narama et al., 2010;Severskiy et al., 2016). These changes have been attributed largely to increasing air temperatures and also to a negative anomaly in precipitation observed in the 1970-1980s (Severskiy et al., 2006(Severskiy et al., , 2016Shahgedanova et al., 2018), which was forced by changes in the atmospheric circulation and affected most of Tien Shan (Cao, 1998).
The changes in area of the glaciers of the Tuyuksu group correlated well with the changes of the total glacier area in the Ile and Kungey Alatau system (Severskiy et al., 2016). The cumulative mass balance of the Central Tuyuksu remained negative and the annual mass balance was mostly negative, with the exception of 15 years, mostly in the 1960s and in the 21st century (Severskiy et al., 2016).  Tuyuksu glacier are shown for 1958, 1998, and 2016. Pléiades image is used as background.
Frontiers in Earth Science | www.frontiersin.org  (Simon et al., 1961;Eder et al., 2002Eder et al., , 2005. The maps covered nine glaciers, including the Central Tuyuksu itself. Both maps were scanned with a pixel resolution of 1,200 dpi and automatically digitized using ArcScan module of Arc GIS. The centerline method of automatic vectorization was used. Following vectorization, both maps were re-projected to WGS 84, zone 43 N projection. Two DEMs were developed from the digitized and the georeferenced maps using ArcGIS 10.5 software. The elevations were interpolated on regular grids with a resolution of 10 m using the kriging method. A network of 36 GCP (including five GCP established in 1958) obtained in August 2018 and September 2019 using Topcon and Leica DGPS, respectively, was used for co-registration and correction of DEMs. All GCP were located on the firm terrain, which did not experience a change in elevation due to ice melt or ground movement. Planimetric correction was not required. The elevation differences between the DEMs and the GCP are shown in Table 1 and these are assumed to be due to errors in DEMs. On average, GCP elevation was 2.9 m higher than the elevation of the corresponding pixels in the 1958 DEM, while the 1998 DEM showed a very close match.
Vertical displacement is a function of aspect and terrain slope (Nuth and Kääb, 2011) and it is desirable to assess and, if required, remove spatially varying elevation errors. In this study, the DEMs covered a relatively small area located on the northern slope of Ile Alatau. Most glaciers have northerly and westerly aspects. A small number of GCP from higher elevations did not allow us to produce a comprehensive assessment of errors due to slope steepness relative to GCP. Therefore, vertical bias was removed from both DEMs by subtracting the mean values of vertical shift registered off ice from all DEM pixels (Table 1).

Pléiades DEM From 2016
The Pléiades 1A and 1B twin satellites, launched in 2011 and 2012, respectively, generate stereo images with a pixel size of 0.5 m for the panchromatic channel and 2.5 m for the multispectral blue, green, red, and near-infrared channels 2 . For DEM generation, the multispectral data are preferred due to better point matching afforded by multi-channel images (Gim and Shin, 2015). A pair of multispectral images, obtained at the end of the ablation season on 27 August 2016, was used in this study (Figures 1, 2A).
The Pléiades DEM was generated using Leica Photogrammetry Suite software using the Automatic Terrain Extraction processing engine, with the rational polynomial coefficients distributed as metadata with the imagery. These data describe the geometric analytic model between image space and ground space and are based on the direct scene orientation derived from a combination of GNSS, star sensors, and giros (Jacobsen and Topan, 2015). We refined this model using 36 GCP distributed across the image area. A 2.5-m-resolution DEM was derived from the point cloud data derived from automatic image matching. The obtained DEM was co-registered to 36 GCPs, and statistics characterizing the uncertainties of the DEM on stable terrain were calculated following the co-registration ( Table 1).
The Pléiades DEM was resampled from 2.5 m to a coarser resolution of 10 m using bilinear interpolation method in order to match the spatial resolution of the earlier DEM.

Co-registration of DEM and Accuracy Assessment
Changes in surface elevation were calculated for three time periods as defined by the DEM acquisition dates. Three pairs of co-registered DEM were produced using the re-sampled Pléiades DEMs (1958-1998, 1998-2016, and 1958-2016) to calculate changes in elevation of the terrain by subtracting a later DEM from an earlier DEM for each pair.
At this step, DEMs that was adjusted for bias following the co-registration with GCP were used. There were no planimetric shifts between the DEMs relative to each other. The accuracy of vertical fit between DEMs in each pair was evaluated by computing the pixel-by-pixel difference in elevation (bias) for the ice-free areas where no elevation change was expected (Berthier et al., 2007;Agarwal et al., 2017). As errors in DEM increase with increasing slope, the accuracy assessments are often restricted to the areas with slope less than 30 • (e.g., Pieczonka et al., 2011). An analysis of error distribution with elevation in the study area confirmed this assumption and 30 • slope value was used as a cutoff point for selecting pixels for accuracy assessment.
The error statistics for stable, ice-free terrain are shown in Table 2. Despite the initial correction for elevation bias using GCP, large elevation differences were still present in each pair of DEMs. Following Gardelle et al. (2013), we considered the elevation differences outside the range of mean ± three standard deviations as outliers. The outliers were removed and the relative vertical error statistics were recalculated, showing some improvement in standard deviations ( Table 2).
The difference between the DEMs within stable terrain following the removal of outliers was used to quantify uncertainty in elevation change. We used median absolute deviation (MAD), calculated using Eq. (1), to quantify uncertainty in elevation change. The MAD values, shown in Table 2, were applied to the entire footprints of each pair of DEMs.
where X n is a sample of differences in elevation between two DEMs used in the accuracy assessment.
To identify outliers within the glacierized area, statistics of elevation change were calculated for 50-m altitudinal belts (Gardelle et al., 2013). Pixels, in which elevation change exceeded mean ± three standard deviations as calculated for a given altitudinal belt, were removed. While few outliers were identified in the ablation zone (particularly at the Central Tuyuksu which has a gentle slope), unrealistic elevation differences were present in areas of steep terrain and deep shadow in the 1998-2016 and 1958-2016 pairs of DEMs built from the optical Pléiades imagery.

Assessment of Elevation Change and Calculation of Geodetic Mass Balance
Elevation changes were calculated for the glacierized areas (open parts of glaciers) including the Central Tuyuksu (Table 3) and eight other glaciers of the Tuyuksu group ( Table 4). The supraglacier debris cover is almost absent in the region. We note that the 2015-2016 mass balance year was characterized by high accumulation, resulting primarily from precipitation in April-June, and snow was present on the glaciers at the time of the Pléiades imagery acquisition on 27 August 2016. This could potentially introduce additional uncertainty to the calculation of changes in glacier surface elevation and geodetic mass balance. In order to establish whether snow cover melted between the acquisition date and the end of the mass balance year (and therefore the DEM required adjustment) or persisted into the   indicating the end of the 2015-2016 mass balance year. Figure 2 shows that there were no significant changes in the position of the snow line between the two images, while the snow depth measured in situ on 31 August and 24 September 2016 did not change significantly either. Therefore, no adjustments were made to the Pléiades DEM. The area-averaged specific geodetic mass balance was calculated using Eq. (2). The mean density of 850 ± 60 kg m −3 was used for volume-to-mass conversion for the calculation of geodetic mass balance according to the recommendations by Huss (2013) and the calculation of the geodetic mass balance of other glaciers in Tien Shan (Barandun et al., 2018).
where V is change in volume (m 3 ),Ā is the mean of glacier area at the start and at the end of a time interval (m 2 ), and ρ is density (kg m −3 ). To calculate the mass balance rate (m w.e.a −1 ), this value was divided by the number of years between the two DEMs in each pair. The uncertainty in the calculation of glacier volume change included two terms: relative vertical error and uncertainty of glacier area measurement. The former was quantified using MAD ( Table 2). The uncertainty in glacier area measurements was quantified using the buffer method, whereby an area of a buffer with a width of 5 m (map accuracy and half of the DEM pixel size), drawn around the glacier shapefile, was normalized by the glacier area.
The uncertainty in geodetic mass balance calculations included the vertical and the area uncertainty terms and the uncertainty due to density assumptions (Huss, 2013). The overall uncertainties in volume and mass balance were calculated according to the standard rule of error propagation as the square root of the sum of the squares of the individual uncertainty terms.

Glaciological Mass Balance
Glaciological mass balance has been measured at the Central Tuyuksu since 1957. Currently, a network of over 100 stakes is used. The maximum snow accumulation is measured at the end of the cold season (winter mass balance) and the mass decrease is measured from maximum snow accumulation until the end of the ablation season (summer mass balance). The sum of winter and summer mass balance is reported as annual mass balance. Snow density is measured on a regular basis at several snow pits located at different elevations.
Although the mass balance program has been run without interruptions, a reduced number of stakes and less frequent surveys took place in 1991-1993 and in 2002-2005. During these periods, the number of stakes declined from 140-150 used before 1991 to approximately 60 and increasing to 100-120 between these periods and after 2005. Importantly, despite a dense network of stakes, there is paucity of measurements in the accumulation zone of the glacier due to accessibility problems.
The uncertainty in glaciological mass balance measurements originates from (i) height determination, (ii) density measurement errors, (iii) presence of superimposed ice, (iv) changing location and density of stakes, (v) spatial averaging and interpolation, (vi) under-sampling in the inaccessible areas, and (vii) using a constant glacier area (Zemp et al., 2013). A large number of stakes, regular density measurements at different All changes and mass balance values are negative. The total geodetic mass balance and mean annual rate of geodetic mass balance of all glaciers are calculated from the total changes in area and volume.
Frontiers in Earth Science | www.frontiersin.org elevations, and annual area measurements minimize all sources of uncertainty except (iv) and (vi). Vilesov and Uvarov (2001) estimated that the absolute error of the glaciological method varied in different years between ±0.04 m w.e. and ±0.06 m w.e. in the region of direct measurements. This assessment included error sources (i)-(iii) but did not take into account the reduced number of stakes. To test the sensitivity of glaciological mass balance to a changing number of stakes, we re-calculated the annual mass balance for the whole glacier in the period 2007-2016 using those stakes which were retained in 1991-1993 and in 2002-2005. The root mean square error for annual mass balance was ±0.02 m w.e. This value represents 5% of the mean annual mass balance averaged over the whole record of glaciological measurements (−0.40 m w.e.a −1 ). The impact of under-sampling in the accumulation zone is more difficult to evaluate. Currently, the mass balance above 3,800 m a.s.l. is calculated using regression equations based on 20 years of direct measurements which link mass balance in the 100-m elevation bands to mass balance measured in the 3,700-3,800-m band (Makarevich, 2007).

Changes in Surface Elevation, Volume, and Mass Balance of the Central Tuyuksu
The changes in elevation of the Central Tuyuksu glacier surface are presented in Figures 3, 4 and the mean statistics are summarized in Table 3. Between 1958 and 2016, the glacier surface lowered on average by 23.2 ± 2.2 m. The altitudinal distribution of elevation changes confirms the expected pattern with glacier downwasting in the ablation zone and thickening in the accumulation zone. In the glacier terminus area, downwasting exceeded 60 m. On average, the annual rates of area loss increased slightly in 1998-2016 in comparison with those in 1958-1998 ( Table 3). The mean rate of glacier surface lowering decreased, but the change in the annual rates was close to the uncertainty of measurements. The altitudinal distribution of changes in the elevation of glacier surface shows that high rates of glacier downwasting (reaching 2-2.5 m a −1 ) affected a much larger area of the glacier in the period 1998-2016, extending to 3,600-3,700 m a.s.l. (Figures 3, 4). At the same time, glacier thickening increased in the accumulation zone in the period 1998-2016 (Figures 3-5). Overall, (67.7 ± 6.7) × 10 6 m 3 of ice was lost from the Central Tuyuksu between 1958 and 2016, giving a negative geodetic mass balance of −21.8 ± 2.6 m w.e. The differences between the annual mass balance rates in 1958-1998 and in 1998-2016 were within the accuracy of the measurements ( Table 3).

Changes in Surface Elevation, Volume, and Geodetic Mass Balance of the Glaciers of the Tuyuksu Group
The trends in loss of area and surface lowering of the other glaciers of the Tuyuksu group were similar to those of the Central Tuyuksu (Figures 6, 7 and Table 4). As these glaciers are smaller than the Central Tuyuksu, the mean area loss was higher, averaging 48% between 1958 and 2016 (Glacier No. 2 melted completely) against 23.4% lost by the Central Tuyuksu. The mean lowering of glacier surface was 17.8 ± 2.2 m against 23.2 ± 2.2 m (Tables 3, 4). Similarly to the Central Tuyuksu, the annual rates of surface lowering increased at lower elevations in the second assessment period (Figure 7). Glacier thickening was observed in the accumulation zone in the second assessment period, while it was close to zero in the first assessment period.
The uncertainty of measurements was higher at three glaciers-Molodezhniy, Kosmodemjanskaya, and Ordzhonikidze-in comparison with the Central Tuyuksu. While at the Central Tuyuksu the changes observed during the period 1958-2016 were very close to the sum of changes measured in 1958-1998 and in 1998-2016, at these three glaciers, there were discrepancies between the changes in surface elevation and volume as measured using the 1958 and 2016 pair of DEMs and a sum of changes in the two assessment periods (Table 4). At the other glaciers (Igly Tuyuksu, Mayakovskiy, Mametova, and Glacier No. 1), the changes measured in 1958-2016 were in good agreement with the combined changes in 1958-1998 and in 1998-2016. The quality of the DEMs is the most likely explanation for the observed discrepancies as most GCP were located in proximity to the Central Tuyuksu, which was the primary focus of this study. Both Molodezhniy and Ordzhonikidze are hanging glaciers positioned on steep slopes, where errors in DEM are higher (Pieczonka et al., 2011).

Comparison Between Geodetic and Glaciological Mass Balance of the Central Tuyuksu
The differences between the geodetic mass balance for the periods 1958-1998 and 1958-2016 and the cumulative glaciological mass balance are compared in Table 3 and Figure 8. The two methods were in good agreement in both assessment periods. Overall, the geodetic mass balance was less negative than the glaciological mass balance, but the observed differences were within the uncertainty of measurements. Hagg et al. (2004) and Zemp et al. (2013) provided a comparison of geodetic and glaciological measurements at the glaciers with the well-established mass balance monitoring program, showing that the differences in annual mass balance values were within the 0.01-0.4 m w.e. a −1 range. In the case of the Central Tuyuksu, the annual values of geodetic mass balance were −0.39 ± 0.05 and −0.35 ± 0.18 m w.e.a −1 in the two assessment periods versus the glaciological mass balance of −0.40 m w.e.a −1 in both periods, which is within the range of differences obtained in other studies. Hagg et al. (2004) used the same geodetic surveys of 1958 and 1998 to develop DEMs with a resolution of 20 m and to derive the geodetic mass balance of the Central Tuyuksu. Their value of −12.6 m w.e. showed a 33% difference with the cumulative mass balance, which was attributed to the uncertainties of the glaciological method. In our study, the geodetic mass balance of the Central Tuyuksu was calculated as −15.7 ± 2.1 m w.e. Our comparison showed 5 and 14% difference between the two methods in the 1958-1998  and the 1998-2016 assessment periods, respectively, and 8% difference overall. We attribute this improvement to the coregistration of DEMs with a network of GCP established by DGPS surveys and the subsequent correction of DEMs ( Table 2). These facilities were not available previously, and Hagg et al. (2004) assumed no off-ice difference between the 1958 and the 1998 DEMs, in contrast to this study ( Table 2). We note a slightly higher discrepancy between the glaciological and the geodetic mass balances in the period 1998-2016, which may be due to the lack of direct measurements in the accumulation zone and thus not capturing an increase in accumulation (Figures 3, 4).

DISCUSSION
The record of the Central Tuyuksu glaciological mass balance is the longest continuous record in Central Asia, providing valuable data for the assessment of the impacts of climate change on glaciers and water resources in the region, which depends on glacier melt for water provision more than any other part of the world (Kaser et al., 2010;Pritchard, 2019). However, as with any long-term observations, the record of the glaciological mass balance from the Central Tuyuksu requires verification, and the best way to achieve it is through the estimation of the geodetic mass balance at regular intervals (Cogley, 2009;Hoelzle et al., 2017). The estimation of geodetic mass balance relies on the availability of high-resolution DEM. Currently, various DEMs are available, e.g., the widely used SPOT and ALOS, a more recent High Mountain Asia 8-m DEM derived from cross-track and along-track optical imagery (Shean, 2017), and Pléiades stereo imagery (Berthier et al., 2014), available through the WGMS Pléiades Glacier Observatory. For the Central Tuyuksu, two geodetic surveys, conducted in 1958 and in 1998, provided two historical DEMs which, together with the DEM derived from the 2016 Pléiades stereo imagery, enabled the calculation of geodetic mass balance for three time intervals. The co-registration of all three DEMs with a network of GCP showed the absence of planimetric shifts and helped to reduce vertical bias. Following co-registration and bias correction, the MAD values, which were used to quantify vertical shifts between two DEMs in each pair, were reduced to ±2-3 m on ice-free, stable, exposed terrain with slope less than 30 • ( Table 3). Few outliers were registered in the ablation zones of the Central Tuyuksu and the other glaciers in all three DEMs. The Pléiades DEM had data voids and unrealistic values of elevation change in the 1998-2016 and the 1958-2016 time periods in the accumulation zones of the Central Tuyuksu and the other glaciers in the areas with steep slope and shading. However, the number of affected pixels was small (Figures 3, 4, 6, 7). Overall, all three DEMs showed good performance over the Central Tuyuksu, whereby the sum of changes registered in the periods 1958-1998 (two historic DEMs from photogrammetric surveys) and 1998-2016 (DEMs from 1998 photogrammetric survey and Pléiades) amounted to changes in the period 1958-2016 using DEMs from the 1958 photogrammetric survey and Pléiades (Table 3).
At three out of the other seven glaciers of the Tuyuksu group (Molodezhniy, Kosmodemjanskaya, and Ordzhonikidze), larger errors were observed ( Table 4). There is a discrepancy between the mass balance values calculated for the period 1958-2016 and the two sub-periods, whereby the mass balance of all seven glaciers, calculated from 1958 and 2016 DEMs, is −20.28 ± 2.32 m w.e. and the sum of the mass balances calculated for the two sub-periods is −22.6 ± 2.50 m w.e. ( Table 4). The difference is within the uncertainty margin, and we attribute it to deficiencies in co-registration due to the very steep slopes  characterizing these glaciers. We will aim to resolve this problem by establishing a more even distribution of GCP, most of which are currently located in proximity to the Central Tuyuksu because it was the focus of this study.
The derived negative values of geodetic mass balance indicated a strong downwasting of all the glaciers of the Tuyuksu group (Figures 6, 7), confirming the results of the in situ measurements (Figure 8). While the mean rates of glacier surface lowering, volume loss, and annual rates of geodetic mass balance changed a little from the first to the second assessment period at the Central Tuyuksu, their altitudinal distribution and mass balance gradients did. The period 1998-2016 was characterized by stronger downwasting in the ablation zone and also by a likely increase in accumulation (Figures 4, 5), potentially indicating that mass turnover and water cycle intensified at the Central Tuyuksu in the 21 st century. Similarly, for two other larger glaciers, Molodezhniy and Igly Tuyuksu, the rate of geodetic mass balance did not change, but for the smaller glaciers, whose elevations do not exceed 4,000 m a.s.l., therefore limiting accumulation (e.g., Mametova, Glacier No. 1), the annual rates of geodetic mass balance doubled in 1998-2016 due to a rapid loss of glacier area (Figure 7 and Table 4).
The registered increase in accumulation can be a result of the observed decadal variability in precipitation or a result of using DEM for 2016 when high annual precipitation was observed in Ile Alatau ( Figure 9A). While seasonal snow had melted by late August of 2016 when the Pléiades imagery was obtained and did not contribute to the uncertainty of the DEM (Figure 2), the annual glaciological mass balance in 2015-2016 (0.56 m w.e.) was the second highest positive mass balance on record after 1992-1993 as a result of high accumulation ( Figure 10A). The time stamp of DEM used for the calculation of geodetic mass balance is key to the interpretation of a derived elevation change (Marzeion et al., 2017). In this study, multi-decadal time periods were used, which helped to reduce the influence of the short-term climatic anomalies and reflect low-frequency climatic fluctuations. While no long-term linear trends have been detected in seasonal or in annual precipitation in northern Tien Shan (Unger-Shayesteh et al., 2013;Shahgedanova et al., 2018), a negative anomaly in regional precipitation, which was observed in the 1970-1990s (Severskiy et al., 2016;Shahgedanova et al., 2018), affected the glacier mass balance across Tien Shan (Farinotti et al., 2015). Dyurgerov (2010) assessed the mass balance for the entire Tien Shan and found that the strongest negative mass balance was observed in the 1970s (−0.61 m w.e. a −1 in [1974][1975][1976][1977][1978][1979][1980] and in the mid-1990s (−0.65 m w.e. a −1 in 1994-1997).
A step change in precipitation in the early 1970s (Cao, 1998) was evident from the observations at the Mynzhilky meteorological station located in Kishi Almaty valley at 3,010 m a.s.l. (Figure 9A). Annual precipitation declined from an annual mean of 920 mm in 1958-1970 to 845 mm in 1971-1998 and increased to 892 mm in 1998-2016. This variability in precipitation has been accompanied by an increase in summer temperature at a rate of 0.03 • C a −1 , with a linear trend explaining 42% in temperature variance between 1957 and 2016 ( Figure 9B). While in 1958-1970 the mean ELA and AAR values were 3,742 m a.s.l. and 52%, respectively, during the period 1971-1998, these values were 3,860 m a.s.l. and 35%, reverting to 3,825 m a.s.l. and 41% in 1998-2016 (Figure 10B), which is equal to the all-record means (section "Study Area"). In the latter period, the mass balance was positive in five mass balance years at the Central Tuyuksu ( Figure 10A). This includes the 2015-2016 year when the mass balance was positive both at the Central Tuyuksu and at the Abramov and Golubin glaciers in central Tien Shan (Barandun et al., 2018). Therefore, while the impact of the positive precipitation anomaly of 2016 cannot be excluded, the observed increase in precipitation in the 21st century in comparison with the late 20th century may account for the detected increase in accumulation at higher altitudes (Figures 3, 4, 6, 7). The variability in annual precipitation ( Figure 9A) explains 55% of the total variance in annual mass balance of the Central Tuyuksu, while the variability in air temperature ( Figure 9B) explains 35%.
The values of geodetic glacier mass balance obtained in this study were in agreement with those obtained for the neighboring regions of Tien Shan. Thus, Barandun et al. (2018) reported annual values of geodetic mass balance, derived from the Pléiades and SPOT imagery, for the Golubina, Abramov, and No. 354 glaciers as ranging between −0.30 and −0.42 m w.e. a −1 over the period 2003-2015. Brun et al. (2017 reported annual mass balance values of −0.20 ± 0.08 and −0.40 ± 0.20 m w.e. a −1 , as derived from ASTER DEM, for western and eastern Tien Shan, respectively, for the period 2000-2016. Bolch (2015) reported a geodetic mass balance of −0.42 ± 0.66 m w.e. a −1 derived from multiple imagery, ASTER, and SRTM DEM for the glaciers of the Ala-Archa basin in central Tien Shan. These assessments, except perhaps the less negative mass balance derived by Brun et al. (2017) for western Tien Shan, are close to the rates of change reported for the Central Tuyuksu (−0.35 ± 0.18 m w.e. a −1 ) and for the rest of the Tuyuksu group (−0.43 ± 0.16 m w.e. a −1 ) for 1998-2016, except the very small glaciers positioned at lower elevation (e.g., Mametova and Glacier No. 1) and whose annual geodetic mass balance was more negative at −0.59 ± 0.16 and −0.77 ± 0.16 m w.e.a −1 (Table 4).
Importantly, the estimations of geodetic mass balance derived from the stereo Hexagon KH9 imagery for the longer assessment periods are also in broad agreement with the geodetic mass balance derived for the Tuyuksu group. Bolch (2015) reported a geodetic mass balance of −0.45 ± 0.27 m w.e. a −1 for the Ala-Archa region in 1964-2010. Pieczonka and Bolch (2015) reported an annual mass loss ranging between −0.70 ± 0.24 and −0.36 ± 0.28 m w.e. a −1 for northeastern Tien Shan for 1973-1998. Given the uncertainty ranges, these rates overlap with those derived for the Tuyuksu group where the glaciers were losing mass at a rate of −0.29 ± 0.04 and −0.45 ± 0.04 m w.e. a −1 between 1958 and 2016 (Tables 3, 4). The central values of mass balance estimates (i.e., −0.68 ± 0.40 and −0.70 ± 0.24 m w.e.a −1 ) reported by Pieczonka and Bolch (2015) for eastern Terskey-Alatoo, the region closest to the Tuyuksu group, were more negative than the values derived for the Tuyuksu group possibly because Pieczonka and Bolch (2015) focused on the period with the strongest negative anomalies in precipitation (Figure 9A), resulting in a strongly negative mass balance across Tien Shan (Dyurgerov, 2010).

CONCLUSION
The geodetic mass balance of the glaciers of the Tuyuksu group was calculated for two time steps using two historical DEMs and a high-resolution DEM derived from Pléiades stereo imagery. The results confirmed a continuing mass loss at the Central Tuyuksu, where the long-term glaciological mass balance data were available, and at other the glaciers of the group where there were no direct mass balance measurements. Overall, there was a good agreement between the glaciological and the geodetic mass balance values, and the difference between the two methods was within the uncertainty margin for the Central Tuyuksu. The geodetic mass balance values of −0.35 ± 0.18 and −0.43 ± 0.16 m w.e.a −1 measured at the Central Tuyuksu and the other glaciers of the group in 1998-2016, respectively, were in agreement with those of other studies conducted in Tien Shan. An increase in accumulation in 1998-2016 in comparison with that in 1958-1998 was registered in the accumulation zone of the Central Tuyuksu, where the direct stake measurements are limited and attributed to an increase in precipitation observed in the 21 st century in comparison with the last three decades of the 20 th century. As a result of increasing accumulation, the mean annual geodetic mass balance rates of the Central Tuyuksu and the other glaciers extending over 4,000 m a.s.l. did not change between 1958-1998 and 1998-2016, but they became more negative at the smaller glaciers with a smaller elevation span. Following this application of the Pléiades DEM, subsequent Pléiades imagery and other high-resolution DEMs will be used in the future at regular intervals to assess the regional mass balance and, potentially, the trends in precipitation using glaciers as highelevation gauges.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
VK led the data collection and analysis, assisted by NK and ZU. MS and IS designed the study. MS contributed to the data collection and analysis and produced the initial draft of the manuscript. KW assisted with the production of the digital elevation models and data analysis. All the authors contributed to the discussion of the results and the final draft of the manuscript.