Central Himalayan tree-ring isotopes reveal increasing regional heterogeneity and enhancement in ice mass loss since the 1960s

Tree-ring δ¹⁸O values are a sensitive proxy for regional physical climate, while their δ¹³C values are a strong predictor of local ecohydrology. Utilizing available ice-core and tree-ring δ¹⁸O records from the central Himalaya (CH), we found an increase in east–west climate heterogeneity since the 1960s. Further, δ¹³C records from transitional western glaciated valleys provide a robust basis for reconstructing about 3 centuries of glacier mass balance (GMB) dynamics. We reconstructed annually resolved GMB since 1743 CE based on regionally dominant tree species of diverse plant functional types. Three major phases became apparent: positive GMB up to the mid-19th century, the middle phase (1870–1960) of slightly negative but stable GMB, and an exponential ice mass loss since the 1960s. Reasons for accelerated mass loss are largely attributed to anthropogenic climate change, including concurrent alterations in atmospheric circulations (weakening of the westerlies and the Arabian Sea branch of the Indian summer monsoon). Multi-decadal isotopic and climate coherency analyses specify an eastward declining influence of the westerlies in the monsoon-dominated CH region. Besides, our study provides a long-term context for recent GMB variability, which is essential for its reliable projection and attribution.

1 Introduction

Glaciers in the Himalayan–Tibetan orogen are an important component of the regional hydrological cycle, and a major fraction of regional potable water is stored and provided by them. However, recent climate warming has imposed a serious alteration on the equilibrium of these glaciers (Bandyopadhyay et al., 2019; Bolch et al., 2012; Maurer et al., 2019; Mölg et al., 2014; Yao et al., 2012; Zemp et al., 2019). High uncertainty prevails in future projections of glacier mass balance (GMB), since sound understanding of glacier fluctuations and GMB response to climate change on multi-decadal timescales is poorly understood (e.g. the erroneous statement in the Fourth Assessment Report of the IPCC; Kargel et al., 2011). Reliable projections of future Himalayan ice mass loss require robust observations of glacier response to past and ongoing climate change. Long-term estimation of GMB is also imperative for regional water security. Currently, coupled glacier–climate models do not even agree on the sign of change, and hence projections of GMB are ambiguous (Watanabe et al., 2019; Jury et al., 2019; DCCC, 2018). Nevertheless, a consistent picture emerges of net ice mass loss in recent decades, which is highest in the western and central Himalaya (except the Karakoram and the Pamir Mountains) (Bolch et al., 2012; Brun et al., 2017; Dehecq et al., 2019; Maurer et al., 2019; Mölg et al., 2014; Shekhar et al., 2017; Yao et al., 2012).

The central Himalayan glaciers show a rather homogeneous behaviour (Azam et al., 2018; Bandyopadhyay et al., 2019; Brun et al., 2017; Dehecq et al., 2019; Kääb et al., 2012; Sakai and Fujita, 2017; Yao et al., 2012). In this study, we focus on the transitional climate zone between the western and central Himalaya, where knowledge about multi-decadal glacier dynamics in relation to climate change is still missing. Evidence from tree-ring isotopes and hydroclimatic studies from the central Himalaya suggest that the glacier mass balance behaviour is primarily determined by the conjoint effect of the winter westerlies (WD) and Indian summer monsoon (ISM) (Fig. 1 and references therein). The influence of the ISM declines towards the northwest Himalaya, and the westerlies progressively become dominant. Towards the eastern Himalaya, the influence of the Arabian Sea branch of the ISM declines, and the Bay of Bengal branch dominantly regulates the climate, besides an influence of the East Asian monsoon (Benn and Owen, 1998; Bookhagen and Burbank, 2010; Hochreuther et al., 2016; Liu et al., 2014; Lyu et al., 2019; Mölg et al., 2014; Sano et al., 2017; Yang et al., 2008; Yao et al., 2012). Monsoon-influenced glaciers, particularly those in the transitional climatic zone (such as the western region of the central Himalaya), are more sensitive to climate warming than winter-accumulation type glaciers (Kargel et al., 2011; Sakai and Fujita, 2017). Moreover, a weakening in the moisture-delivering systems (i.e. WD and ISM) since the mid-20th century have had a direct impact on the summer-accumulation type glaciers in the region (Qin et al., 2000; Hunt et al., 2019; Joswiak et al., 2013; Khan et al., 2019; Roxy et al., 2015, 2017; Sano et al., 2012, 2013, 2017; Singh et al., 2019; Xu et al., 2018; Yadav, 2011).

Figure 1 (a) The study region at the transitional climate zone of the western central Himalaya (WCH) showing tree-ring sampling sites, meteorological and aethalometer stations, and four benchmark glaciers of the Uttarakhand Himalaya (b: Dokriani Glacier; c: Chorabari Glacier; d: Tipra Bamak; e: Dunagiri Glacier). Tree-ring sites from the WCH include (1) Manali (Sano et al., 2017), (2) Uttarakashi (Singh et al., 2019), and (3) Jageshwar (Xu et al., 2018). The Dasuopu ice core site in the eastern part of the central Himalaya (ECH) has been indicated with a blue circle with a dot (Thompson et al., 2000), while tree-ring sites of the ECH are (5) Ganesh (Xu et al., 2018) and (6) Bhutan (Sano et al., 2013). The fourth tree-ring site (Humla; Sano et al., 2012) has been indicated between the WCH and ECH.

In the present study, we attempt to reconstruct glacier mass balance history since the end of the Little Ice Age (since 1743 CE) in the transition zone between the ISM and the westerly-dominated climate in the western central Himalaya (Fig. 1). We use tree-ring stable-isotope (δ¹³C and δ¹⁸O) chronologies of three dominant tree species of two different plant functional types and synthesize available regional δ¹⁸O chronologies from different archives, including tree rings and ice cores across the central Himalaya (Fig. 1; Table S1).

Field-based mass balance measurements are logistically challenging in the Himalaya. Nevertheless, since the 1980s, workers have succeeded in monitoring four benchmark glaciers (Pratap et al., 2016; Garg et al., 2017, 2019) of the Uttarakhand Himalaya (Fig. 1). Based on a detailed analysis of hypsometric curves and several morphological and glaciological factors of these four valley glaciers (Dokriani Glacier, DOK; Chorabari Glacier, CHO; Tipra Bamak, TIP; Dunagiri Glacier, DUN; Fig. 1, Table S2) (Bandyopadhyay et al., 2019; Pratap et al., 2016; Garg et al., 2017, 2019) and building upon previous studies (Shekhar et al., 2017; Bandyopadhyay et al., 2019; Azam et al., 2018), we compiled a mean observed glacier mass balance time series starting in the 1980s (1982–2013) (Figs. 1, S1; Table S2).

The length of observed mass balance records principally limits our understanding of the response of Himalayan glaciers to climate change. Remote sensing techniques provide the best alternative in this regard (Garg et al., 2017, 2019; Maurer et al., 2019; Kääb et al., 2012; Dehecq et al., 2019; Brun et al., 2017; Bandyopadhyay et al., 2019; DCCC, 2018). Nevertheless, observation-based mass balance records cannot be extended beyond a few decades. To understand long-term changes, we must rely on proxy data. One valuable proxy is ice-core isotopes, which are yet to be obtained from the Indian Himalaya. Tree rings are another reliable and sensitive proxy allowing the reconstruction of glacier history. They have been widely used in mountainous regions around the world (Bräuning, 2006; Duan et al., 2013; Gou et al., 2006; Hochreuther et al., 2015; Larocque and Smith, 2005; Linderholm et al., 2007; Nicolussi and Patzelt, 1996; Solomina et al., 2016; Tomkins et al., 2008; Watson and Luckman, 2004; Xu et al., 2012; Zhang et al., 2019). Tree-ring width (TRW) is known to bear the influence of non-climatic factors like biological tree age, which overlap climatic signals (Bunn et al., 2018; Zhang et al., 2019). TRW records of Himalayan evergreen conifer trees (including Cedrus deodara, Abies, Larix, and Picea spp.) have been widely used to reconstruct hydroclimate and natural hazards (Borgaonkar et al., 2009, 2018; Cook et al., 2003; Singh et al., 2006; Singh and Yadav, 2013; Yadav and Bhutiyani, 2013; Yadav et al., 2011; Yadav, 2011). Despite numerous TRW-based climate reconstructions, there is only one TRW-derived glacier mass balance reconstruction (Shekhar et al., 2017). Tree-ring stable isotopes are a reliable source of past hydroclimate variability because of their sensitivity to local and regional climate, the coherence in their climatic response, and the well-understood control by the environmental conditions regulating tree physiology (Levesque et al., 2019; Sano et al., 2012, 2013, 2017; Singh et al., 2019; Zeng et al., 2017). Although many studies have reconstructed hydroclimate in the Himalaya using tree-ring isotopes (e.g. Sano et al., 2012, 2013, 2017; Singh et al., 2019; Xu et al., 2018), glacier mass balance reconstructions from tree-ring isotope series are still not available. Higher growing-season temperatures stimulate photosynthesis rates. This facilitates a stable-isotope fractionation process, resulting in decreased intercellular CO₂ concentration and fractionation against the heavier ¹³C isotope, which leads to increased δ¹³C. Thus, climate (temperature) is a bridge that indirectly connects the mass balance of a glacier with tree-ring C isotope ratios. Hence, the mass balance can be reconstructed using stable-carbon-isotope chronologies from trees growing in the proximity of the glacier (Zhang et al., 2019).

In this study, therefore, we utilized tree-ring carbon isotope (δ¹³C) chronologies of three different species belonging to two plant functional types (PFTs) that differ in their annual phenological cycle (Table S1). One PFT includes evergreen conifers (Abies pindrow and Picea smithiana) that have been widely used to reconstruct hydroclimatic regimes across the Himalaya and are abundantly found in the moist valley regions of the high central Himalaya. Moreover, we investigated the dendroclimatological potential of a major dominant species of broadleaf deciduous trees (Aesculus indica), as its growth period (April to September) coincides with the warm–wet phase during the Indian summer monsoon (ISM) season (and lacks westerlies' signal due to complete winter dormancy). Additionally, we utilized available Dasuopu ice core and tree-ring isotope chronologies from six different locations in the central Himalaya to substantiate our results (Fig. 1). This study aims (i) to identify the tree-ring parameter (TRW and/or δ¹³C, δ¹⁸O) of deciduous and evergreen tree species containing the strongest mass balance signature; (ii) to reconstruct annual glacier mass balance using the identified best proxy chronology; and (iii) to analyse the reconstructed mass balance in relation to regional climate, local forcing factors, and large-scale atmospheric circulation.

2 Materials and methods

2.1 Study region and climate

The present study focuses on the glaciers of the Uttarakhand Himalaya (Indian central Himalaya) in the transitional climate zone between the western and the central Himalaya. The study region extends from latitude 30.15 to 31.03^∘ N and longitude 78.78 to 80.73^∘ E. The highly glaciated upper reaches of the Uttarakhand Himalaya encompass about 1000 glaciers of varying size, which cover a total area about 2850 km² (Raina and Srivastava, 2008; Bandyopadhyay et al., 2019). Previous studies have found that glaciers in the region are mainly fed by the ISM (Azam et al., 2018; Sakai and Fujita, 2017). Therefore, glaciers in the region are usually classified as summer-accumulation type glaciers (Azam et al., 2018; Sakai and Fujita, 2017). Snow accumulation during winter (December–March) is influenced by a precipitation regime driven by mid-latitude westerlies (WD). Decadal-scale meteorological records from the study region are available only for lower elevations (< 2000 m a.s.l.) (Fig. 1). Therefore, long-term gridded temperature and precipitation datasets were obtained through the Climatic Research Unit (CRU TS v. 3.22, 0.5 latitude × 0.5 longitude, 1901–2015) (Harris et al., 2014). However, given the limitations of the CRU precipitation dataset in high-altitude regions, we complemented it with existing meteorological records (Singh et al., 2019). Analyses of meteorological records indicate that mean annual precipitation is ∼ 800 mm, of which the warm–wet summer months (April–September) receive about 80 %. Mean annual temperature varies around 4 ^∘C, with a minimum (−7 ^∘C) in January and maximum (13.5 ^∘C) in June (Fig. S2). As previously mentioned, studies indicate an increasing trend in temperature, while precipitation during both ISM and WD seasons show a recent declining trend. As local climate is the common factor influencing both glacier mass balance and the physiology and growth of the trees growing in the valleys, we first analysed tree growth–climate relations (Fig. S3). Averaged monthly temperature and total monthly precipitation over a period spanning October of the preceding year to September of the current growth year were correlated with δ¹³C chronologies of the species (as δ¹³C is a strong predictor of local climatology) (Fig. S3). In addition, various climatological indices (sea surface temperature, SST; the El Niño–Southern Oscillation, ENSO; ENSO Modoki; the Pacific Decadal Oscillation, PDO; and the Indian Ocean Dipole, IOD) were obtained from the archive (http://climexp.knmi.nl, last access: 23 January 2020) to analyse large-scale climatic relations.

2.2 Glacier mass balance data

Four benchmark valley glaciers that are distributed across the Uttarakhand Himalaya, namely Dokriani (DOK), Chorabari (CHO), Tipra Bamak (TIP) and Dunagiri (DUN) glaciers, have been individually monitored for their mass balance over the time period from 1982 to 2013 (with a few gap years; Fig. 1, Table S2) (Pratap et al., 2016; Garg et al., 2017, 2019). Two glaciers (DOK and CHO) lie in the Garhwal Himalaya, while two (DUN and TIP) are located in the Kumaun Himalaya (India). We used the available glaciological mass balance records to produce the best possible time series (Table S2). We built a mean mass balance time series based upon previous studies (Shekhar et al., 2017; Bandyopadhyay et al., 2019; Borgoankar et al., 2009) and the premises that the glaciers show homogeneous behaviour in response to changes in the regional climate and that their mass balances show a high inter-correlation (Azam et al., 2018). Given the same climatic forcing owing to their location in the humid central Himalaya and their similar geomorphological characteristics (Garg et al., 2017, 2019) (Table S2), the four valley glaciers should generally correlate well (Solomina et al., 2016). Ideally, land-terminated valley type glaciers of a simple configuration and hypsometry are most suitable for inferring palaeoclimatic information (Solomina et al., 2016). The four studied glaciers are of this kind and are therefore designated by workers as “benchmark glaciers”. The glaciers are of a similar small size (except DUN), have simple shapes, and have moderate elevation ranges, and thus their recent deglaciation responded to changes in the climate equivalently (Table S2). A surface ice velocity study (Table S2) confirmed a similar response time of these glaciers (Garg et al., 2019). Geomorphological features determining glacier hypsometry primarily regulate the response of an individual glacier to changes in the climate. Thus, we prepared hypsometric curves of these four benchmark glaciers (Fig. S1), which were convex for the DOK and CHO glaciers and concave for the TIP and DUN glaciers. These curves and other similar glacio-morphological indices (Garg et al., 2017, 2019) (Table S2) provide a reasonable basis to compile these individual mass balance time series (Shekhar et al., 2017; Bandyopadhyay et al., 2019).

2.3 Tree-ring data

The three tree species utilized in this study are ubiquitously distributed throughout the western central Himalaya and have been used earlier to reconstruct past climatic variations. The collected species occupy both aspects of a glacier valley, stretching from near the treeline to downslope towards the valley. About 20 to 30 increment core samples of healthy tree individuals (to minimize the influence of non-climatic factors on growth) were collected from higher-elevation sites (2500–3800 m a.s.l.) from two representative glacial valleys (Dokriani – DOK, Pindar – PIN) encompassing the basin in the Indian central Himalaya (Fig. 1). Tree-ring isotopic signals (particularly δ¹⁸O) are coherent over a large area, as shown by high inter-species and inter-site correlations in the Himalaya (Sano et al., 2012, 2013, 2017; Xu et al., 2018; Grießinger et al., 2019; Table S1). Recently, Singh et al. (2019) reconstructed regional ISM precipitation derived from δ¹⁸O chronologies of the same three tree species from the DOK valley and found high inter-species and inter-site correlations. Due to the unavailability of long-term tree-ring stable-isotope records from PIN valley, we utilized δ¹³C chronologies of the three tree species from the DOK valley.

In summary, we adopted a methodology to reconstruct glacier mass balance history from tree rings that had been successfully tested before (Nicolussi and Patzelt, 1996; Duan et al., 2013; Zhang et al., 2019). We used standard dendrochronological methods and techniques to develop tree-ring-width (TRW) and stable-isotope chronologies. Then, we calibrated the climate response using linear-regression models and tested the reliability of our reconstruction. The leave-one-out cross-validation method (LOOCV; Michaelsen, 1987; Yadav and Bhutiyani, 2013; Duan et al., 2013; Zhang et al., 2019) was used to verify our reconstruction, given the relative shortness of glacier mass balance data (23 years, 1982–2013 after omitting gap years of 1990–1992, 1996, 1997, 2001–2003, and 2011) (Table S3).

2.3.1 Tree-ring-width chronology development

Two core samples per tree were collected at breast height using 5.15 mm diameter increment corers. Standard dendrochronological procedures (Fritts, 1976; Holmes, 1983) such as mounting and surface smoothing were applied to render the ring boundaries clearly visible. TRW of the samples was measured at a resolution of 0.001 mm using a LINTAB^™ system interfaced with a computer. Cross-dating was performed by matching variations in ring widths in all cores to determine the absolute age of each ring. Dating and ring-width-measurement quality control was conducted using the COFECHA computer program (Holmes, 1983). A coherent growth pattern between species revealed a common regional climate signal affecting growth of the trees. However, TRW chronologies of all species only showed a weak correlation with glacier mass balance data and hence were not used for any reconstruction efforts.

2.3.2 Stable-carbon-isotope chronology development

Five trees per species were selected for stable-isotope analyses based on best TRW inter-series matches. Each year's growth rings were dissected with precaution with a sharp scalpel under the microscope. To remove any possible juvenile effect, the innermost approximately 30–40 rings of each tree core were omitted from the analyses. Wood samples were grounded using an ultracentrifuge mill (Retsch ZM1). Extraction of cellulose from whole wood and carbon isotope analysis was carried out at the Institute of Geography, University of Erlangen–Nuremberg, Germany. Cellulose was extracted from the wood samples using the method of Wieloch et al. (2011). Isolated cellulose was homogenized using an ultrasonic method (Laumer et al., 2009) and freeze-dried. Before the pooling procedure, we checked co-relatedness in all five individual time series at 20-year intervals across the entire chronology. About 270 µg of cellulose was weighed into tin capsules with a microbalance (ME36S, Sartorius, Germany). The carbon isotope analyses were performed with an elemental analyser (NC 2500, Carlo Erba, Italy) linked to an isotope ratio mass spectrometer (IRMS; DELTAplus, Thermo Finnigan, Germany). Prior to isotope analyses, samples were thoroughly dried in a vacuum-drying cabinet at 60 ^∘C. Isotope values were calibrated with international (IAEA-CH-7, USGS41) and laboratory standards (peptone). The analytical precision was equal to or better than 0.2 ‰. Carbon (δ¹³C) values of samples were calculated by comparison with isotope-ratio-predetermined peptone and cellulose lab standards and certified international isotope standards (IAEA-601, IAEA-602, IAEA-CH-7, USGS41), which were inserted frequently in the course of sample measurements. Isotope ratios are presented in the common δ notation against PDB as

$${\mathit{\delta }}^{\mathrm{13}}\text{C}=\left[{\displaystyle \frac{\left(\frac{{\mathrm{13}}_{\mathrm{C}}}{{\mathrm{12}}_{\mathrm{C}}}\right)\mathrm{sample}}{\left(\frac{{\mathrm{13}}_{\mathrm{C}}}{{\mathrm{12}}_{\mathrm{C}}}\right)\mathrm{PDB}}}-\mathrm{1}\right]\times \mathrm{1000}\left(\mathrm{‰}\right).$$

The final tree-ring carbon isotope chronology was corrected to incorporate isotopically light carbon released by the burning of fossil fuels and increasing CO₂ concentration, as proposed by McCarroll and Loader (2004) and McCarroll et al. (2009). The correction procedure applied here has the advantage of being objective, as it effectively removes any declining trend in the δ¹³C series post 1850 CE, which is attributed to the physiological response to increased atmospheric CO₂ concentrations (McCarroll et al., 2009).

2.4 Statistical analyses

Tree growth–climate relationships were analysed applying simple Pearson correlation analysis. The relationship with glacier mass balance was analysed using the pair-plot correlation package “PerformanceAnalytics” (Carl et al., 2010) in R (Tables S4, S5; Fig. S3). Correlations were computed between monthly temperature and rainfall data and the tree-ring δ¹³C chronology for a window from the previous year's October to the current year's September.

To better visualize the comparison between our reconstructed glacier mass balance time series and other hydroclimatic reconstruction series from the monsoon-dominated central Himalayan region (Fig. 1), data time series were standardized using Z scores and smoothed with 11-year fast Fourier transform to highlight common climate signals.

Based on correlation analysis, a linear-regression model (Briffa and Jones, 1990; Cook et al., 1994) was used to perform the reconstruction using the lm module in R (ggplot2 package; Wickham, 2016). The leave-one-out cross-validation method (LOOCV; Michaelsen, 1987) was used for the entire calibration period (1982–2013) and to verify the reconstruction (Table S3). This method is most suitable when the length of observed records is short (Shah et al., 2013; Shekhar et al., 2017; Yadav and Bhutiyani, 2013; Vehtari et al., 2017; Zhang et al., 2019). In this method, each observation is successively withdrawn, a model is estimated on the remaining observations, and a prediction is made for the omitted observation. The LOOCV analysis was performed using the package “caret” (Kuhn et al., 2015). Rigorous statistics, including the sign test, reduction in error (RE), and correlation coefficients, were calculated to evaluate the similarity between observed and estimated values. The sign test measures the degree of association between two series by counting the number of agreements and disagreements. The series are highly correlated if the number of similarities is significantly larger than the number of dissimilarities. The RE statistic provides a rigorous test of the association between actual and estimated series. Any positive value indicates the predictive capability of the model. A positive RE is evidence of a valid regression model (Fritts, 1976). In addition, other rigorous statistics, viz. the root mean square error, coefficient of efficiency (CE), and Durbin–Watson (DW) test were carried out to evaluate the linear-regression model (Table S3).

Spatial and temporal correlations (moving correlations) were used to identify the coherence between reconstructed mass balance and gridded (${\mathrm{0.5}}^{\circ }\times {\mathrm{0.5}}^{\circ }$) temperature, precipitation, SST, ENSO, ENSO-Modoki, PDO, and IOD for the studied region in the central Himalaya. The variability in the climate and glacier mass balance reconstruction in the frequency domain was investigated using the multi-taper method (MTM) of spectral analysis (Mann and Lees, 1996) and wavelet transform (Grinsted et al., 2004) to identify periodicities and their temporal variability in the reconstructed data.

3 Results and discussion

3.1 Annual glacier mass balance reconstruction

The detailed results on the analysed relationship between δ¹³C and compiled glacier mass balance as well as of the relationship between δ¹³C and available climate datasets are presented in Tables S4 and S5. Descriptive statistics of δ¹³C chronologies and inter-species correlation (for the common period) are shown in Table S1. Uncorrected and corrected δ¹³C chronologies (McCarroll and Loader, 2004) of the PFTs are illustrated in Fig. S6. The mean difference in the isotopic composition of the two conifer species is ∼ 0.53 ‰ (Abies pindrow, 1743–2015, $-\mathrm{22.19}\pm \mathrm{0.7}$ ‰; Picea smithiana, 1920–2015, $-\mathrm{22.72}\pm \mathrm{0.65}$ ‰). The mean δ¹³C value of broadleaf deciduous species (Aesculus indica) is $-\mathrm{24.31}\pm \mathrm{0.82}$ ‰. The difference (1.6 ‰–2.1 ‰) between δ¹³C time series of two PFTs indicates a higher level of isotope discrimination in the broadleaf species relative to conifers. Further, break-point analysis indicated 1954 as a year of change in the isotopic composition. Irrespective of the PFTs, a mean ∼ 5 % decline in isotopic composition was noted after 1954 relative to the pre-1950s level (Fig. S6). For the broadleaf species, slopes of the trends for the periods prior and after 1954 are −0.001 and −0.032 ‰ yr⁻¹, respectively. For the conifers, they are −0.0007 and −0.025 ‰ yr⁻¹, respectively (P<0.001). Δ¹³C time series of broadleaf species showed an increasing trend from 1954, with an average rate of 0.15 ‰ yr⁻¹, while it remained stable for the conifers (0.003 ‰ yr⁻¹) (Fig. S6). Nevertheless, a significant and positive inter-species correlation exists between them. This indicates an influence of common and coherent climatic factors on the physiological processes. Climate-response functions indicate that for both PFTs, the correlation strength between δ¹³C and mean annual temperature (MAT) or mean annual precipitation (MAP) remained similar (MAT – Abies pindrow, $r=-\mathrm{0.215}$, P<0.05, n=66; Aesculus indica, $r=-\mathrm{0.281}$, P<0.05, n=66; Picea smithiana, $r=-\mathrm{0.405}$, P<0.001, n=66), (MAP – Abies pindrow, r=0.243, P<0.05, n=66; Aesculus indica, r=0.115, P<0.05, n=66; Picea smithiana, r=0.185, P<0.05, n=66). Results on monthly and seasonal climate-response-function analyses are provided in the Supplement (Fig. S3, Table S5).

A weak correlation between TRW and glacier mass balance could arise due to sensitivity issues (Bunn et al., 2018). At high elevations in a valley environment, temperature and precipitation signals mix and make TRW records difficult to interpret (Bunn et al., 2019). In view of this, we resorted to tree-ring isotopes (δ¹³C) of the species that are known to be sensitive to climate. Interestingly, we found a strong correlation only between δ¹³C of our studied conifer species and compiled glacier mass balance (Abies pindrow, r=0.596, P<0.001, n=23; Picea smithiana, r=0.631, P<0.001, n=23). The correlations are strong enough to establish a significant calibration model (Shah et al., 2013; Shekhar et al., 2017; Yadav and Bhutiyani, 2013; Vehtari et al., 2017; Zhang et al., 2019). In contrast, correlation with the deciduous species (Aesculus indica, r=0.343, P=0.1101, n=23) is weak during summer (April–September) and non-significant during winter. This could be due to a lack of storage of wintertime climate signals, as the species remains physiologically active only during the wet–warm period from April to September. Therefore, we utilized both evergreen conifer species for our reconstruction (Table S4). A linear-regression model was employed for the reconstruction of annual glacier mass balance (GMB) over the past 273 years, and the corresponding empirical equation is

$$\text{GMB}=\mathrm{8.3984}+\mathrm{0.3816}\times {\mathit{\delta }}^{\mathrm{13}}{\text{C}}_{\text{conifers}}.$$

Here, δ¹³C_conifers is the mean chronology of Abies pindrow and Picea smithiana and GMB is annual glacier mass balance (metre water equivalent, m w.e.). The detailed model statistics are presented in Table S3. Model and calibration–verification statistics indicate the reliability and strength of our reconstruction model (Table S3, Fig. S4; Fritts, 1976; Vehtari et al., 2017). Validation tests including the number of sign agreements between the reconstructed series and observed mass balance records and the cross-correlation between reconstruction and measurements are significant (P<0.001) (Fig. S4). However, the error estimates are based on measured mass balance data of only 23 years (1982–2013, with a gap of a few years), so a possibility of uncertainty still exists in the reconstruction, particularly during the pre-observation period. Nevertheless, the use of a sensitive isotope chronology (δ¹³C) and the combination of two conifer species in this study may help to minimize several factors responsible for higher uncertainty, such as sensitivity issues, decreasing sample size, and temporal variability, which are unavoidable with TRW.

3.2 Three major phases in the mass balance dynamics

Historical records of glacier change are rare from the Himalaya. However, available studies provide evidence that the monsoon-influenced southeast Tibetan Plateau (TP) glaciers and those of the central Himalaya have responded synchronously to the change in climate. In contrast, glaciers in the western regions of the Himalaya–Tibet orogen have behaved asynchronously (Solomina et al., 2016; Owen et al., 2008; Kaspari et al., 2008; Liu et al., 2013; Hochreuther et al., 2015; Xu et al., 2012; Bräuning et al., 2006). Studies highlight the relative importance of two moisture delivery systems (westerlies and Asian summer monsoon) in driving regional glacier fluctuations; how variations in moisture delivery systems have changed at millennial to glacial–interglacial timescales and their impact is even perceptible at interannual to decadal timescales (Hou et al., 2017, Mölg et al., 2014).

In a gross regional sense, the current state of the knowledge suggests that since the last glacier advance (LIA 1500–1850 CE) glaciers are in a general state of retreat, whereby the role of regional climate in regulating regional glacier fluctuations has increased gradually (Mayewski and Jeschke, 1979; Mayewski et al., 1980). Therefore, across the Himalaya this period (1850s onwards) is characterized by the dominancy of retreat, advance, and/or standstill regimes (Mayewski and Jeschke, 1979; Mayewski et al., 1980; Bolch et al., 2012; Rowan, 2017; Solomina et al., 2016). Available records from the central Himalaya indicate a state of glacier retreat since the 1850s, regardless of the glacier type (Mayewski and Jeschke, 1979), with a recent acceleration in ice mass loss as indicated by remote sensing studies (Maurer et al., 2019; Bandyopadhyay et al., 2019).

Figure 2 (a) Reconstructed glacier mass balance for the western central Himalaya (WCH) since 1743 CE (coloured inset indicates annual and pentadal mass change rate since the 1960s with SEs). (b) Annual snow accumulation record derived from Dasuopu ice core from the eastern central Himalaya (Thompson et al., 2000). Dark lines indicate 11-year moving averages. Vertical dashed lines are results of break-point analyses that reveal three major phases in the reconstruction.

Our 273-year-long GMB reconstruction agrees with the general retreat of glaciers since the mid-19th century (Fig. 2). Smoothing of the reconstruction with an 11-year moving average and break-point analyses revealed three distinguishable main phases: (1) a phase of positive mass balance up to ∼ 1870 CE, coincident with the LIA in the central Himalaya, (2) reorganization of Northern Hemisphere atmospheric circulation following the LIA and a progressively increasing influence of regional climate could have resulted in a phase of near-zero or stable mass balance that lasted up to the 1960s, and (3) a phase of accelerated ice mass loss since the 1960s. The latter phase corresponds to a global glacial retreat and can be attributed to increasing temperatures, combined with a decline in ISM and westerly circulations and with anthropogenic climate change (e.g. Bollasina et al., 2011).

Hydroclimatic evidence such as regional and composite tree-ring δ¹⁸O chronologies from the monsoon-influenced region (from southeast TP to central Himalaya) (An et al., 2014; Singh et al., 2019; Grießinger et al., 2011; Wernicke et al., 2017; Liu et al., 2014; Sano et al., 2012, 2013, 2017; Shi et al., 2012; Xu et al., 2018) and speleothems and the Dasuopu ice core record from the central Himalaya (Thompson et al., 2000; Kaspari et al., 2008; Denniston et al., 2000; Kotlia et al., 2012, 2015; Liang et al., 2015) show that regional climate has changed since the 1860s, with a reorganization of hemispheric atmospheric circulation. Consequently, the regional hydroclimate abruptly shifted towards a drier phase concurrent with the changes in atmospheric circulation. These records show a spatially coherent signal and serve as a validation test of the accuracy of our GMB reconstruction (Fig. 2). Particularly, tree-ring and speleothem δ¹⁸O records (Singh et al., 2019, Sano et al., 2012, 2013, 2017; Xu et al., 2018; Kotlia et al., 2015; Liang et al., 2015) from our study region (Fig. 1) indicate increased westerly precipitation prior to 1860–1870 CE (Murari et al., 2014; Yang et al., 2008), which resulted in a phase of positive mass balance. Increased winter snowfall could be anticipated from a southward shift in southwesterly winds due to reduced Northern Hemisphere air temperatures during the LIA (Rowen, 2017). The stronger influence of mid-latitude westerlies in driving glacier variability in monsoonal high Asia (Mölg et al., 2014) may have led to higher snowfall and positive mass balance prior to the 1870s as observed in our reconstruction (Fig. 2). Regional temperature reconstructions (Zhu et al., 2011; Borgaonkar et al., 2018; Yadav et al., 2011) also suggest a cold and cloudy climate prior to 1850, which was followed by a warmer and sunnier climate thereafter (Liu et al., 2014; Xu et al., 2012). Indeed, a brief phase of a negative GMB trend during the 1770s to 1790s could be the result of warmer winter temperatures (Huang et al., 2019a) and below-average snow accumulation during mega-droughts as observed around this period (Cook et al., 2010; Thompson et al., 2000; Kaspari et al., 2008) (Fig. 2). The close correspondence between reconstructed mass balance and regional hydroclimate reconstructions supports the notion of hemispheric synchronicity to climate change prior to 1860–1870 CE (Solomina et al., 2016). However, gradual changes in regional oscillations in temperature and moisture delivery sources in the Himalaya–Tibet orogen, particularly a change in the interplay of the ISM and westerlies could have induced regionally distinguished climate zones with specific behaviour of glacier dynamics (Sect. 3.3).

Glacial history during the late Holocene and LIA in our studied region (Murari et al., 2014; Saha et al., 2019), and the above-mentioned proxy records and climate dynamics studies (Mölg et al., 2014; Khan et al., 2019) indicate an out-of-phase relationship between the westerlies and ISM influence. Therefore, regardless of the strength of the westerlies (Joswiak et al., 2013), ISM weakening since the 1860s has tended to shift the regional hydroclimate towards a drier and a warmer climate. This could have favoured ice mass loss that led to a slightly negative mean mass balance ($-\mathrm{0.046}\pm \mathrm{0.134}$ m w.e., ±SD) during the middle phase (1870–1959) in our reconstruction. The timing of our reconstructed glacier advancement or recession is consistent with records of cold periods during the early 19th (1815–1825) and 20th (1900–1910) century (Zech et al., 2003). A brief period of declining mass balance from 1865 to 1885 could be ascribed to warm winter temperatures observed during 1848 to 1891 (Huang et al., 2019a) and the late-Victorian intense drought during 1875–1878 (Singh et al., 2018; Cook et al., 2010). The slightly negative but stable mean GMB observed from 1920 to 1960 is consistent with reported mass balances from the Himalaya (Bolch et al., 2012) and pluvial conditions from enhanced ISM activity during 1920 to 1960 (Sinha et al., 2015).

Severe glacial mass loss ($-\mathrm{0.437}\pm \mathrm{0.191}$ m w.e.) since the 1960s corresponds to a further and prominent shift towards drier conditions arising out of the weakening of both the ISM and westerlies (Qin et al., 2000; Hunt et al., 2019; Basistha et al., 2009; Singh et al., 2019; Sano et al., 2012, 2013, 2017). Increasing regional and global temperatures coupled with progressive regional industrialization may have exacerbated the retreat rates. We observed a doubling of average ice mass loss rate ($-\mathrm{0.577}\pm \mathrm{0.022}$ m w.e. yr⁻¹) (±SE) in the last 30 years compared to 1960–1985 ($-\mathrm{0.275}\pm \mathrm{0.022}$ m w.e. yr⁻¹). An extensive remote-sensing-based study encompassing the last 40 years observed a consistent and similar trend of glacier loss across the Himalayan transect (Maurer et al., 2019). For our study region, Bandyopadhyay et al. (2019) computed mean weighted mass balance as $-\mathrm{0.61}\pm \mathrm{0.04}$ m w.e. yr⁻¹ (2000–2014), while we observed $-\mathrm{0.65}\pm \mathrm{0.02}$ m w.e. yr⁻¹ (excluding the very large-sized Gangotri Glacier). Similarly, several geodetic mass balance estimates from the western central Himalayan region are comparable (−0.4 to −0.8 m w.e. yr⁻¹) to our results and are within the error limits of earlier studies (DCCC, 2018; Bandyopadhyay et al., 2019; Maurer et al., 2019).

In recent decades, along with the decreasing intensity of moisture influx (ISM and WD), the importance of regional temperature in regulating mass balance behaviour has been gradually increasing (Fig. 3). Several studies from the eastern Himalaya–Tibet orogen have emphasized the increased influence of regional temperature in determining mass balance behaviour (Bräuning, 2006; Yang et al., 2008). Our results confirm historical observations made in the early 19th century (without the availability of geochronological data to constrain the timing) that the central Himalayan glaciers have been receding since 1850 CE (Mayewski and Jeschke, 1979; Rowen, 2017). The rate of ice mass loss has accelerated since the 1960s and almost doubled in the last 30 years with respect to pre-1985. Here, our results contradict the notion of UNEP (2009) that “despite the widespread shrinkage of the Himalayan glaciers in area and thickness, the nature of shrinkage has not changed significantly over the last 100 years” (UNEP, 2009).

Figure 3 (a–c) Coherence in the anomalies (Z scores) of glacier mass balance (a), summer monsoon precipitation (Singh et al., 2019) (b), and CRU temperature (c). Dark lines indicate 11-year moving averages. (d) Enhanced correlation between glacier mass balance and temperature after the 1960s and the effect of the global warming hiatus of the late 1990s. (e, f) Correlations of glacier mass balance with precipitation (Singh et al., 2019) (e) and temperature (HadCRUT4) (f) running over 51 years. (g–l) Anomalies (Z scores) of six central Himalayan tree-ring δ¹⁸O series (Manali: Sano et al., 2017; Uttarakashi: Singh et al., 2019; Jageshwar: Xu et al., 2018; Humla: Sano et al., 2012; Ganesh: Xu et al., 2018; Bhutan: Sano et al., 2013). Dark lines denote 21-year moving averages. (m–r) Corresponding right panels indicate low-frequency temporal correlations (51-year running correlations) with reconstructed glacier mass balance. Different behaviour of regional δ¹⁸O chronologies, particularly a phase shift in WCH chronologies (Manali, Uttarakashi, and Jageshwar), is prominent after the mid-20th century (dashed vertical line).

3.3 Increasing regional heterogeneity

Substantial evidence exists from the monsoon-influenced Himalaya–Tibet orogen (Fig. 1) suggesting that LIA glacial fluctuations were more or less synchronous prior to 1850s (Hochreuther et al., 2015; Bräuning, 2006; Yang et al., 2008; Xu et al., 2012; Mayewski and Jeschke, 1979; Rowen, 2017; Owen et al., 2008). Following the LIA, a global readjustment in atmospheric circulation resulted in a southward shift of the Intertropical Convergence Zone (ITCZ). It broadly affected the Asian summer monsoon, and the region has progressively moved towards a drier phase, subsequent to a breakdown in the ITCZ and Northern Hemisphere temperature relation since the mid-19th century (Ridely et al., 2015; Xu et al., 2012). Nevertheless, the strength of the westerlies in the region remained unaffected until the mid-20th century (Joswaik et al., 2013; Hunt et al., 2019; Khan et al., 2019). Thus, depending upon the geographical location and proximity to the oceans, a distinct regionality developed, with a decline in summer monsoon precipitation. The middle phase (mid-19th to mid-20th century) in our reconstructed GMB and in the annual snow accumulation recorded in the Dasuopu ice core supports this point of view (Fig. 2). During this period, cryospheric mass balance turned slightly negative, with a gradual decline in summer precipitation, which was more prominent towards the western part (Figs. 2, 3).

Dendroglaciological and palaeoclimate studies suggest that on a centennial timescale, temperature changes rather than precipitation changes remain the prime factor for glacier fluctuations (Yang et al., 2008; Wang et al., 2019). After the mid-20th century (1960s), the role of temperature increased in determining mass balance behaviour, concurrent with a further decline in moisture influx both from the ISM and WD (Roxy et al., 2015, 2017; Qin et al., 2000; Hunt et al., 2019; Basistha et al., 2009; Singh et al., 2019; Sano et al., 2012, 2013, 2017). Correlation analysis indicates a high coherence of GMB dynamics with CRU temperature data after the 1960s ($r=-\mathrm{0.78}$, P<0.001). Interestingly, reconstructed GMB even sensitively responds to the effect of the slowdown in global temperature increase since the late 1990s (warming hiatus) (Fig. 3a). In contrast, correlations with gridded precipitation including that of northern India and ISM rainfall indicate a low association. However, we observed a high coherence (at both low and high frequencies) with local δ¹⁸O-reconstructed summer monsoon precipitation derived from several regionally dominant tree species (Singh et al., 2019) (Fig. 3a). Moreover, we also found a tight correspondence with regional (Fig. 1) tree-ring δ¹⁸O chronologies; the strength of which declined towards the eastern central Himalaya (Fig. 3b, Table S6). Results indicate an abrupt change in hydroclimate after the mid-20th century, which is particularly prominent in the western part of the central Himalaya. Moving correlations (51-year) between GMB and tree-ring δ¹⁸O chronologies indicate a phase shift in the western part, which is more prominent relative to the eastern central Himalaya after the 1960s (Fig. 3b).

Regional tree-ring δ¹⁸O chronologies constitute sensitive records that exemplify regional hydroclimatic heterogeneity and show that the drivers of glacier fluctuations in the western part of the central Himalaya (WCH) are somewhat different from those in the eastern central Himalaya (ECH) (Fig. 3b). Here, we contend that caution should be applied when referring to the glaciers in the WCH as summer-accumulation type glaciers. We showed that over the WCH (compared to the ECH), the westerlies still have a significant impact on annual mass balance behaviour. Earlier studies confirm that glaciers in the ECH and further east are mainly fed by summer monsoon precipitation and thus are undisputedly classified as summer-accumulation type glaciers (Sakai and Fujita, 2017). However, our results reveal a contrasting hydroclimate relation between the WCH and ECH, with an opposite behaviour of moving correlation patterns (51-year) between records of regional glacier–snow accumulation and tree-ring δ¹⁸O series (Fig. 4). When comparing the δ¹⁸O chronology of the deciduous species (Aesculus indica, growing in the WCH), we find the correlation pattern is strikingly similar to that of the ECH, rather than to the WCH. In contrast to conifer species that take up meltwater containing winter precipitation isotopic signals in the early growing season (Huang et al., 2019b), the broadleaf deciduous species takes up soil water mostly after leaf flush and full canopy development in May (Negi, 2006). At this time of the year, the snowmelt signal may already be gone, so the species does not reflect the winter precipitation signal. These results support our interpretation that the winter westerlies have a strong impact on the annual mass balance of glaciers in the WCH compared to in the ECH. The correlation patterns indicate that there has been an abrupt phase change in the WCH since the 1960s, while correlation in the ECH has remained stable, suggesting a super-positioning role of temperature and differential precipitation seasonality derived from the westerlies and the two branches of the ISM (AS and BoB) (Figs. 2, 4).

Figure 4 Low-frequency correlations between regional cryospheric dynamics and hydroclimate. Blue line indicates 51-year moving correlations between reconstructed glacier mass balance (GMB) and regional tree-ring δ¹⁸O (Manali, Uttarakashi, and Jageshwar) from western part of the central Himalaya (WCH). The abrupt phase shift since the 1960s is peculiar for the WCH. Green line indicates 51-year moving correlations in the eastern central Himalaya (ECH) between annual snow accumulation from the Dasuopu ice core (Thompson et al., 2000) and regional δ¹⁸O (Ganesh and Bhutan), while the red line indicates 51-year moving correlations between GMB and tree-ring δ¹⁸O of a dominant deciduous species (Aesculus indica) growing in the WCH. Its similarity in coherence to that of the ECH is remarkable and illustrative of the summer monsoon influence during the annual growing period (April–September).

Analyses of tree-ring δ¹⁸O records and water stable-isotope ratios from the central Himalaya specify a greater influence of the Arabian Sea branch of the ISM in the WCH, the strength of which declines towards the eastern Himalaya (Sano et al., 2017), while the climate in the ECH is dominantly regulated by the Bay of Bengal branch, along with some influence of the East Asian summer monsoon (Sano et al., 2013, 2017; Sinha et al., 2015; Liu et al., 2014). At the northwestern periphery of ISM incursions in the WCH, moisture flux from the Arabian Sea competes with the influence from the Bay of Bengal (Sano et al., 2017). Thus, the role of the Arabian Sea branch increases in the WCH, particularly when the Bay of Bengal branch is relatively weak (Sano et al., 2017).

It is pertinent here to mention the differential influence and behaviour of the two branches of the ISM. Although the Arabian Sea (AS) and the Bay of Bengal (BoB) are part of the Indian Ocean, crucial ocean–atmosphere interaction experiments have shown that they possess distinctly different features (Saikranthi et al., 2019). They are strongly different in terms of sea surface temperature (SST) and background atmosphere and related precipitating systems. The monsoonal winds and low-level Findlater jet are stronger over the AS than over the BoB. Tropospheric thermal inversions are more frequent and stronger over the AS than over the BoB. The variability in SST is larger over the AS than over the BoB. The SST in the AS cools between 10 and 20^∘ N during the monsoon season, whereas warming occurs in all oceans between the same latitudes (Misra et al., 2019; Saikranthi et al., 2019; Roxy et al., 2015). Thus, the AS plays a predominant role in regulating the ISM rainfall variability, which appears to be particularly impactful for the prevailing climate in the WCH.

Compelling evidence suggests there has been increased warming over the Indian Ocean since the mid-20th century. However, warming over the AS has been monotonous for more than a century, at a rate faster than in any other region of the tropical oceans (Roxy et al., 2015, 2017; Misra et al., 2019). Currently, the AS is the largest contributor to the trend in global mean SST. The abnormal warming over the AS tends to weaken the land–ocean thermal contrast and hence has been implicated in a weakening of the ISM since the mid-20th century (Roxy et al., 2015, 2017). Besides, warming over the AS has also been shown to reduce the magnitude of the westerlies, as their interaction enhances the moisture convergence over the AS, leading to decreased precipitation from the westerlies (Misra et al., 2019). This mechanism may explain the decline in precipitation frequency, snowfall, and total precipitation amount from the westerlies since the mid-20th century (Hunt et al., 2019; Shekhar et al., 2010; Kumar et al., 2015; Khan et al., 2019). Tree-ring δ¹⁸O records from Bhutan to the southeast TP (Sano et al., 2013; Hochreuther et al., 2016; Lyu et al., 2019) show more or less unaltered conditions during the 20th century. Thus, weakening of the BoB branch as indicated in some studies (Roxy et al., 2017) needs reappraisal. However, irrespective of the strength or variability in the Bay of Bengal monsoon branch, accelerated ice mass loss since the 1960s could possibly be ascribed to the combined effect of a decline in moisture influx from the Arabian Sea and the westerlies.

Figure 5 (a–c) Temporal correlation behaviour at a multi-decadal timescale (51 years) between tree-ring δ¹⁸O chronologies available from the western central Himalaya (Manali: Sano et al., 2017; Uttarakashi: Singh et al., 2019; Jageshwar: Xu et al., 2018), Niño 4 SST (violet line) (a), ENSO Modoki (orange line) (b), and the Pacific Decadal Oscillation (PDO; red line) (c). (d) Strengthening decadal correlations between glacier mass balance and oceanic indices after the 1960s. (e) Spectral-analysis plot of glacier mass balance (GMB) time series reflecting low-frequency variations preserved in GMB reconstruction.

Given the importance of the westerlies (Mölg et al., 2014) and their weakening, as reflected in the decline in the indices of the Arctic Oscillation and North Atlantic Oscillation (Hunt et al., 2019; Peings et al., 2019), we focus on the El Niño–Southern Oscillation (ENSO), the Indian Ocean Dipole (IOD), and the Pacific Decadal Oscillation (PDO), which are all known to strongly modulate the summer monsoon precipitation. Hydroclimatic studies from the monsoon-dominated Himalaya–Tibet orogen indicate a pervasive but temporally unstable coherence. Low-frequency coherence between tree-ring δ¹⁸O-derived regional hydroclimatic reconstructions and the above indices of atmosphere–ocean interaction appears to have strengthened in recent decades (Sano et al., 2012, 2013, 2017; Singh et al., 2019; Hochreuther et al., 2016; Lyu et al., 2019). In accordance, our results too show a strong correlation at a multi-decadal timescale between tree-ring δ¹⁸O chronologies available from the WCH and the above indices (Fig. 5).

Tree-ring δ¹⁸O chronologies from the monsoon-dominated Himalaya generally contain a 2–5-year cyclic signal, which is coherent with the ENSO cycle (Xu et al., 2018, Singh et al., 2019; Sano et al., 2012, 2013, 2017). However, spectral analysis of our GMB time series showed significant periodicities (P<0.05) at decadal and multi-decadal scales, which reflect low-frequency variations preserved in the reconstruction (Fig. 5e). The lack of high-frequency periodicities indicates a minor role of ENSO in modulating the glacial hydroclimate on interannual or shorter timescales. Low periodicity was also found in a regional reconstruction of mass balance (Shekhar et al., 2017). Cross-correlations between reconstructed GMB and Niño 4 SSTs revealed a weak negative correlation ($r=-\mathrm{0.25}$, P<0.05, 1901–2015). Interestingly, weak and non-significant negative correlations between a δ¹⁸O-based summer monsoon rainfall reconstruction for our study region and SST (Singh et al., 2019) during 1909–1926 and 1957–2015 coincide with a negative GMB trend (Fig. 2), while significantly strong negative correlations observed during 1870–1908 and 1927–1956 coincide with positive GMB trends (Fig. 2). This indicates a clear influence of SSTs on mass balance behaviour particularly for the WCH at low-frequency timescales. Wavelet coherence analyses of reconstructed GMB with annual Niño 4 SST revealed an unstable common frequency composition in the different series over time, indicating a wobbly linkage with the change in SSTs (Fig. S5), as observed in previous studies of regional hydroclimates (Sano et al., 2012, 2013, 2017; Singh et al., 2019; Hochreuther et al., 2016; Lyu et al., 2019). Similarly, we noted a weak correlation with other indices (Fig. S5). However, decadal correlations between GMB and these indices strengthen after the 1960s (Fig. 5), indicating a negative influence of monotonous ocean warming on mass balance dynamics through precipitation feedback.

A lack of strong coherence in the high-frequency bands between GMB dynamics and SST indices is possible because once the monsoon circulation sets in, the local weather begins to decouple from the circulation. There is a complete absence of a monsoon footprint on regional glaciers both from atmospheric drivers and from mass balance perspectives (Mölg et al., 2014). Point-scale micrometeorological experimentation and surface energy balance analyses of a typical glacier from the WCH (Pindari Glacier, Fig. 1) also indicate a tight land–atmosphere coupling during the summer-accumulation season as well as in the winter-accumulation season. The coupling turned negative and remained weak only during the transition phases (pre- and post-monsoon), making it amenable to atmospheric perturbances (Singh et al., 2020). Thus, any repeated atmospheric alterations during the summer monsoon season such as active or break periods may not be able to strongly influence the mass balance dynamics. However, critical transition phases, particularly the pre-monsoon melt season, have the potential to shape mass balance on decadal or even longer timescales (Mölg et al., 2014; Singh et al., 2020).

The pre-monsoon season is also important for the onset of the summer-accumulation season (during the ISM). Numerical experiments and simulations in the Himalaya–Tibet orogen have identified seasonal changes during the pre-monsoon season that lead to the onset of the ISM. The Tibetan Plateau maintains a large-scale thermally driven vertical circulation, which initially remains unconnected to the monsoon circulation. Rising motion exists only on the far western side of the plateau during winter. Gradually, it extends towards the eastern side as the pre-monsoon season progresses. The ISM onset over the region is thus an interaction process between the plateau-induced circulation and northward-migrating rain-bearing depressions from the AS and BoB. It is thus highly possible that warming-induced decline in winter–spring snow cover, particularly over the western TP, could alter this interaction process and monsoon onset. The weakening of the relationship between snow cover over the western TP and ISM rainfall could be a possible manifestation of this (Zhang et al., 2019). Studies have affirmed that atmospheric conditions during the pre-monsoon season are pivotal in determining mass balance even at decadal or longer timescales (Mölg et al., 2014). These studies also suggest anti-correlated precipitation variability between the ISM and mid-latitude circulation during the pre-monsoon season. The interaction between these large-scale circulations in this region is largely determined by the position of thermally driven subtropical jet streams (Peings et al., 2019; Hunt et al., 2019). Thus, a possibility emerges that anthropogenic factors influencing tropospheric warming could alter the position of jet streams and annual mass balance behaviour, particularly in the transitional climate zone in the WCH.

3.4 Local forcing factors

Our study region has experienced an enhanced warming trend during the pre-monsoon season since the 1960s that has been amplified in the last 2–3 decades (Lau et al., 2010; Gautam et al., 2009, 2010; Xu et al., 2009; Prasad et al., 2009). This is most likely due to increased aerosol loading, enhanced atmospheric water vapour content (Mortin et al., 2016), and anthropogenic greenhouse gas (CO₂) emissions. These factors increase radiative heating rates, aid net atmospheric warming, and create a positive feedback loop with tropospheric warming. Enhanced tropospheric warming has the potential to alter the position of the subtropical westerly jet and, as described earlier, determines the ISM strength by influencing the AS branch and ultimately the glacier mass balance.

Of the local factors influencing tropospheric heating, greenhouse forcing of CO₂ emissions is pervasive. Studies suggest an inverse geological relation between glaciation events and atmospheric CO₂ concentration (Willeit et al., 2019; Macdonald et al., 2019). Causes of enhanced CO₂ concentration over the Himalaya include deforestation, increased frequency of forest fires, regional industrialization, and emission and transport of greenhouse gases. Besides major contributions from anthropogenic activities, sub-glacial carbonate weathering, especially the sulfide oxidation process coupled to carbonate carbonation, is a potential source of CO₂ to the atmosphere. Fast kinetics associated with these processes in glacial environments are also responsible for present-day atmospheric CO₂ concentrations (Shukla et al., 2018). Our corrected chronologies (Δ¹³C) (McCarroll and Loader, 2004) showed an effect of rising atmospheric CO₂ concentrations on tree physiology, which was prominent for the deciduous trees (Fig. S6). The current regional annual radiative forcing due to CO₂ is about 1.7 W m⁻², which is almost 5-fold that of the early 20th century. A peak in the transport of dust and greenhouse gases from the Indo-Gangetic Plain occurs during the pre-monsoon season, and we anticipate enhanced seasonal CO₂ radiative forcing during this critical transitional phase.

Climate sensitivity studies suggest that surface temperature change is higher for aerosols than for greenhouse gases due to their strong longwave effect (Chakraborty and Lee, 2019). The severe lack of aerosol data precludes the drawing of any conclusions on black-carbon-induced ice and snow melting in the Himalaya. Some space-borne studies utilizing ground radiometric measurements indicate that absorbing aerosols have favoured a localized warming over the western central Himalaya (0.26±0.09 ^∘C per decade) (Gautam et al., 2010). Related studies from Tibetan glaciers suggest that black carbon aerosol deposition is a significant factor contributing to the observed rapid glacier retreat (Xu et al., 2009; Zhao et al., 2017).

To obtain an estimate of radiative forcing due to ambient concentrations of black carbon aerosols over the glacier environment, we utilized data of two existing aethalometer stations from our study region (Fig. 1; Gangotri Glacier, 30^∘58^′48^′′ N, 79^∘4^′48^′′ E; 3600 m a.s.l.; Negi et al., 2019; and Dokriani (DOK) Glacier, 30^∘51^′6.86^′′ N, 78^∘45^′8.55^′′ E; 3900 m a.s.l.). Annual mean radiative forcing (at atmosphere) due to black carbon aerosols (computed as described by Dumka et al., 2016) was $+\mathrm{10.1}\pm \mathrm{3.0}$ W m⁻² for the Dokriani Glacier and $+\mathrm{7.4}\pm \mathrm{2.2}$ W m⁻² for the Gangotri Glacier, respectively. The lower value for the Gangotri Glacier could be due to its location on the leeward side. In support of these results, we found similar mean annual values and diurnal to seasonal dynamics in aerosol radiative forcing to those of sites on the southeastern TP (Zhao et al., 2017; Li et al., 2018; Zhang et al., 2017). Although we were not able to compute seasonal radiative forcing arising out of black carbon aerosols and/or dust, we anticipate manyfold-enhanced pre-monsoonal radiative forcing. These observations indicate the significant contribution of the local forcing factors associated with anthropogenic climate change.

4 Conclusions

We present a 273-year-long ice mass loss record from the Uttarakhand Himalaya, which is the only record of annual mass balance variability from the Himalaya so far. Utilizing coherent analyses of tree-ring stable-isotope chronologies of regionally dominant tree species and synthesizing oxygen isotope chronologies from different regional archives, we show that ice mass loss in the region has accelerated since the mid-20th century to its highest levels during the past 273 years. Stable-isotope chronologies of three species of diverse plant functional types and well-calibrated δ¹³C-derived mass balance reconstruction indicate a mean rate of mass loss of $-\mathrm{0.437}\pm \mathrm{0.025}$ m w.e. yr⁻¹ since 1960. The reconstructed mass balance reveals three major phases: mass balance remained positive prior to the 1860s to 1870s; there was a slightly negative but stable mass balance up to the 1960s; and there has been a highly negative mass balance since then. Further, our results suggest a doubling of average ice mass loss rate ($-\mathrm{0.577}\pm \mathrm{0.021}$ m w.e. yr⁻¹) in the last 30 years as compared to 1960–1985 ($-\mathrm{0.275}\pm \mathrm{0.022}$ m w.e. yr⁻¹). Isotopic and climate coherency analyses indicate that reconstruction is consistent with regional climatic variability and indicate a significant influence of the westerlies in this transitional climate zone of the central Himalaya. The current trend in ice mass loss goes hand in hand with increased anthropogenic aerosols and CO₂ loading and with a concurrent change in climatological factors, including a decline in the westerlies and in the strength of the Arabian Sea branch of the Indian summer monsoon. These results present observational support to calibrate and validate coupled regional climate–glacier models. Improvements in glacier mass change assessments are still possible and necessary. We suggest that incorporating growth increments of shrubs growing in the alpine zone could complement this work to provide robust estimates of glacier responses to future climate scenarios and modelling studies of glacier contributions to regional runoff.