Projected climate change and its impacts on glaciers and water resources in the headwaters of the Tarim River, NW China/Kyrgyzstan

Glacierised river catchments are highly sensitive to climate change, while large populations may depend on their water resources. The irrigation agriculture and the communities along the Tarim River, NW China, strongly depend on the discharge from the glacierised catchments surrounding the Taklamakan Desert. While recent increasing discharge has been beneficial for the agricultural sector, future runoff under climate change is uncertain. We assess three climate change scenarios by forcing two glacio-hydrological models with output of eight general circulation models. The models have different glaciological modelling approaches but were both calibrated to discharge and glacier mass balance observations. Projected changes in climate, glacier cover and river discharge are examined over the twenty-first century and generally point to warmer and wetter conditions. The model ensemble projects median temperature and precipitation increases of + 1.9–5.3 °C and + 9–24%, respectively, until the end of the century compared to the 1971–2000 reference period. Glacier area is projected to shrink by 15–73% (model medians, range over scenarios), depending on the catchment. River discharge is projected to first increase by about 20% in the Aksu River catchments with subsequent decreases of up to 20%. In contrast, discharge in the drier Hotan and Yarkant catchments is projected to increase by 15–60% towards the end of the century. The large uncertainties mainly relate to the climate model ensemble and the limited observations to constrain the glacio-hydrological models. Sustainable water resource management will be key to avert the risks associated with the projected changes and their uncertainties.


Introduction
The Tarim River in the Xinjiang Uighur Autonomous Region of NW China is the country's largest endorheic river basin and is home to approx. 10 Million people. Runoff is generated in the glacierised Tien Shan (Aksu River), Pamir, Karakoram (Yarkant R.) and Kunlun Mountains (Hotan R.) encircling the Taklamakan Desert and is vital to the river oases that live from irrigation agriculture. A decline in upstream runoff may have severe consequences for the population, agriculture and ecosystems downstream, threatening the very existence of the river as it has done in the past (Thevs 2011). As most mountainous regions, the highly glacierised and largely uninhabited headwaters of the Tarim River are vulnerable to climate change, as higher temperatures endanger the glaciers as natural reservoirs for downstream communities (Barnett et al. 2005;Immerzeel et al. 2020;Pritchard 2019). At the same time, competition for river water has dramatically increased over the past 30 years with the rapid expansion of irrigated agriculture for the production of cotton, fruit trees and subsistence farming. While an increase in headwater discharge has been observed Krysanova et al. 2015), downstream discharge has declined in line with a four-fold increase in irrigation area over the past 40 years (Thevs 2011).
The continental climate produces scarce precipitation mainly in summer when it coincides with the glacier melt peak (Aizen et al. 1995). Meteorological measurements are sparse and mainly in valleys or at the fringes of the Taklamakan Desert (Tao et al. 2011;Krysanova et al. 2015), which introduces an underestimation in catchment-wide precipitation interpolations. Comparing various datasets and using glacio-hydrological modelling, Wortmann et al. (2018) show that the best available dataset (APHRODITE, Yatagai et al. 2012) has to be corrected by factors of 1.2-4.3 at the catchment average (more at higher elevations) to be consistent with discharge records and glacier mass balances estimates. Similar findings are reported for other areas in the Tarim basin and its surroundings (Duethmann et al. 2013;Immerzeel et al. 2015;Tong et al. 2014).
Several studies about the Tarim River basin investigate the observed trends in discharge and its anticipated climatic (Xu et al. 2004;Chen et al. 2006;Yaning et al. 2009;Li et al. 2020), as well as water management drivers (Song et al. 2002;Liu and Chen 2006;Tang et al. 2007;Hao et al. 2008;Wu et al. 2010). A small number of contributions have addressed changes of glaciers or glacier runoff under projected climate change, but only at a global or regional scale (Kraaijenbrink et al. 2017;Huss and Hock 2018;Rounce et al. 2020) or in parts of the Tarim Basin. Duethmann et al. (2016) examine three climate change scenarios in the Aksu headwaters and project glacier area to decrease by 32-90% over the twenty-first century with an initial increase and an eventual decrease in discharge. Zhang et al. (2012) assess glacier runoff, mass balance and area over the first half of the twenty-first century using a simple degree-day glacier melt model and three delta-change scenarios. However, no systematic assessment of future changes under common climate change scenarios, considering the extensive glacier cover and their heterogenous response to climate change in the headwaters, has been conducted to date. This is largely because of a complex hydrology and the data scarcity or restrictions on the limited data availability.
Assessments of glacier area and mass balance in the Tarim headwaters are crucial for adequate hydrological modelling. Advances have been made; for example, different glacier inventories exist covering large parts (2 nd Chinese glacier inventory, Guo et al. 2015) or the whole basin (the Randolph Glacier Inventory consisting of different sources, RGI consortium 2017). Investigations of observed glaciological changes at the catchment scale have been conducted for the Aksu headwaters by Pieczonka and Bolch (2015), using multitemporal satellite images and digital terrain models. Although the existing mass balance assessments for the Yarkant and Hotan catchments only cover the period after 2000 (Brun et al. 2017;Shean et al. 2020;Hugonnet et al. 2021) or not the entire catchment area (Holzer et al. 2015;Zhou et al. 2018;Bhattacharya et al. 2021), they still provide evidence for the calibration of glacio-hydrological models.
This study aims to provide a first systematic climate change impact assessment for the twenty-first century of the Tarim River headwaters simulating glacier and hydrologic changes based on two glacio-hydrological models (SWIM-G and WASA). The use of an ensemble of global circulation model (GCM) results, as well as two hydrological models, is meant to expose the inherent uncertainty associated with both sources. Huang et al. (2018) use the projected discharge and consider the water management in the river oases downstream under the same climate change scenarios evaluated here.

Study region: the Tarim River headwaters
The Tarim River is one of the largest endorheic rivers in the world with a topographical catchment size of about 800,000 km 2 and a mainstream length of 600-800 km depending on discharge and water abstractions (Tao et al. 2011). The river is fed by three large tributaries with their confluence at the northern edge of the Taklamakan Desert: The Aksu River originating in the central Tien Shan to the north, the Hotan River originating in the Kunlun Shan to the south and the Yarkant River originating in the eastern Pamir and the Karakoram (see Fig. 1). The desert climate in the lower parts of these rivers and the Tarim produces virtually no river runoff (except for rare extreme rain events) with annual potential evapotranspiration exceeding precipitation by factors of 30-50 (Zhao et al. 2012). The vast majority of river discharge is generated in the mountainous and glacierised headwaters. This study focuses on the five catchments of the gauging stations that are situated at the boundary of the Taklamakan Desert, i.e. before significant river abstractions and transmission losses occur (Table 1 and Fig. 1).
The climate of the study region is highly continental with a strong seasonality governed by the Westerlies, with parts of the Hotan and Yarkant catchment also influenced by the Indian monsoon (Aizen et al. 1995;Maussion et al. 2014). As a consequence, river regimes exhibit a strong peak in summer when snow and glacier melt coincide with the precipitation peak. Mean summer (winter) temperatures fall between 3 and 20℃ (− 19 to − 5℃) with about 75% of precipitation falling between the months of April and September (in the 1971-2000 reference period).
The mountain ranges surrounding the Taklamakan Desert comprise steep, high-altitude terrain that has given rise to an extensive glacier cover accounting for significant proportions of the catchment areas considered here (Fig. 1, Table 1). The two Aksu headwaters comprise a total glacier area of 3410 ± 118 km 2 (2008, Pieczonka and Bolch 2015;Osmonov et al. 2013). A unique glaciological feature of the Aksu headwaters is the icedammed Merzbacher Lake that sends near-annually recurring subglacial outburst floods (jökulhlaups) downstream with consequences for communities and modelling efforts (Glazirin 2010;Wortmann et al. 2014). The two catchments of the Hotan River span the north-western edge of the Tibetan Plateau and have a glacier cover of some 5880 km 2 (Shangguan et al. 2007). The Yarkant headwater has a total glacier cover of about 5600 km 2 at extremely high and steep altitudes (Dyurgerov 2010).
Observed signs of a changing climate in the headwaters have been widely discussed. An increasing trend in river discharge over the past 50 years of the Aksu River, the most important tributary to the Tarim River, was reported to be as much as 30% between 1957 and 2004 relative to mean discharge (Wang et al. 2008;Krysanova et al. 2015). This was attributable to increasing trends in air temperature and precipitation, with a higher contribution of air temperature (leading to higher ice melt) in the Xiehela catchment, as found by a modelling and data-based study ). An increasing trend could, however, not be confirmed for the discharge of the Hotan and Yarkant catchments (Tao et al. 2011). The observed climate in Xinjiang has experienced a trend towards warmer and wetter conditions since the 1970s (Shi et al. 2006). Statistically significant increasing trends were found for temperature, precipitation and vapour pressure over nearly the entire Tarim catchment (at 1% significance level, Tao et al. 2011). Due to the poor observation density in the Hotan and Yarkant headwaters, however, those results must be interpreted with caution. Long-term, high-altitude observations for these catchments do not exist, rendering the hydrological observations the only longterm glimpse of these catchments' water balance and even those are interrupted. Glaciological studies in the Aksu catchment have received most attention, especially over the last decade. Pieczonka and Bolch (2015) have found a heterogeneous glacier mass and area loss over the catchment, but the average mass balance of 0.35 ± 0.34 m weq. a −1  is comparable to global values. Other studies have confirmed and contributed to this assessment Pieczonka et al. 2013;Osmonov et al. 2013). Investigations in the Karakoram based on declassified satellite imagery from the 1970s indicate stable conditions or only slight mass loss: Bolch et al. (2017) found a mass balance of − 0.01 ± 0.09 m weq. a −1  for the Hunza catchment and Zhou et al. (2017) found a mass balance of − 0.04 ± 0.05 m weq. a −1  for the central Karakoram. However, both studies focus on the part draining south into the Indus. Studies of mass changes since 2000 show a positive or stable balance in the Tarim part of the Karakoram, the Kunlun Mountains and eastern Pamir, e.g. + 0.04 ± 0.05 m weq. a −1 (2000, + 0.05 ± 0.07 m weq. Table 1 Catchment statistics for all considered stations with station elevation, annual and June-August mean discharge (Q, 1961-1989, with some gaps), glacier cover of 2008 Duethmann et al., 2016), mean annual precipitation (P, 1971-2000) of the APHRODITE dataset (Yatagai et al. 2012), corrected median values of precipitation Ps by Wortmann et al. (2018)

Data and methods
This climate change impact assessment considers future hydrological and glaciological changes by assessing multiple scenarios and climate model results by means of two glacio-hydrological models. For this study, the SWIM-G (Soil and Water Integrated Model -Glacier Dynamics; Wortmann et al. 2019) and WASA (Model of Water Availability in Semi-Arid Environments; Güntner 2002) models were as follows: (1) implemented in the five headwater catchments, (2) calibrated and validated to discharge and glacier observations, and then (3) driven by results of eight GCMs. What has become a common practice of climate impact assessment in many river catchments of the world (Fowler et al. 2007;Beniston 2003; Teutschbein and Seibert 2012) is a task marred with difficulties in the complex hydrology of the Tarim River; the stark contrast between desert and mountain climate, severe data scarcity, heterogenous past glacier changes , substantial river abstractions downstream (Tao et al. 2011) and regular large glacial lake outbursts (Wortmann et al. 2014) to name but a few. The validated models were then run with climate scenarios over the reference period and the twenty-first century. The reference period is defined as 1971-2000 and the near (2011-40), medium (2041-70) and far (2071-2100) future periods are investigated.

Input data
To overcome the precarious data availability, mostly preprocessed (homogenised, synthesised) datasets are used in this study (Table 2). They offer best-practice interpolations of Outlines from generated by Pieczonka and Bolch (2015) and Osmonov et al. (2013) and GlabTop simulated thickness Linsbauer et al., 2012) for model initialisation/calibration Discharge Daily river discharge from Chinese hydrological yearbooks for model calibration at stations: Xiehela (1964Xiehela ( -1987, Shaliguilanke , Wuluwati , 3 years missing), Tonggizuluoke , 5 years missing), Kaqun (1964Kaqun ( -1987  the sparse observations (homogenised also beyond national boundaries), while also circumventing restrictive data sharing policies of China. Precipitation is the most crucial driving variable and the renown APHRODITE dataset (Yatagai et al. 2012), a gridded precipitation interpolation of the densest gauge network in Asia, was used. Although it is the best data available, the station density is still extremely poor, e.g. the Hotan and Yarkant catchments are nearly devoid of any observations with the closest stations located at the edge of the Taklamakan Desert or at significantly lower-laying locations (Tao et al. 2011;Wortmann et al. 2018) leading to a negative bias. The necessary precipitation correction used is described in Duethmann et al. (2013, WASA) and Wortmann et al. (2018, SWIM-G).
Other driving data include temperature (daily mean, minimum, maximum), radiation and relative humidity, which are provided by the WATCH Forcing Data (v2) (Weedon et al. 2011). Although it suffers from the same or worse station sparsity as the APHRO-DITE data, the climate variables are significantly more stable over space and time . The variability of temperature with elevation (lapse rate) is also more stable-especially in dry climates-and is therefore parameterised in the snow and glacier melt components of the models.
Daily discharge observations for model calibration and validation were available for the five outlet gauges in the period 1961-1988 (with some missing years, see Table 2). Highresolution ice thickness simulations (GlabTob, Duethmann et al. 2015;Linsbauer et al. 2012) and regional glacier mass balances for the period ~ 1975-1999 were used for model initialisation and calibration with variations between the two glacio-hydrological models (see next section).

SWIM-G
The glacio-hydrological model SWIM-G is used to simulate both the glacier and catchment hydrology in a tightly integrated approach (Wortmann et al. 2019;. It was developed from the widely applied and tested semi-distributed ecohydrological model SWIM (Krysanova et al. 1998) that was recently extended by a glacier dynamics module and was specifically developed for long-term climate change impact assessments for medium to large river basins. Krysanova et al. (2015) provide an overview of the hydrological processes considered. The glacier dynamics extension includes all important glacier processes at the glacier scale, including ice flow, avalanching, sublimation and debris evolution (see Annex 1.1 for a brief description and Wortmann et al. 2019). It was developed to overcome the conflict of scale between individual glaciers and the catchment hydrology when modelling larger, highly glacierised, typically data-scarce catchments. It relies on computational units that disaggregate complex terrain into parts of similar elevation, aspect and hydrological subbasin as opposed to gridded (fully distributed) or empirical approaches (e.g. Immerzeel et al. 2011;Huss et al. 2010). They are conceptually similar and computationally equal to the hydrological response units (HRUs) implemented in SWIM, enabling a tight integration between glaciological and hydrological modelling over large catchments.
The model was calibrated and validated in the period 1961-1987 as in situ observations were available, using a multi-objective Pareto optimisation with three objective functions (calibration ensemble median ranges in parentheses): (1) to the observed discharge at the five outlet gauge stations using a combined Nash-Sutcliffe efficiency (0.73-0.94) and the root mean square error of mean annual discharge as objective functions (13-21%), (2) to observed glacier area using the hypsometry matching metric X 2 as objective function (8-15%) and (3) to mass balances using an error probability function (0.001-0.023). The calibration also includes the initialisation of the glacier cover by running the glacier dynamics driven by the climate of 1960-1975 repeatedly over 300 years to yield a glacier cover and thickness adapted to both the driving data and the model structure. The model initialisation and calibration is described in more detail in Wortmann et al. (2018) and the median performance values are given in Table A1 (Annex 1.2). For the Hotan and Yarkant catchments, long-term mass balance measurements were not available and the most plausible scenario was chosen instead. The model calibration yielded several 'non-dominated' parameter sets, all of which are considered equally valid. The parameter set that offers the best trade-off between all objective functions was chosen for the scenario runs presented in the results (see description and Figure A2 in Annex 1.3). For an uncertainty assessment, the scenarios were also run with the best/worst parameter sets and an analysis of variance (ANOVA, Gottschalk 2006) was conducted, comparing the variance between climate model ensemble, scenario and parameter set ( Figure A8 in Annex 2).

WASA
WASA is a semi-distributed, process-based hydrological model. It was first applied to semi-arid catchments in Brazil (Güntner 2002;Güntner and Bronstert 2004) and Spain (Francke et al. 2008;Mueller et al. 2010) and more recently to snow and glacier dominated high mountain catchments in Central Asia (Duethmann et al. 2013;2015). Detailed model descriptions including equations can be found in Güntner and Bronstert (2004) and Duethmann et al. (2015) with a brief description provided in Annex 1.4.
The present study uses a daily time step and the spatial structure is organized based on catchments, subcatchments and 200 m elevation zones. For calculating glacier area and volume changes, glaciers are considered individually and a finer discretization of 50 m elevation bands is used. The initial glacier ice thickness distribution of each glacier was derived from a spatially distributed ice-thickness model (GlabTop2, Linsbauer et al. 2012, Frey et al 2014. Based on the simulated mass balance of the glacier, thickness changes in each elevation band are calculated using predefined functions derived from observed glacier thickness changes of the Ak-Shirak massif in the period 1977-1999 (Surazakov and Aizen 2006). For the Hotan and Yarkant catchments, the model was calibrated to daily discharge data over a 10-year calibration period from 1979 to 1988. The period 1972-1978 was used for model evaluation. Additional 2-year periods were used for model initialization. Glacier mass balances for this region were not available for the period 1979-1988. It was assumed that glacier mass balances for the Hotan and Yarkant basins during the calibration period were close to balance, and the mass loss was constrained to 0 ± 0.1 m w.eq. a −1 Shean et al. 2020). For the Aksu catchment, the model was calibrated in a multi-objective way using daily discharge variations, interannual variations of seasonal discharge, discharge trends, correlation to the observed annual glacier mass balance series and cumulative glacier mass change based on geodetic estimates . The calibration period was defined from 1976-1999 and the periods 1957-1975 and 2000-2004 were used for model evaluation. From the set of Pareto optimal solutions, one solution was selected for further evaluations in this study. For details of the WASA model calibration in the Aksu catchment, please refer to Duethmann et al. (2015). Daily calibration (validation) Nash-Sutcliffe efficiency ranged between 0.60 and 0.85 (0.71-0.84). All performance results of the calibration are provided in Table A2 (Annex 1.4).

Climate change scenarios
We use three well-established scenarios from the IPCC Representative Concentration Pathways (RCP; IPCC 2014): (a) RCP 2.6 (atmospheric greenhouse gas concentration peaking around 2040 at 490 ppm CO 2 -eq and eventual decline), (b) RCP 4.5 (stabilisation at 650 ppm CO 2 -eq at the end of the twenty-first century) and (c) RCP 8.5 (continuous rise above 1370 ppm CO 2 -eq after 2100). The results of eight GCMs from the coupled model intercomparison project phase 5 (CMIP5) were used ( Table A3 in Annex 1.6, Taylor et al. 2012). The models were selected to cover the full range of precipitation and temperature change signals from the CMIP5 ensemble over the Tarim headwater catchments, i.e. the moderate and strong cases of wetter-warmer, wetter-colder, drier-warmer, drier-colder signals. The GCMs have spatial resolutions ranging from 1.5-3° and all data was provided at a daily temporal resolution.
GCM results were bias-adjusted relative to the baseline period 1971-2000 to the original driving data that the glacio-hydrological models were calibrated to , Duethmann et al. 2016. A non-parametric quantile mapping approach with trend preservation was chosen, as was previously used by Hempel et al. (2013). Despite concerns of the validity of using a bias-adjustment for climate change impact assessments (Ehret et al. 2012), it was considered necessary in the Tarim headwaters because of considerable deviations in precipitation between calibration and scenario driving data as well as the great sensitivity of the glacier cover to even slight differences in climatic conditions. The application was paramount to ensure plausible reference conditions for the ensemble assessment.

Changes in temperature and precipitation
Increases in both temperature and precipitation are projected under most climate scenarios, with significant changes in all headwater catchments relative to the reference period 1971-2000 (Fig. 2 for annual and period changes, Figures A4-A5 in the Annex for monthly regimes of all periods and catchments). All ensemble median values in the three projection periods increase across all regions. Only ensemble minimum signals indicate possible negative changes in precipitation, especially in the Aksu. Strong increases in (median) temperature and precipitation from the reference period to the near future are projected (often stronger than changes in subsequent periods).
Temperature increases are similar across all regions. In the near future, they range between 0.5 and 2.4℃ with only small differences between the scenarios. Those differences become more striking in the two later periods, when ensemble median changes increase with RCP scenario. In the medium-term future, they peak for the RCP 2.6 scenario at about 2.5℃ and slightly decrease thereafter, while for the higher RCP scenarios, they continue to increase up to 7℃ (ensemble max.) in the far future for RCP 8.5. The ensemble variability generally increases with time, leading to ranges of about 2-4℃ in the far future. As expected, snowfall is more confined to the winter months, as time and emission scenario progresses ( Figure A6 for monthly snow fractions in the Annex).
Changes in precipitation are less uniform across regions and projection periods and exhibit an even greater ensemble variability, a pattern common in other regions (Thompson et al. 2013;Vetter et al. 2013). Relative changes in the near future are modest and similar across the regions and scenarios with median values of about 6-8%. Scenario differences become more pronounced in the two later periods, as ensemble median changes vary with RCP scenario (with the exception of the medium-term in the Fig. 2 Climate change scenarios for the Tarim headwaters for the twenty-first century, including the reference period . Ensemble maximum, median and minimum values are shown as 10-year running means and signals averaged over the near (2011-2040), medium (2041-2070) and far (2071-2100) future. Absolute values are given in the left vertical axis and changes relative to the reference period average along the right axis. Note that these are bias-adjusted GCM results without the static precipitation bias correction applied by the glacio-hydrological models. Monthly regimes for all periods, scenarios and catchments are provided in Figure A4-A5 in the Annex Aksu catchment). Changes in the Aksu catchment range from decreases of up to 25% to increases of the same magnitude, but the majority of models indicate increases in precipitation with median changes of 11-18%. These increases are similar in the Hotan and Yarkant catchments, but the spread is large and mostly positive (− 15 to + 54%).

Changes in glacier area and volume
A receding trend over the twenty-first century is evident in all catchments for the ensemble medians of both models and it is strengthening with RCP scenario (Fig. 3, period mean changes in Table A4 in the Annex). For the ensemble medians, glaciers are projected to  Table A4 in the Annex shrink by up to 35% in the near future, 2-64% in the medium future and 8-89% in the far future compared to the reference period across the three scenarios. and large ranges are due to differences between two impact models and RCP scenarios. In the Xiehela and Shaliguilanke (Aksu) catchments, area shrinkage for the high-end RCP 8.5 scenario steadily rises to 55-80% and 80-95% by the end of the century, respectively, compared to 2010. Both lower scenarios show lower levels of 32-49% and 51-74%, respectively (Fig. 3).
For the SWIM-G projections, median shrinkages are the highest for Shaliguilanke (Aksu) and lowest for Tongguziluoke (Hotan) catchment. Glacier changes in this highly glacierised catchment are similar to those in the Xiehela (Aksu) catchment (both with a catchment glacier cover ≈ 20%), i.e. reaching about − 50% until 2100 under RCP8.5. The long-term changes of the less glacierised Wuluwati (Hotan) and Kaqun (Yarkant) catchments are similar to the Shaliguilanke (Aksu) catchment (7, 12 and 5% glacier cover, respectively), i.e. close to − 75% under RCP8.5. Differences between the lower two scenarios, however, are larger in the Hotan and Yarkant catchments than in the Aksu, and they show a slower recession in the first half of the century.
The median shrinkages of the WASA projections are also the highest for Shaliguilanke and the lowest for Tongguziluoke. Both Wuluwati and Kaqun show lower levels than Xiehela (all periods) but higher than Tongguziluoke (from the mid-century). In the near future period, practically no changes are projected for the Yarkant and Hotan. The changes for the end of the century are also the highest for Shaliguilanke (reaching almost − 100% under RCP8.5), and the lowest for Tongguziluoke (> − 25%). The final changes by the end of the century are similar for Xiehela and Kaqun (reaching approx. − 75% at RCP8.5), and moderate at Wuluwati.
Comparing projections between the two impact models, it is evident that WASA simulates a much stronger decline of glacier areas in the Aksu basin than SWIM-G from the middle of the century. However, SWIM-G generates higher losses than WASA in the Hotan and Yarkant catchments, though by the end of the century both models come to similar results for the Kaqun catchment (Yarkant).
In line with area changes, glacier volume (in water equivalents) is also projected to decrease under all scenarios (Table A4, Figure A6 in the Annex). Projected mass loss is 4-25% in the near, 12-75% in the medium and 18-97% in the far future across all catchments and scenarios for the ensemble median. Losses are consistently the greatest in the Shaliguilanke (Aksu) catchment, lower in the Xiehela (Aksu), Wuluwati (Hotan) and Kaqun (Yarkant) catchments, and they are the lowest in the Tongguziluoke (Hotan) catchment.
Glacier mass balances give an indication of the glacier imbalance, the hydrological impact thereof and they may be better compared to studies of past glacier evolution. Decadal mean annual mass balances are provided in Figure A7 (Annex). The glacier recession described above is associated with negative mass balances that maintain or exceed the negative rates of the past in the first half of the twenty-first century with a recovery for the RCP 2.6 and 4.5 scenarios and an acceleration for the RCP 8.5 scenario by the end of the century. The near future  shows similar mass balances to those observed in the past across the region -about − 0.4 to 0 m weq. a −1 -with little differences between scenarios. The largest negative mass balances for the two lower-end scenarios are projected for the middle of the century, while they continue to grow more negative in the high-end scenario until they reach − 1.2 to − 0.8 m weq. a −1 by the end of the century. A recovery to stable or even positive mass balances is only projected under the RCP 2.6 and 4.5 scenarios in the far future analogues to the temperature projections (Fig. 2).

Changes in river discharge
The discharge projections of the multi-model analysis are summarised in Fig. 4 for each catchment and period. In the two Aksu subcatchments, the increase in mean annual discharge is most pronounced in the near future but recedes in subsequent periods except for the RCP 8.5 scenario simulated by SWIM-G in the Xiehela catchment. The initial increase is 10 to 25% for the ensemble median relative to the 1971-2000 reference period (higher values from WASA compared to SWIM-G) with a slight increase relative to the RCP forcing. In the medium and long-term, the increase in the Aksu catchment reduces progressively and turns to no or negative changes compared to the reference period, especially under the RCP 2.6 scenario in the Xiehela catchment and under the RCP 8.5 scenario in the Shaliguilanke catchment.
The impacts are different in the Hotan and Yarkant catchments. There is a general trend of increasing mean annual discharge apparent at all stations and nearly all time periods. Changes in mean annual discharge in these catchments show a uniformly increasing trend at greater magnitudes as in the Aksu catchment. The initial increase in the near future is around 25 to 40% (ensemble median). Subsequent periods also show increases, especially under RCP 8.5. The RCP 4.5, simulations peak in the medium-term at 28-60% compared to the baseline period, while the RCP 2.6 scenario shows lower or similar increases in the later periods than in the first.
The changes in discharge are predominantly concentrated in the summer months (Fig. 5). April to June discharge predominantly increases in the Aksu catchments for all scenarios (mostly higher increase according to WASA), July-September discharge tends to increase in the near future but decreases in the far future. The Hotan and Yarkant catchments show increases in the summer months (June to September) that progress with period and RCP scenario. Discharge only decreases under the RCP 2.6 and RCP 4.5 scenarios in Fig. 5 Absolute regime changes in the three main headwaters of the Tarim River (sum of both headwaters in the Aksu and Hotan). The ensemble median is shown for each scenario and glacio-hydrological model. Rather than relative changes for each month in comparison with the reference period, the absolute monthly changes are shown to account for the highly seasonal flows August-September in the Yarkant in the far future. Winter discharge changes in all catchments are mostly positive but represent only small fractions of the increase in the annual discharge. Changes in discharge according to SWIM-G are mostly lower than WASA in the Aksu catchment but higher in the Hotan and Yarkant.
The origin of these changes can be traced by investigating the changes in the runoff generating water inputs, i.e. rain, glacier and snow melt (Fig. 6). As would be expected under warmer and wetter climate conditions, the rain is increasing in all catchments under most scenarios and the glacier melt is mostly decreasing. An exception is the highly glacierised Tongguziluoke (Hotan) catchment, where both rain and glacier melt are increasing, indicating an increased redistribution of ice into the ablation zone and explaining the sharp projected increases in discharge. The snowmelt is nearly constant or slightly increasing, except in the Shaliguilanke (Aksu) catchment according to WASA under the RCP 8.5 scenario. The changes in discharge components reveal the peak in glacier melt predominantly in the first half of the twenty-first century for the Hotan and Yarkant catchments (Fig. 6). Under the RCP 8.5 scenario, the peak is shifted further into the second half of the century, which is especially evident in the Tongguziluoke (Hotan) catchment.

Discussion
Results indicate a marked increase in river discharge, which is in-line with observations over the past decades (Tao et al. 2011;Li et al. 2020). Given the arid climate combined with the vast ice reserves, higher temperatures leading to negative glacier mass balances are able to drastically change river discharge in the short-and medium-term until 'peak meltwater' is reached Huss and Hock 2018). The two existing studies estimate the peak of the glacier meltwater to be reached in about 2030 in the low emission scenario and about 2060 or 2070 in the high emission scenario for the whole Tarim basin (Huss and Hock 2018;Rounce et al. 2020). Our study showed that these projected peaks for the whole basin can mainly be attributed to the Asku basin as Hotan and Yarkand show a distinctly different behaviour. The peak water seems to be in sight under all scenarios in the Aksu catchments as both glacio-hydrological models indicate the largest increases in the near future (Fig. 4, 5). In the Hotan and Yarkant catchments, this peak is only evident in the medium to far future mainly under the low and moderate emissions scenarios. Projected significant increases in precipitation are also driving the strong increases in discharge and are able to compensate for melt water losses. However, precipitation is quite variable from year to year while glacier melt is a more reliable water source Pritchard 2019). Moreover, precipitation projections are subject to high uncertainties.
The use of two glacio-hydrological models with different representations and processes of glaciers allows a more robust understanding of model uncertainties. Both models generally agree on the direction of glacial and discharge changes with the exception of discharge in the Aksu catchment in the far future, where ensemble ranges point in both directions but the medians disagree. Disagreements in change magnitude are most evident in the glacier recession in the Aksu (WASA greater reductions) and the Hotan catchments (SWIM-G greater reduction) with WASA indicating a stronger increase in discharge in the short-and medium-term in the Xiehela catchment and SWIM-G predicting slightly stronger median increases in the Hotan catchments. Ensemble ranges of discharge changes (Fig. 4) are generally greater for WASA with some exceptions in the Tongguziluoke catchment.
While generally in agreement with previous studies (Fang et al. 2018;Liu et al. 2011;Wang et al. 2021;Pieczonka and Bolch 2015), water balance components also show some strong differences between the models, as WASA generally shows higher values of rain and snow melt and SWIM-G larger values of glacier melt in the Hotan and Yarkant catchments. The differences between the glacio-hydrological models is predominantly caused by different assumptions in the precipitation correction, a key calibration term and uncertainty in these catchments (cf. P in Table 1), as well as by differences in the calculation of evapotranspiration. SWIM-G simulates lower rates of rain and snowmelt and in the Hotan and Yarkant catchments in the first half of the twenty-first century also higher rates of glacier melt, since WASA assumes a stronger precipitation correction. While the uncertainties between the two models with respect to precipitation correction appear large, these are actually small when compared to the precipitation estimates based on observation products or to simulated precipitation by GCMs. Thus, glacio-hydrological modelling in combination with observed discharge and glacier mass balances or glacier area distribution can be a suitable alternative to derive precipitation estimates. An ensemble of glacio-hydrological models is useful to assess the uncertainties of such precipitation estimates.
The differences in flow components may explain some of the differences in the changes in annual discharge. Higher contribution of rain makes WASA more susceptible to the higher ensemble variability in precipitation leading to greater uncertainty ranges. Higher shares of glacier melt in SWIM-G in the Tongguziluoke catchment lead to greater sensitivity to temperature increases, producing greater increases in discharge. Differences in glacial changes are likely attributable to the different representations of ice in the models. For example, the glacier recession in WASA is steeper in the first half of the twenty-first century in the Aksu catchments leading to a pronounced glacier melt peak, while SWIM-G produces more gradual changes. This may be caused by the different approaches for calculating glacier area and volume changes, and processes that are considered by the SWIM-G model but not by the WASA model, such as glacier dynamics, sublimation and debris accumulation (Wortmann et al. 2019).
Model parameter uncertainty was quantified using the calibration and climate change scenario ensemble via an analysis of variance (for SWIM-G in Figure A5 in the Annex and for WASA in Duethmann et al. (2016) Fig. 4). For discharge, the variance related to different parameter sets amounts to less than 20% of the total, while the differences related to the climate model ensemble makes up more than 50%. Reductions in the different sources of uncertainties will ultimately only come from better observations to further constrain parameter ranges and improve modelled process representations . For example, high elevation precipitation measurements could help constrain the precipitation correction and help improve estimates of elevation-dependent warming (MRIEDWWG 2015) or stable isotope measurements could reveal the share of runoff from rain, snow and glacier melt (Kumar et al. 2011;Ohlanders et al. 2013;He et al. 2019). Lastly, improvements in the GCM projections would likely have the largest impact on uncertainties in the glacier and discharge projections presented here, as the use of a GCM ensemble has shown. This study may serve as a justification for future investments in the improvement of these broader deficiencies.

Conclusions
This study represents the first systematic and comprehensive climate change impact assessment of the Tarim River headwater catchments for the twenty-first century, using state-of-the-art GCM climate projections and scenarios as well as two large-scale, glacio-hydrological models. Together, changes in these catchments will largely determine the water resources of the Tarim River and the communities living in and around the Taklamakan Desert. The assessment is the product of significant research efforts to simulate the complex water resources of these remote and largely ungauged mountain ranges. It required the construction of two large-scale glacio-hydrological models (Wortmann et al. 2019;Duethmann et al. 2016) and the correction of precipitation datasets as well as the subsequent model calibration under data-scarce conditions Duethmann et al. 2013). These efforts have enabled a robust impact assessment based not only on a single impact model but two that can expose the modelinherent uncertainties.
In line with regional and global trends, results indicate a warmer and wetter climate with substantial consequences for the high-mountain glacier cover and the meltwaterdriven rivers. Depending on the low, medium or high emission scenarios, temperatures are projected to rise by about 1.9 °C, 3.2 °C or 5.3 °C in the ensemble mean by the end of this century (2071-2100) compared to the 1971-2000 reference period, while precipitation may intensify by about 9%, 14% or 24%, respectively. The two calibrated and specifically adapted glacio-hydrological models, SWIM-G and WASA, project these climatic changes to result in a shrinkage of the glacier cover by around 45%, 52% or 73% (median of both models and catchments) in the Aksu headwaters and by about 21%, 28% or 40% in the Hotan and Yarkant headwaters for the low, medium or high emissions scenarios, respectively. Similarly, the models project river discharge to mainly increase in the Hotan and Yarkant catchments by the end of the century (by about 15%, 30% or 60%), while in the Aksu headwater, it is projected to first increase (by about 20% regardless of scenario) to return to or even decrease below reference discharge by the end of the century. Discharge increases are the greatest in the spring and summer months (April-September) and decreases in the Aksu catchment mainly occur in August-September due the diminished glacier meltwater.
Despite the shrinking glacier cover, precipitation increases drive the projected changes in river discharge and at least partially compensate the losses in meltwater. The receding trend in the Aksu headwaters is important, however, as they currently dominate the discharge of the Tarim River. With the projected increases in the Hotan and Yarkant headwaters, this could change, depending on the agricultural water abstractions in the river oases lining these tributaries (Huang et al. 2018). Future research should focus on reducing the large uncertainties that remain in the simulation of glacio-hydrological changes in this region that are mainly driven by the climate ensemble, the glacio-hydrological models and the limited observations available for constraining the models during the reference period.
Author contribution M.W., D.D. and V.K. designed the study. M.W. performed the SWIM-G simulations and drafted the manuscript. D.D. carried out the WASA simulations and contributed to the draft. C.M. processed and analysed the climate input data. T.B. provided the glacier input and validation data. S.H., J.T. and Z.K. were involved in data acquisition and research design. All authors contributed to the final version of the manuscript.
Funding This study was conducted within the project SuMaRiO (Sustainable Management of River Oases along the Tarim River; http:// www. sumar io. de/), funded by the German Federal Ministry of Education and Research (BMBF grants 01LL0918J, 01LL0918I and 01LL0918B). T. Bolch acknowledges funding by Deutsche Forschungsgemeinschaft (DFG, BO3199/2-1). Open Access funding enabled and organized by Projekt DEAL.