Estimating the longevity of glaciers in the Xinjiang region of the Tian Shan through observations of glacier area change since the Little Ice Age using high-resolution imagery

Abstract

Glacial retreat in response to warming climates in the arid Xinjiang region of northwestern China directly impacts downstream water resources available for local communities. We used high-resolution satellite imagery from 1969 to 2014 to delineate spatial changes in 54 active glaciers in the upper Kaidu River Basin in the Tian Shan as well as their past expanses during the Little Ice Age (LIA). We manually delineated their boundaries based on the interpretation of glacial, geomorphic and topographic features. From the total glacier surface area, we estimated glacier volume and mass. From 1969 to 2014, glacier area decreased by 10.1 ± 1.0 km² (relative loss of 34.2 ± 3.5%) and mass by 1.025 ± 0.108 Gt (relative loss of 43 ± 4.6%). From the LIA maximum (est. 1586 CE) to 1969, relative losses were less (25.7 ± 4.3% area loss and 33.1 ± 5.7% mass loss). Our results indicate that glacier recession is accelerating over time and that the glaciers are currently losing over 1.5 times more relative area than elsewhere in the Tian Shan. Using linear and non-linear projections, we estimate that these glaciers may disappear between 2050 and 2150 CE if climatic warming continues at the same pace.

1. Introduction

The majority of fresh water resources for communities in the Xinjiang region of Northwest China originate as glacial meltwater in the Tian Shan range (Li and others, 2016a). The Xinjiang region is one of the most arid regions in the world, receiving <200 mm of total annual precipitation (Chen, 2014). Studies in the Tian Shan indicate that most glaciers have been retreating rapidly since the end of the Little Ice Age (LIA) and even more rapidly since the 1970s (Aizen and others, 2006; Wang and others, 2011). Mountain glaciers, especially smaller glaciers, react sensitively to changes in climate (Oerlemans, 1994; Haeberli and others, 2007; Wang and others, 2012; Chen, 2014; Li and others, 2014a). Air temperature is the primary factor controlling the rates of glacial melt in the Tian Shan (Aizen and others, 2006; Liu and others, 2006; Chen, 2014). Xinjiang's regional temperature increase of 0.34°C/decade from 1960 to 2010 surpassed global temperature increases over this period (Chen, 2014).

Rising temperatures impact glacial meltwater production, which accounts for over 10% of the river runoff in Xinjiang and has the greatest impact on water resource stability in this region (Li and others, 2010; Li and others, 2014a; Chen, 2014). In the mountains, glacial meltwater is a direct fresh water resource (Immerzeel and others, 2010; Mupenzi and others, 2011; Bolch and others, 2012; Chen, 2014). In the arid lowlands, glacial meltwater recharges aquifers and regulates river runoff, decreasing during colder, wetter periods and increasing during warm, dry periods, thereby supplying much needed fresh water to desert lowlands during the dry season (Immerzeel and others, 2010; Bolch and others, 2012). The continued temperature rise predicted for the region poses an urgent threat to the stability of water resources for Xinjiang communities (Li and others, 2014a).

We investigate past and present glacier changes in the Kaidu River basin of the central Tian Shan mountains. The Kaidu River, the third largest in the Tian Shan, flows into Bosten Lake (Fig. 1), the largest inland fresh water lake in China, which provides water for flood irrigation and protects two populated oases within Xinjiang, directly influencing the Tarim Basin's ecosystem and environment (Chen, 2014; Chipman and others, 2016). The continued thinning and recession of these glaciers will have significant impacts on the hydrology of the drainage basin and negative ramifications for water availability and resource management, especially in the lower reaches of the river system (Mupenzi and others, 2011). The Kaidu River also feeds the Yanqi Basin, which is home to 450 000 people as of 2003 (Li, 2009). Agricultural production is their main source of income, responsible for over four times the amount of income from industrial production (Li, 2009). As a result, sustainable agriculture and water resource management are of major interest in the area.

Fig. 1. (a) Location of the study area (red square) in Xinjiang, location and elevation of meteorological stations (pink triangles), and location of other regional glacier change studies (orange circles) referred to in the discussion and Figure 5: (1) Ala Archa; (2) Akshiirak range; (3) Sary-Jaz and Aksu basins; (4) Tekes watershed; (5) Mt. Bogda. Map credits: Esri, HERE, Garmin, FOA, NOAA, USGS and OpenStreetMap contributors. (b) Study area inset showing the three satellite images used for delineation. The 2014 Pléiades image (topmost layer), displayed in true color, covers 35 glaciers in our study area. The 2012 Worldview-2 image (second layer), displayed in false color, and the 1969 panchromatic Corona image (bottom layer) cover all 54 glaciers in our study area. Yellow stars indicate the location of seven LIA end moraines identified.

Despite the significant influence of Tian Shan glaciers on water resources in Xinjiang, field data and long-term monitoring of the glaciers are extremely limited due to the difficulty of access (Wang and others, 2012; Li and others, 2014a). Thus, most studies that investigate glaciers in this region rely on change detection using remotely-sensed data to characterize glacial change and calculate melt rates. These studies are limited in accuracy, especially on shorter timescales, because they are mostly based on the imagery of low or medium resolution (30 m or coarser). In order to address this inadequacy, we investigate the glacier changes in the warming Kaidu River basin using high-resolution (<3 m resolution) satellite imagery.

We examine temporal changes in the surface area, volume and mass of 54 glaciers in the upper reaches of the Kaidu River based on the delineation of glacier margins using high-resolution images from 1969, 2012 and 2014. In addition to active glacier margins, we also delineate the maximum ice extent during the LIA using moraines and other glacial landforms to provide a longer-term record of glacier changes. In the Tian Shan, the LIA glacial advances began in the 15th century and lasted until the beginning of the 20th century (Liu and others, 2003). We use a range of LIA dates centered around ~430 years ago (~1586 CE) from surface exposure dating of Tian Shan glacier moraines from Li and others (2014a) and Li and others (2016b) to calculate change rates from this time period. The outermost moraines associated with the LIA advance found at the Urumqi River headwaters have been dated using cosmogenic nuclide surface exposure dating of ¹⁰Be based on 31 moraine samples, producing ages ~430 ± 100 years ago (Li and Li, 2014; Li and others, 2016b). These agree with the dates of 412 ± 146 and 428 ± 127 years ago from ¹⁴C dating of the inorganic carbonate coating the moraines (Yi and others, 2004; Li and others, 2011c). Using our measured glacier surface area changes from the LIA maximum through 2014, we estimate when these glaciers will disappear if the current trends in climate continue over the next century.

2. Study area

The Tian Shan is a 2500 km long mountain range formed during the Late Paleozoic Alpine orogeny (Savoskul, 1997; Chen, 2014) that traverses Central Asia, dividing Xinjiang into a northern and southern half (Chen, 2014). In Xinjiang, these mountains span 1700 km and have an average summit elevation of 5000 m a.s.l. with the most prominent summit elevations of 7435 m a.s.l. at Tuomuer Peak and 6995 m a.s.l. at Khan Tengri Peak (Chen, 2014). The orographic effect of the mountain range results in a precipitation gradient from 1500 to 2000 mm a⁻¹ in the northwest to <100 mm a⁻¹ in the southeast (Aizen and others, 1997), where our glaciers are located. Temperatures in Xinjiang range from an average of −20°C in the winter to 23°C in the summer, with plentiful solar radiation (54.4 × 10⁸ to 71.2 × 10⁸ J m⁻² a⁻¹) (Chen, 2014). Thirteen meteorological stations in the region (Fig. 1) have been recording temperature and precipitation data since the 1950s through 2014 and the data were available as monthly averages. Out of these stations, the Bayinbuluke Station is the most proximal in horizontal distance (~65 km) and elevation (2459.4 m a.s.l.) to our study site.

We investigated the changes in 54 glaciers within a 27 km by 18 km area in the eastern central Tian Shan mountains, occurring at a latitude of ~42°N (Fig. 1). These glaciers range in elevation from 3450 m a.s.l. at their termini to 4600 m a.s.l. at their headwalls and had an average glacier area of 0.53 km² in 2012. More than half (57.4%) of the glaciers in the study area are north-facing and the rest are mostly east- or west-facing.

3. Methods

3.1 Data sources and image processing

High-resolution satellite images acquired during the summer months and digital elevation models (DEMs) of the study area are scarce. We obtained our data through a variety of sources including the United States Geological Survey (USGS) EarthExplorer, Digital Globe, Spot Image and Airbus (Table 1). The satellite images used for our visual interpretation have better than 3 m resolution. The three images used for our delineations were a panchromatic Corona image from August 1969, a multispectral Worldview-2 image from August 2012, and a multispectral Pléiades image from September 2014 (Table 1, Fig. 1). Meanwhile, the Shuttle Radar Topography Mission (SRTM) and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEMs have 30 m pixel spacing. A SPOT-7 DEM was created photogrammetrically from high-resolution oblique imagery and has 1.8 m resolution.

Table 1. Image properties for the satellite images and DEMs used in the analysis

We orthorectified the 2012 Worldview-2 image using the SRTM DEM with metadata describing the sensor's viewing geometry. Bilinear interpolation was used when resampling the SRTM DEM to match the high-resolution Worldview-2 image. The 2014 Pléiades image was downloaded as an orthorectified image product from Airbus with 6.5 m CE90 (circular error 90%) geolocation accuracy. The orthorectification was performed using the worldwide Reference3D dataset (http://www.intelligence-airbusds.com/en/4572-pleiades-technical-documents). The 1969 Corona image was scanned from the original panoramic film material at the EROS Data Center at the USGS. Because a model of the sensor's viewing geometry was not available, an empirical Direct Linear Transform (DLT) projection was used to orthorectify the image, using the previously-orthorectified Worldview-2 image as a base. We subset the 1969 Corona image to just the area immediately surrounding the study site to reduce geometric errors. DLT model parameters were computed based on 16 pairs of manually-selected Ground Control Points (GCPs) between the Corona image and the 2012 Worldview-2 image (RMS = 1.38 m). The final orthorectification resulted from the DLT projection model and the SRTM DEM.

The SPOT-7 DEM was created photogrammetrically using stereoscopic panchromatic images acquired off-nadir on 13 July 2015. While the nominal resolution of the SPOT-7 panchromatic band is 1.5 m, the images used here had an effective resolution of 1.8 m due to the oblique viewing geometry, and this 1.8 m spacing was used when producing the DEM. The photogrammetric processing was done using the Automatic Terrain Extraction (eATE) algorithm with dense point matching in ERDAS Imagine, using 31 control points.

3.2 Delineating glacier area

Our glacier delineations are based on the Global Land Ice Measurements from Space (GLIMS) definition of glaciers, which complies with GLIMS standardized methods and the World Glacier Monitoring Service standards (Raup and Khalsa, 2007; Racoviteanu and others, 2009). In accordance with the definitions provided by the GLIMS workshop (Racoviteanu and others, 2009), we included the sections of ice above the bergschrund only when they remained partly connected to the rest of the body. We excluded snowfields and detached stagnant ice as well as nunataks and any exposed ground (Raup and Khalsa, 2007). Our methods comply with the GLIMS recommendations for reducing the amount of error from glacier delineations used for change detection, which include (1) choosing images acquired around the same time of year; (2) calculating area change by determining glacier areas in each analyzed year and subtracting them separately rather than subtracting one map from another using GIS software; and (3) keeping the consistency of upper glacier boundaries and location of internal rocks (Racoviteanu and others, 2009; Paul and others, 2013). We chose images acquired in August and September, at the end of the melt season, to keep the snow cover consistent and minimal.

We visually interpreted the boundaries based on glacial features on and near the glaciers. We differentiated flowing ice and debris-covered ice from snowfields, disconnected stagnant ice, debris and exposed rock (Raup and Khalsa, 2007; Racoviteanu and others, 2009). Crevasses indicating the locations of ice movement and breakage were visible on a few glaciers but were not consistently present, so we relied on the presence of flow bands and color contrasts to differentiate between glacier ice and non-glacial features. When viewed in a false color combination, flow bands are visible as arcuate white or cyan features (Fig. 2a). Our delineation process differed slightly between years of analysis due to the differences in resolution, spectral bands and spatial extent of each of the images. For example, geomorphic features such as boulders and peaks of exposed bedrock in the middle of glaciers were highly visible in the 2014 Pléiades image but were indistinguishable in the 2012 Worldview-2 image. Visual analysis of the panchromatic Corona image from 1969 relied on the interpretations of contrast and intensity. Due to difference in spatial extent of the 2014 Pleiades image compared to the other images used, only 35 of the 54 glaciers were delineated for 2014. The glaciers are referred to by a Glacier ID consisting of a letter or combination of letters (e.g. Glacier ZA). They will be referred to by their Glacier IDs from here onwards.

Fig. 2. Geomorphic and glacial features visible in the Glacier A valley (42.538° N, 83.690° E) in the (a) 2012 Worldview-2 Image with the 2012 outline in purple and (b) the 3-D perspective view of the same area. (1) Three characteristic end moraines corresponding to LIA ice extent. 1A indicates the LIA maximum position used for delineation. (2) Trim lines aligning with the most distal LIA end moraine that informed LIA lateral boundary placement. (3) Flow bands. (4) Nunatak. (5) Topography used to determine the direction of ice flow near ice divides.

3.2.1 Delineations on 1969 Corona imagery

Panchromatic satellite images cover a wider range of wavelengths than most multispectral bands, which allows them to capture a wide dynamic range. Flow patterns in the ice and boundaries between ice and rock were more prominent than in the other images, due to the contrast provided by this dynamic range. Cloud cover challenged our interpretation of the 1969 Corona image, obscuring several glacier boundaries that are visible for the other years. If the ice boundary of a glacier was too obscured by clouds to be delineated, the glacier was not given a Glacier ID and was excluded from the change detection analysis.

3.2.2 Delineations on 2012 Worldview-2 imagery

The August 2012 Worldview-2 image was viewed in false color. We kept the upper boundaries of the 1969 outlines consistent for the 2012 delineations and adjusted them to account for minor changes in the upper extents and any inconsistencies in registration in accordance with GLIMS standards. Then, we modified the lower reaches of the glaciers (where there is the most retreat) to match their 2012 ice extents based on the glacial features visible in the imagery. We supplemented the visual interpretation of ice extent in 2012 with three-dimensional (3-D) models of the glaciers and their valleys. The 3-D visualization of the topography aided interpretation near the glacier termini and at flow divides, where the accumulation zones of neighboring glaciers coincided (Fig. 2b).

Moraines and other depositional features indicative of the LIA maximum were easily identified in the 2012 Worldview-2 image, so we primarily used this image to determine the LIA glacier extents. We placed the past glacier boundaries at the location of depositional features such as lateral and end moraines and outwash fans. Liu and others (2003) identified three characteristic end moraines that distinguish the LIA extents from other glaciation events. Our delineations of the LIA glacier boundaries are based on the location of the LIA moraines that are clearly visible in the 2012 and 2014 images (Fig. 2a). These end moraines were identified for Glaciers A, D, O, N, ZC, ZI and ZN (locations shown in Fig. 1b). These moraines occurred at the trim line in the 2012 Worldview-2 image (Fig. 2a), so we used the trim lines to inform our placement of the LIA outlines where moraines were not found.

3.2.3 Delineations on 2014 Pléiades imagery

The Pléiades image from September 2014 is of higher resolution (0.5 m resolution) than the other images. Smaller features such as boulders on top of the ice and glacial erratics were identifiable in the Pléiades image and informed our placement of the glacier boundaries. We maintained the upper limits of the 2012 extents for the 2014 boundaries, since most of the glaciers were snow-covered at their upper reaches. Then, we re-delineated the lower limits of the boundaries to match the 2014 glacier extents. Since snow covered the glacier surfaces but was absent from non-glacial terrain, we placed the lower limits of the boundaries at the snow extent and then adjusted them to remove debris and exposed rock.

3.3 Determining glacier surface area changes and uncertainty

The total glacier area is calculated as the sum of the delineated glacier areas. The total area and relative area changes, as well as annual rates of area change, were calculated for the following time periods: LIA to 1969, 1969–2014 and 2012–2014 for the 35 glaciers measured across all three images. In addition, we calculated the change from 1969 to 2012 for all 54 glaciers. While image processing errors (from georeferencing, orthorectification, etc.) influence the uncertainty of glacier delineations (Racoviteanu and others, 2009; Paul and others, 2013, 2017), the errors from the human interpretation of the features in the images will be much greater than any error from image processing. The primary source of error for our calculated glacier areas is the uncertainty from analyst interpretation.

We conducted a multiple digitization experiment to estimate the uncertainty of our glacier areas following methods described by previous work that estimates uncertainties from human error (Hall and others, 2003; Paul and others, 2013; Frey and others, 2014; White and Copland, 2018). To quantify the error from interpretation, we identified a representative sample of six glaciers (A, D, P, I, ZN, ZL) in our study area which varied in size from 0.2 to 1.8 km². A single analyst delineated the glaciers five times (Fig. 3). Additional outlines from the repeat delineations are shown in the Supplementary Material (Fig. S1). We extracted the glacier areas from each digitization experiment and calculated the 95% confidence interval for the sum of the glacier areas in each experiment, which represents the uncertainty in our total surface area measurements. These uncertainties were propagated through all calculations using the areas following standard error propagation procedures.

Fig. 3. Repeat delineations for the uncertainty analysis shown for (a) Glacier A in 2012 (1σ = 0.0231 km²), (b) Glacier ZN in 2012 (1σ = 0.0046 km²), (c) Glacier A in 2014 (1σ = 0.0235 km²) and (d) Glacier P in 2014 (1σ = 0.0706 km²). Sigma values represent the std dev. in surface areas calculated from the five repeat delineations for each glacier.

3.4 Glacier slope, aspect and hypsometry

All glacier elevations, slopes and aspects were derived from the highest-resolution DEM, the 1.8 m resolution SPOT-7 DEM. The values within each glacier's bounds were extracted using the 1969 glacier outlines. Following Jiskoot and others (2009), the aspect values were binned into eight 45° intervals corresponding to the compass rose directions. The dominant aspect for each glacier shown in Table 2 corresponds to the bin with the greatest percentage of counts. Following a similar approach, the slopes in each glacier's boundaries were binned into 100 intervals. The most representative slope for each glacier was the bin with the greatest percentage of counts. We also report the maximum elevation and minimum elevations in the LIA, 1969, 2012 and 2014, as well as the change in minimum glacier elevation from the LIA to 2014 (Table 2).

Table 2. Surface area, direction of the dominant glacier aspect, minimum and maximum elevations (Z _min and Z _max), terminus elevation change (ΔZ _min), Hypsometric Index (HI) and class, and area change data for the 35 glaciers measured from the LIA to 2014

The HI classes are: very-bottom heavy (VBH), bottom-heavy (BH), equidimensional (EQ), top-heavy (TH) and very top-heavy (VTH). For glaciers that split into multiple smaller glaciers over time (e.g. Glaciers B and C), the ΔZ _min from LIA to 2014 is calculated by subtracting the average of all 2014 Z _min values by the Z _min in the LIA. This approach for split glaciers is used for relative area change calculations as well.

We determined the hypsometric curves for each of the glacier using the distribution of elevations in each glacier's boundaries, binned into 5000 equal intervals. The cumulative sum of the surface area represented by each of the elevation intervals was used to generate the curves of elevation with respect to the normalized area. Additionally, we calculated the Hypsometric Index (HI) using the expression from Jiskoot and others (2009) (Eq. 1) in order to investigate the area changes based on the hypsometric classification of the glaciers.

(1)$${\rm HI} = {\rm \;}\displaystyle{{Z_{{\rm max}}-{\rm \;}Z_{{\rm med}}} \over {Z_{{\rm med}}-{\rm \;}Z_{{\rm min}}}}\,{\rm and}\,{\rm HI} = {\rm \;}\displaystyle{{-1} \over {{\rm HI}}}\;\,{\rm if}\,0 \lt {\rm HI} \lt 1.$$

These HI values were used to differentiate the glaciers into five hypsometric categories: very top-heavy (HI < −1.5), top-heavy (−1.5 < HI < −1.2), equidimensional (−1.2 < HI < 1.2), bottom-heavy (HI > 1.2) and very-bottom heavy (HI > 1.5), as done in Jiskoot and others (2009). Table 2 shows the HI values as well as the hypsometric categories (i.e. HI class). Uncertainties reported in Table 2 are the 95% confidence interval for the values determined from the sensitivity of the measured values to different source DEMs and the number of bins used in histogram calculations.

3.5 Total glacier area to volume and mass

Volume–area scaling relationships have been widely used in glaciology and have a strong theoretical basis (Bahr and others, 1997, 2015). The power law scaling relationship for glaciers was derived from the models of typical glacier geometry (Bahr and others, 2015):

(2)$$V{\rm \;} = {\rm \;}cA^{\rm g}\comma \;$$

where g is ~1.36 (Frey and others, 2014). We use the constants determined by Chen and Ohmura (1990), c = 28.5 and g = 1.357, to calculate our glacier volume. These coefficients were determined from a regression of measurements of 63 glaciers in different regions which yielded a correlation coefficient of r ² = 0.96 (Chen and Ohmura, 1990).

From glacier volume, we use a glacier ice density value, 850 ± 60 kg m⁻³ (Huss, 2013), to convert volume to mass. This glacier density value accounts for the variation in density vertically throughout the glacier based on a firn compaction model, so it is lower than the density of pure ice. Wei and others (2015) used this value to estimate glacier mass balance from volume estimates of glaciers in Bangong Co Basin on the Tibetan Plateau. The same density value for conversion was used by Berthier and others (2016) and Brun and others (2017) in studies determining glacier mass balances in the Mont-Blanc area and High Mountain Asia, respectively.

4. Results

4.1 Glacier delineations and individual glacier changes

Figure 4 shows the glacier delineations from the LIA through 2014. Additional figures showing the close-ups on outlines from a single glacier or group of glaciers are shown in the Supplementary Material (Figs S2–S6). The surface area, dominant aspect and minimum and maximum elevation for each of the 35 glaciers delineated across all dates are shown in Table 2, along with minimum elevation changes and relative area changes.

Fig. 4. LIA, 1969, 2012 and 2014 glacier outlines with Glacier ID labels. Note the variation in glacier shape and orientation as well as the retreat in their ablation zones. Glaciers without 2014 delineations lay outside the spatial extent of the 2014 Pléiades image.

4.2 Total glacier area, volume and mass changes

We report area changes from the LIA-2014 that include the total area from the 35 glaciers in Table 3. We did, however, calculate the surface area, volume and mass change for all 54 glaciers for the period 1969–2012. For the 54 glaciers, the total area change was −12.2 ± 1.0 km², which equated to a relative area change of −30.3 ± 2.5% and a rate of change of −0.22 ± 0.02 km² a⁻¹ (−0.55 ± 0.02% a⁻¹) from 1969 to 2012. This relative change is comparable to the relative change we determined for the group of 35 glaciers from 1969 to 2014 (−34.2 ± 3.5%). For the 35 glaciers, the percent glacier area, volume and mass changes were greater over the 45-year period from 1969 to 2014 than over the ~400-year period from the LIA to 1969. The percent change over just 2 years from 2012 to 2014 is one-fifth of the percent change over 45 years from 1969 to 2014. Rates of change in area, volume and mass increase over time and the average rate of change in 2012–14 (−2.1 ± 1.1% a⁻¹ area loss) is more than 20 times faster than it was from the LIA to 1969 (−0.05 ± 0.008% a⁻¹ area loss). Table 4 lists the total surface area, volume and mass of the glaciers delineated in the LIA, 1969, 2012 and 2014.

Table 3. Glacier surface area, volume, and mass changes for 35 glaciers from the LIA to 1969, 1969–2014 and 2012–14

ΔS is the surface area change, ΔV is the volume change and ΔM is the mass change. The percent change in mass is the same as the percent change in volume, so it is not duplicated.

Table 4. Number of glaciers delineated (N) and the total glacier surface area, volume and mass calculated for glaciers in the LIA, 1969, 2012 and 2014

5. Discussion

5.1 Glacier delineations

Processes that deviated from the GLIMS recommended methods introduced additional uncertainty to the analysis. GLIMS recommends using the same type of satellite image for each analyzed year and using the same analysis methods for each delineation set (Racoviteanu and others, 2009). Due to changes in the suite of available high-resolution satellites over time, we relied on three different types of satellite images. Visual similarities between snow and flowing ice in the imagery complicated the identification of snowfields adjacent to the glacier, snow on the glacier and snow on top of exposed rock in the middle of glaciers. Supraglacial debris also contributed to the ambiguity of certain glacier outlines. In addition, delineation of past glacial extents was limited by the presence of moraines found in the glacial valleys. While the three characteristic LIA end moraines were identified for seven glaciers, a majority of the glaciers exhibited only one or two of the moraines in the sequence. In order to determine the LIA extents for those remaining glaciers, we relied on the location of the trim lines, which coincided with the most distal LIA end moraines. While this allowed us to draw outlines connecting end and lateral moraines as well as infer the boundary when there were no nearby moraines, this contributes to the uncertainty of our LIA outlines. The moraine sequences identified correspond to the LIA maximum date of 1586 CE ± 100 years from Li and Li (2014) and Li and others (2016b). To account for the possibility that these features could correspond to a more recent LIA advance such as 1797–1865 CE (Liu and others, 2003), the range of the LIA dates that could be associated with these delineations is from 1486 to 1865 CE. This larger range of dates was used to project retreat rates and estimate the timing of these glaciers' disappearance. These uncertainties in delineation are reflected in our multiple digitization experiment. Interpretation of the LIA extent yielded the greatest error.

Despite these challenges, the uncertainty of all our delineated glacier areas from analyst interpretation is small compared to the glacier areas calculated. The 95% confidence interval in area calculated from multiple delineations is ±3.7% of the surface area. The std dev. of these repeat delineations ranged from 2.2 to 4.2%. This std dev. is less than the error determined by the GLIMS remote-sensing glacier digitization workshop, which yielded std dev. of 2.6–5.7% for manual delineations (Paul and others, 2013).

5.2 Rates of glacier change

Our rates of glacier surface area loss increase over time, suggesting that glacier shrinkage is accelerating. The rate of area loss from the LIA to 1969 period was from 0.03 ± 0.005 to 0.15 ± 0.005 km² a⁻¹, calculated using LIA dates ranging from 1486 to 1865 CE. From 1969 to 2014, the rate of area loss was 0.22 ± 0.02 km² a⁻¹ and from 2012 to 2014 it was 0.60 ± 0.33 km² a⁻¹, at least four times and up to 20 times faster than the rate during the LIA to 1969 period (Table 3). Aizen and others (2006) observed a similar trend in the acceleration of surface area loss for glaciers in the Akshiirak glacierized area in the central Tian Shan, where the rates of area loss of 0.53 ± 0.0004 km² a⁻¹ for the 1943–1977 period increased to a rate more than twice as fast (1.35 ± 0.001 km² a⁻¹) during the second period 1977–2003 (Aizen and others, 2006). We also provide the estimates of glacier volume and mass changes, as they respond more directly to climate forcing than surface area (Bolch and others, 2012; Pieczonka and Bolch, 2015). Following the trend in the rates of surface area loss, glacier volume and mass loss are also accelerating over time. The rate of glacier volume loss from the LIA to 1969 (0.0036 ± 0.0007 km³ a⁻¹) increased to 0.0268 ± 0.0028 km³ a⁻¹ for the 1969–2014 period and was more than twice as fast (0.0653 ± 0.0364 km³ a⁻¹) for the recent biennial period from 2012 to 2014. Regional climate changes, mainly increases in air temperature, may be driving these accelerations in glacier shrinkage.

From 1957 to 2014, the air temperature record from the most proximal meteorological station to our study site (Bayinbuluke) indicates a rate of warming of 0.22°C/decade and a summer temperature increase of 0.26°C/decade. The average rate of temperature change for 13 stations in Xinjiang (see Fig. 1) over that same time period is 0.23°C/decade overall and 0.14°C/decade for summer temperatures. This rate of change in Xinjiang is lower than the other reported climate trends in all of northwestern China (0.32°C/decade, Chen, 2014), but on par with air temperature trends for China as a whole since the 1980s (0.22°C/decade, Li and others, 2011b). The warming in northwestern China is approximately twice the global rate of warming of 0.26°C/decade since 1979 (IPCC, 2013). If the rise of regional air temperatures continues, glacier melt at high elevations may accelerate further and hasten desertification processes across the northwest (Aizen and others, 2006).

5.3 Factors influencing rapid surface area loss

The annual rate of surface area loss of all glaciers in our study area combined is 0.76 ± 0.08% a⁻¹ from 1969 to 2014 and is higher than the range of 0.1–0.3% a⁻¹ from the 1970s–2000s reported for glaciers with a similar size range elsewhere in the Tian Shan (Li and others, 2011a). The relative glacier area loss in our study area from 1969 to 2014, 34.2 ± 3.5%, is also greater than the relative area changes determined for other glaciers in the Tian Shan (e.g. Aizen and others (2006), Li and others (2011a), Faronotti and others (2015), Osmonov and others (2013), Xu and others (2015), etc.). Reported values for relative glacier area loss elsewhere in the Tian Shan and calculated over a similar time period were lower, ranging from 10.5 to 26% (Fig. 5). Our calculated relative change in glacier area is 1.5 times greater than the values determined by Xu and others (2015) and is even greater than the glacier area losses calculated in the other studies, suggesting that the glaciers in our study area are losing more relative area than other glaciers in the Tian Shan (Fig. 5). Differences in average glacier size, elevation and aspect in this study area versus other study areas could also be responsible for our greater relative area losses. We examine the effect of these factors on our calculated percent area losses.

Fig. 5. Reported relative surface area changes (red points) and the time interval (grey bars) of glacier change studies in the Tian Shan. Darker red points represent greater relative area losses. Each study is labeled by author(s), study region and number of glaciers (n) in the sample. The black error bars represent the uncertainty in glacier area change for the three studies that report an uncertainty value.

5.3.1 Initial size

Average glacier size between study areas is not likely to be a primary factor influencing the inter-study differences in relative area loss. Smaller glaciers tend to lose more relative area as they recede faster (Aizen and others, 2007; Kehrl and others, 2014; Li and Li, 2014; Li and others, 2014b; O'Neel and others, 2019). In the 1960s, average glacier size in different regions of the Tian Shan ranged from 0.43 to 1.03 km² (Li and others, 2011a), while our average glacier size was 1.79 km² in 1969 and thus larger on average than other glaciers around the Tian Shan during the 1960s. Therefore, we would expect the larger glaciers in our study site to have lost less relative area than the smaller glaciers elsewhere in the Tian Shan. However, they receded more than other Tian Shan glaciers. While average glacier size is not likely to be responsible for the inter-study differences in relative area loss, initial size does affect area loss on a glacier-by-glacier basis. For the glaciers in this study, the mean percent area loss for glaciers larger than 1 km² (28.9%) was significantly different at the p = 0.005 level from the mean percent area loss for glaciers smaller than 1 km² (49.8%), based on a one-way analysis of variance (ANOVA) test. The smaller glaciers experienced 20.9% more relative area loss on average than the larger glaciers (Fig. 6a). Li and others (2014b) and Li and Li (2014) similarly found an inverse relationship between initial glacier area and relative area loss for other glaciers in the Tian Shan.

Fig. 6. (a) Initial glacier size versus percent area loss from 1969 to 2014. Dashed line indicates the separation of smaller glaciers (<1 km² in size) versus larger glaciers (>1 km² in size) for which the one-way ANOVA test was performed. (b) Dominant aspects versus percent area losses.

5.3.2 Aspect

Glacier aspect may be influencing our glaciers' exceptional percent area losses. At mid-latitudes in the Northern Hemisphere, north-facing glaciers receive less solar radiation than south-facing ones (Sakai and others, 2002). As a result, south-facing glaciers in High Mountain Asia lose more relative area than north-facing glaciers (Nainwal and others, 2008). Li and Li (2014) found that 90.3% of 487 glaciers in the central Tian Shan are north-facing. In our study area, only 57.4% of the glaciers are north-facing and thus nearly half are exposed to greater solar radiation. However, further work must be done to connect solar radiation changes to total glacier area loss in this region. On a glacier-by-glacier basis, Li and Li (2014) found that glacier aspect was significantly correlated (at the p = 0.001 level) with area change. We similarly found a statistically significant relationship between our glaciers' dominant aspects and their area losses since 1969. The mean area loss for north-facing glaciers was significantly different at the p = 0.01 level from the mean area loss for south-facing glaciers based on a one-way ANOVA test. South-facing (S, SE and SW) glaciers underwent significantly different (p = 0.0061) losses than north-facing (N, NE and NW) glaciers (Fig. 6b). This is independent of influence from glacier size, as south-facing glaciers and north-facing glaciers did not have significantly different mean glacier sizes (p = 0.9). Mean area loss from 1969 to 2014 for south-facing glaciers was 60.8 versus 37.0% for north-facing glaciers. The difference in solar radiation experienced by glaciers with dominantly south-facing aspect is likely influencing area losses for individual glaciers as well as greater relative area loss at the basin-level.

5.3.3 Average elevation, hypsometry and surface slope

Studies in the Tian Shan, and worldwide, indicate that alpine glaciers recede more rapidly at lower elevations (Li and others, 2011c; Li and Li, 2014). Regional meteorological stations also indicate greater rates of summer temperature warming at lower elevations over the last half-century. From 1957 to 2014, summer air temperatures from the six stations located below 1000 m a.s.l. (Jinghe, Wusu, Yining, Kumishi, Luntai and Kuerle) increased by 0.45°C/decade on average while the Bayinbuluke station, located at 2459.4 m a.s.l., recorded a smaller 0.26°C/decade summer air temperature increase over the same period. The elevations of the glaciers in this study range from 3600 to 4200 m a.s.l., which is lower than many of the other glaciers in the Tian Shan studies shown in Figure 5. The glaciers in the Akshiirak glacierized massif and the Ala Archa glacierized basin in the Tian Shan have higher altitudinal extents of 3600–5000 and 3300–4800 m a.s.l., respectively, and lost less surface area (15.2%, Aizen and others, 2006) compared to our glaciers. In addition, the HI values calculated (Table 2) indicate that a majority (79%) of the glaciers in our study are bottom-heavy or very bottom-heavy, which means that more of their area is distributed at lower elevations where there are warmer temperatures and greater rates of warming, potentially contributing to their large relative area loss compared to the other studies.

While very bottom-heavy glaciers underwent greater percent area losses than very top-heavy glaciers in our study, this difference is not statistically significant. Differences in their relative area loss may be driven by other factors, such as their differences in size or surface slope, as the very top-heavy glaciers are also smaller and have steeper surface slopes (Fig. 7). Recent work has shown that surface slope is related to response time to climactic forcing, with alpine glaciers with steeper slopes responding more rapidly to climatic changes (Zekollari and others, 2020). Analysis of relative area changes, sizes and surface slopes indicate that there is no significant correlation between dominant surface slope and size for the glaciers in this study. Meanwhile, there is a weak positive correlation between surface slope and area loss (Fig. 7d). The steeper glaciers in the study may be more sensitive to regional climatic warming and may lose their mass sooner.

Fig. 7. (a) The 35 glaciers shown by their Hypsometric Index (HI) class: very bottom-heavy (VBH), bottom-heavy (BH), equidimensional (EQ), top-heavy (TH) or very top-heavy (VTH). (b) Hypsometric curves based on normalized area for a glacier from each of the following three classes: VBH, EQ and VTH. (c) Boxplot of the average relative area changes for glaciers in each of the HI classes. (d) Surface slope versus relative area change for each glacier showing a weak positive correlation between the two variables.

5.4 Calculating a time to future disappearance

5.4.1 Linear projections

We calculated an estimate of the time for the glaciers to disappear completely using our glacier areas from the LIA through 2014. First, we calculate a disappearance year based on modern (1969–2014) glacier area change rates. We fit linear equations to the total glacier areas from 1969, 2012 and 2014 to extrapolate a conservative estimate for the end dates of our glaciers (Fig. 8). This model excludes the LIA data point, taking only the modern rate of change to determine the disappearance date. In order to capture the range of possible disappearance years due to uncertainty of the area change measurements, we performed a Monte Carlo simulation of 2000 possible linear fits based on randomly generated area values within one std dev. of the observed value (Fig. 8). The best-fit line is then determined by averaging the coefficients calculated for the 2000 simulated best-fit lines. Our linear model has limited resolving power because it is based on just 3 years of analysis.

Fig. 8. Linear fits to total glacier areas for 35 glaciers from 1969, 2012 and 2014. Gray lines are linear fits generated using the 2000 variations of glacier surface area values, yielding a range of disappearance years from 2078 to 2154 CE. The best-fit line (R ² = 0.9964) shown in red yields a disappearance year of 2106 CE.

The linear fits from the simulations yield a range of disappearance years from 2078 to 2154 CE, with a mean disappearance year of 2106 CE, based on the 1969–2014 rate of melt (−0.22 ± 0.02 km² a⁻¹). These estimates are in accordance with the previous estimates of the future disappearance of Himalayan and Karakoram alpine glaciers which indicate that 3000–13 000 glaciers (64–73% of the total glaciers in 1985 CE) in High Mountain Asia may disappear as early as 2035 CE if the melt rate remains at current levels (Cogley, 2011). To account for potential increases in melt rate due to increasing regional temperatures, we also calculate a disappearance year based on a doubled constant melt rate. With an area loss rate that is double the modern rate of area loss from 1969 to 2014, the glaciers will likely disappear in 2057 CE, nearly 50 years before the predicted disappearance year based on current melt rates (Fig. 9).

Fig. 9. The best-fit non-linear (R ² = 0.9996) and linear models for glacier disappearance based on total glacier surface area for 35 glaciers. The linear model based on the 1969–2014 (modern) rate of area loss indicates a disappearance year of 2106 CE. The various non-linear curves generated from points within the confidence intervals for the data yield disappearance years between 2076 and 2057 CE. The latter coincides with the projection using a doubled linear melt rate, which also indicates glacier loss in ~2057 CE.

To explore the difference in potential disappearance years for the smaller glaciers (<1 km²) in our study area versus larger glaciers (>1 km²), which experience significantly greater area losses, we also calculate potential disappearance years for each of these subgroups, projecting their current rates of change into the future. For glaciers smaller than 1 km², the rate of change in surface area from 1969 to 2014 was −1.1% a⁻¹, yielding a disappearance year of 2048 CE. For glaciers larger than 1 km², which lost surface area at a lower rate of −0.64% a⁻¹ from 1969 to 2014, the disappearance is calculated to occur in 2130 CE, much later.

5.4.2 Non-linear projections

If we consider changes in ice dynamics over time, a linear model does not capture long-term changes in the rates of melt. Previous work, as well as ours, indicates that the Tian Shan glaciers have begun to retreat faster in the last four decades (Aizen and others, 2006). The inclusion of the LIA data point is necessary to improve the accuracy and resolving power for a non-linear fit. The best non-linear fit to the glacier areas yields the same disappearance year estimate as the doubled linear melt rate, indicating that the glaciers may disappear as early as 2057 CE, <40 years from the present (Fig. 9). The disappearance years for the glaciers in our study from all projections range from 2057 to 2154 CE, with the glaciers smaller than 1 km² potentially disappearing even earlier in 2048 CE. Smaller glaciers respond more rapidly to climatic warming, so glaciers will melt at faster rates as they thin and their total volumes decrease (Kehrl and others, 2014; Li and others, 2014b). These estimates do not account for the possibility of glacier stabilization at lower volumes as the glaciers' ablation zones retreat to higher elevations (DeBeer and Sharp, 2009; Huss and others, 2012). While models that incorporate future climatic forcing do account for this relationship (Huss and Hock, 2015, 2018; IPCC, 2019; Zekollari and others, 2019; Rounce and others, 2020), climate modeling is beyond the scope of our study. Our projections represent the lower bounds on the timing of glacier disappearance.

5.4.3 Hydrological response

In the case that the glaciers in our study area melt fully, a total of 2.62 ± 0.04 km³ of glacier ice (as of 2012) would be lost. The annual decrease in water volume stored in ice is equivalent to annual excess water discharge from glacier meltwater (Brun and others, 2017), which increases the base flow of the Kaidu River annually. Glacier meltwater is a dominant source of the Kaidu River's runoff, comprising 41–52% of the total runoff (Haiyan and others, 2018). The annual discharge from meltwater will vary over time, increasing initially with glacier loss, and then will begin to decrease at a critical point (peak water) as the glaciers shrink in area and volume (Braun and others, 2000; Hock and others, 2005; IPCC, 2019). Our projections of glacier disappearance between 2057 and 2106 CE aligns with the projections of the timing of peak water. In High Mountain Asia, the timing peak water for most smaller glaciers (0–10 km² in size) such as the glaciers in our study is projected to be between 2030 and 2040 CE for both RCP2.6 and RCP8.5 scenarios, according to a global glacier runoff model based on 14 General Circulation Models (Huss and Hock, 2018; IPCC, 2019). This projection of timing of peak water to occur within the next two decades is a major concern for the nearby communities. The decline in runoff following peak water could ultimately result in water shortages during the dry season and considerably impact longer-term, water resource management decisions (Chen, 2014). The decline in the Kaidu River's runoff may push communities to rely more on groundwater resources. The glacial meltwater in the upper Kaidu River basin directly affects the availability of water downstream seasonally for its potential uses, such as for human consumption, irrigation, industry, etc. Losses in glacier volume also influence the seasonal timing and variability of streamflow (Shen and others, 2018). These hydrological changes will influence future water management decisions in the region regarding intra- and inter-basin water transfers, irrigation and dam construction (Chipman and others, 2016). Continued monitoring of glacier changes will be critical for determining future impacts on hydrology.

6. Conclusion

Delineations of glacier surface areas in the upper Kaidu River Basin indicate that the glaciers in this region are losing more relative area (34.2 ± 3.5%) than other glaciers in the Tian Shan (12–27%). This greater area loss may be due in part to local or regional factors including average glacier orientation and the relatively low elevations of the glaciers in the study area. Initial glacier size and dominant glacier orientation influence area losses on a glacier-by-glacier basis, with smaller glaciers and south-facing glaciers losing 20.9 and 23.8% more surface area than their counterparts, respectively. Glacial retreat is accelerating across our study area. The recent rate of glacier surface area loss (−0.60 ± 0.33 km² a⁻¹ from 2012 to 2014) is nearly triple the rate of surface area loss over the last half century (−0.22 ± 0.02 km² a⁻¹ from 1969 to 2014). Using two simple forward projections, one linear and one non-linear, we project that the glaciers may disappear as early as 2050 CE. The timeline of total glacier melt will impact water security for the communities residing in the lowlands of the Kaidu River watershed.