Multi-decadal monsoon characteristics and glacier response in High Mountain Asia

The present study focuses on three debris-free, mountain glaciers distributed across distinct climatic sub-regions of HMA (figure 1). Yala Glacier (28.237^∘ N 85.619^∘ E) is a small, 1.5 km² glacier situated in Langtang Valley, Nepal which has been the subject of long-term mass balance campaigns (e.g. Fujita et al 1998, Baral et al 2014, Stumm et al 2021, Sunako et al 2020) as well as the focus of research on energy balance and surface processes (Stigter et al 2018, Litt et al 2019). Parlung Glacier Number 4 (29.245^∘ N 96.928^∘ E, hereafter 'Parlung 4' or 'PAR4') is a 11.7 km², spring-type accumulation glacier in the southeastern Tibetan Plateau. The glacier has been the focus of several recent studies that have explored meteorological processes (Yang et al 2011, Shaw et al 2021), glacier dynamics (Yang et al 2020) and energy/mass balance (Zhu et al 2015, 2018, Jouberton et al 2022). Mugagangqiong Glacier (34.248^∘ N 87.490^∘ E, hereafter 'MUGA' for tables and figures) is a 2 km² glacier in the arid central Tibetan Plateau. The glacier has been the focus of few studies, typically regarding ice cores and historic atmospheric dust concentrations (Feng et al 2020), though is also one of several investigation sites established across the Tibetan Plateau in the last decade (Yang et al 2022). The site specific information of each glacier is given in table 1.

We utilise hourly meteorological data from off- and on-glacier automatic weather stations (AWSs) at each study site (figure 1) to bias-correct the long-term forcing data. AWS data are available for the period 2012–2019 for Parlung 4 and Yala (figure S7), and 2015–2019 for Mugagangqiong. At each site, distributed ablation stake data are available, spanning partial (Parlung 4) or full elevation ranges (Yala/Mugagangqiong) of the glaciers.

Gridded atmospheric data for HMA (1980–2019) are extracted from simulations carried out with the weather research and forecasting (WRF) model at 9 km horizontal resolution (Ou et al 2020, Sun et al 2021). The WRF model (v. 3.7.1) is driven by three-hourly ERA5 pressure and surface level variables. Spectral nudging is applied to geopotential, horizontal winds and temperatures above the approximate planetary boundary layer top as described in Ou et al (2020). The model applies no convection scheme, as found by Ou et al (2020) to produce the most appropriate frequency and initiation timing for short (1–3 h) and long ($\gt$6 h) precipitation events compared with observations. The full model details are provided in table S1.

WRF data are extracted from the nearest corresponding point to each off-glacier AWS (figure 1) which provides a longer series of un-interrupted data for bias correction. Data gaps in the AWS records are ignored such that only hours with all energy balance variables available (air temperature, relative humidity, radiative fluxes, wind speed, air pressure and precipitation) are considered for bias correction. We apply a multivariate bias-correction following Cannon (2018) that combines quantile delta mapping with random orthogonal rotations to match the multivariate distributions of both the WRF and AWS data. This approach has been shown to outperform univariate bias correction methods when used in hydrological modelling applications (Meyer et al 2019, Faghih et al 2021) and we consider it suitable for correction of cross-correlated variables in complex, high mountain regions. Using at least 5 years of AWS data at each site, we generate quantile correction factors that are applied to the full time series of WRF data at the extracted pixel of interest. The variable correlations and long-term trends are preserved in the bias-corrected data and found to be consistent with long-term climate stations for the different sub-regions (see SI).

With the bias-corrected time series of WRF data, we extrapolate the variables of interest to the elevations of the glaciers (table 1) using empirically determined lapses rates. Multiple station observations available in Langtang Valley (Steiner et al 2021) were used to generate hourly-monthly temperature lapse rates, and temperatures were adjusted for the boundary layer cooling effect of Yala using the on-glacier AWS temperatures. For Parlung 4, we consider the off-glacier temperature lapse rates derived from the data of Shaw et al (2021) and as applied by Jouberton et al (2022). Precipitation is distributed similarly to Jouberton et al (2022) who applied a correction factor of 1.75 and gradient of 14% 100 m⁻¹ to account for the distinct precipitation regime between valley and the glacier. Given the differences in model forcing and structure, we recalibrate these gradients for the current study against ablation stakes and an upper elevation ice core record (table S3). No on-glacier AWS exists for Mugagangqiong, though the off-glacier AWS is very close to the glacier and is considered to experience a partial cooling effect associated with the glacier boundary layer (e.g. Shaw et al 2021). For this site, we distribute temperature considering the lapse rates derived for the Mugagangqiong basin (Yang et al 2022). Precipitation gradients and correction factors were calibrated against multi-year ablation stake observations. Use of a precipitation correction factor implicitly accounts for unknown quantities of AWS pluviometer under-catch that were not represented by the bias correction of WRF data.

We simulate hourly glacier energy and mass balance using the snow, ice and hydrological components of the Tethys-Chloris (T&C) land surface model (Fatichi et al 2012, 2021, Mastrotheodoros et al 2020, Botter et al 2021, Fyffe et al 2021, Fugger et al 2022) applied in a semi-distributed manner using 50 m intervals for the full range of elevations at each site (table 1). T&C's cryospheric component constitutes a fully-physical energy balance model that simulates the energy fluxes of the ice/snow surface and subsurface, including sublimation. Precipitation phase partitions are dynamically considered based upon the air temperature, relative humidity and air pressure (Ding et al 2014) and utilised directly in the albedo model developed by Ding et al (2017). We adjust the mean incoming radiative fluxes and snow avalanching according to the topography of each elevation band and the modelled elevational mass balances are area-weighted by time-evolving hypsometries from co-registered DEMs and glacier outlines since the 1970s to produce a continuous time-series of glacier mass balance (figures S16–S18). Avalanching on Parlung 4 Glacier was modelled given a fully-distributed TOPKAPI-ETH simulation of the catchment by Jouberton et al (2022) and used to update end-of-year mass balances for the appropriate elevation band. The model is run for 1980–2019 and the results analysed for 1981–2019, taking 1980 as a spin-up year. Temperature and precipitation distribution parameters are perturbed in 1000 Monte Carlo simulations to provide an estimate of mass balance uncertainty. More details about the model are provided in the supplementary information section 2.3.

As an estimate of the regional onset of the monsoon, we consider the horizontal wind shear index (HWSI) of Prasad and Hayashi (2005), as previously considered in studies of this type (Mölg et al 2012, Li et al 2018). Due to the large scale feature of the monsoon, we utilise ERA5 850 hPa zonal wind data to derive this index. Alternative means to classify local monsoon timing (e.g. Bombardi et al 2020, Brunello et al 2020) resulted in spatial-temporal inconsistencies between study sites for certain years and were thus not applied here. In this study, we analysed the results of the glacier energy- and mass-balance in relation to the monsoon (i) onset dates; (ii) cessation dates; (iii) duration and; (iv) intensity (figure S27). Moreover, we analyse the May–June precipitation amount and phase during the pre-monsoon to monsoon transition and examine its potential influence on the glacier energy and mass balance during the summer. The relevance of the monsoon characteristics is then placed into the context of the prevailing meteorological conditions for the summer. Early (late) onset/cessation dates are considered as those which are less (greater) than one standard deviation from the mean of the whole period and monsoon intensity is taken as the normalised mean HWSI for June–September.

The regional monsoon index reveals a tendency towards an earlier monsoon onset (−3.3 days dec⁻¹) and later cessation (+11.6 days dec⁻¹) between the 1980s and the early 2000s (figures 1(b) and S27), coinciding with a period of relatively strong monsoon intensity (figure 1(c)). Since 2005 the monsoon has arrived later and withdrawn sooner, shortening its total duration and becoming weaker in its intensity (figures 1(c) and S27). All sites see an increase in total precipitation coinciding with this expansion of the monsoon duration (figures 2(c)–(h)), especially for Parlung 4, which also demonstrates a period of relatively cool temperatures and increased cloudiness during the early 2000s, interrupting the overall warming trend (figure 2(d)). The entire HMA reveals a significant warming trend over the past four decades (figure S28), though with markedly different patterns at the study sites: summer (JJAS) temperatures at Yala continue almost unabated at a rate of 0.17 ^∘C dec⁻¹; the stronger warming at Parlung 4 (0.23 ^∘C dec⁻¹) combines a period of relative cooling and no clear trend pre-2000s, with a much stronger warming pattern since 2005, and; Mugagangqiong demonstrates a coincident increase in mean annual air temperature and precipitation during the early millennium (figure 2(e)).

Zoom In Zoom Out Reset image size

A higher intensity monsoon produces increased total precipitation for Parlung 4 and Yala, but has no significant correlation with precipitation at Mugagangqiong (figure 2(a)). An early monsoon arrival relates to increased precipitation (negative correlation) for much of southeastern and south central Tibet during May–June (figures 2(b) and S48), but shows no clear correlation for the Nepalese Himalaya (figures 2(b) and (f)). The increase in monsoon duration is governed by a tendency towards later cessation (figure S27) and has a more notable influence on summer temperature for the eastern Tibetan Plateau, though with few statistically significant correlations (figure S36). The patterns of radiative fluxes, cloudiness and humidity equally show few statistically significant trends across HMA for the full 40 year period (figures S28–S34), though at the glacier sites, the inter-decadal and long-term patterns suggest an increasing cloudiness and humidity alongside reduced shortwave radiation (figure S35).

The long term glacier mass balance of Mugagangqiong is highly distinct to the lower elevation, southerly sites (figure 3(a)) with a near-neutral, negative mass balance that has changed little since the 1980s. This is characteristic of its low temperatures, dry environment and low mass turnover rates (figure S44). Yala and Parlung 4, meanwhile, exhibit a much larger total mass loss over the past 40 years, despite significant differences in their inter-annual and inter-decadal mass balances (figure 3(a)). The mass balances of Yala and Parlung 4 closely reflect the inter-annual variability of air temperature and total precipitation (figure 2). Yala, however, experiences an almost continuous mass loss since 1980 with few positive mass balance years (figure 3(c)). Conversely, the mass loss at Parlung 4 is a combination of largely negative mass balance years interspersed by extended periods of positive or near-neutral mass balance (figure 3(d)) resulting in more mass loss during the drier 1980s (mean of −0.35 ± 0.11 m w.e. a⁻¹ compared to −0.24 ± 0.20 m w.e. a⁻¹ of Yala) followed by minimal mass loss from 2000 to 2005 (mean of −0.03 ± 0.13 m w.e. a⁻¹ compared to −0.28 ± 0.17 m w.e. a⁻¹ of Yala) and finally a period of more rapid mass loss continuing until 2019 (−1.00 ± 0.12 m w.e. a⁻¹ vs. −1.01 ± 0.19 m w.e. a⁻¹ at Yala). Patterns of mass balance at Mugagangqiong indicate no extended periods of negative or positive mass balance, but rather year to year fluctuations between weak positive and negative mass balances (figure 3(e)). On average, modelled glacier mass balances were −0.59 ± 0.20, −0.55 ± 0.12 and −0.04 ± 0.09 m w.e. a⁻¹ for Yala, Parlung 4 and Mugagangqiong, respectively.

Zoom In Zoom Out Reset image size

A key driver of the energy budget at all sites is the pre-monsoon and monsoon snowfall (figures 4 and S43), which dictates the increases in early monsoon albedo and reduction in net shortwave radiation (figure S47). While all sites behave similarly in this regard, the mass balance of Mugagangqiong shows a stronger correlation to both pre-monsoon and monsoon snowfall (figure S43), as its low melt rates and notable sublimation (constituting 20% of total ablation vs 8% at Yala and 3% at Parlung 4) are highly controlled by radiative fluxes under cold temperatures.

Zoom In Zoom Out Reset image size

Figure 4 emphasises the distinct patterns of snowfall and precipitation phase during the pre-monsoon to monsoon transition in May–June. At Yala, an early monsoon onset coincides with earlier snowfall at the glacier, though the total amount of snowfall is larger (∼37 mm on average) when the monsoon arrives late (figure 4(a)), which appears to be strongly linked to the higher fraction of solid precipitation for those years (figure 4(a)). Despite this average increase in May–June snowfall to the glacier, a late monsoon onset does not significantly alter the mass balance of Yala for the remainder of the year (figure 4(b)).

At Parlung 4, a later monsoon produces slightly cooler average conditions (figure S47), delaying liquid precipitation events during May compared to an early onset year (figure 4(c)). However, because of the substantial decrease in precipitation (∼130 mm on average), a late monsoon onset can sizeably influence the summer mass balance of Parlung 4 (figure 4), because the onset timing itself can explain 32(27)% of the variability in pre-monsoon (monsoon) snowfall. Due to the reduced surface albedo (figure S47—which can explain up to 87% of glacier mass balance variability) and heightened net shortwave radiation, mass balances are almost 0.5 m w.e. more negative on average than early onset years (figures 4 and 5), despite that only a weak correlation (−0.18, p = 0.2) to the monsoon onset itself exists (figure S43). While a late monsoon onset produces a consistent mass balance response for Parlung 4, early onset years result in a wide range of potential mass balances for the summer months (figure 4(d)). The summer and annual glacier mass balances are, however, significantly distinct between early and late onset years (p = 0.05, figure 4(d) boxplots).

Zoom In Zoom Out Reset image size

Interestingly, for years where the monsoon arrived late, May–June precipitation at Mugagangqiong was augmented (40 mm (58%) increase—figure 4(e)), suggesting a stronger role of westerly storm events in the absence of the monsoon precipitation (figure 1(a)—correlation = 0.27, p = 0.1). Given its elevation, however, no notable differences in precipitation phase are apparent for Mugagangqiong during this transition period (figure 4(e)). Average differences in cumulative mass balance are only apparent until late August, however, and are typically ∼0.04 m w.e. Variability in onset timing produces negligible correlations (−0.15) with the amount total monsoon snowfall at this site. As later monsoon onset years were typical for the cooler period of the simulation (figure S40), we also tested early vs late monsoon years only for the pre-2000s period and found a consistent pattern of May–June precipitation amount and phase for all glaciers as described above (figure S45).

Annual mass balances of Parlung 4 are generally affected less by the monsoon cessation date (figure 5). The mass balances of Yala and Mugagangiong are, however, more responsive to the cessation timing of the monsoon, but for different reasons: a later cessation at Yala relates to prolonged warmer air temperatures (figure S48(a)) which drives more melting and more negative mass balances; for Mugagangqiong, a later cessation date provides a longer period with precipitation (up to 200 mm on average compared to early cessation years—figure S48(f)), and thus produces a more positive mass balance. However, a more positive post-monsoon mass balance at Mugagangqiong does not significantly impact to the spring mass balance in the following year (not shown).

Monsoon duration itself balances many of the aforementioned processes to reveal no clear relationship to annual or summer mass balance (figure 5). The years with the weakest monsoon intensity, however, coincide with recent, warmer conditions (Yao et al 2022), producing a consistently more negative mass balance at all sites: up to 1 m w.e. more mass loss on average at Yala; 0.45 m w.e. greater mass loss at Parlung 4 and; 0.13 m w.e. more mass loss at Mugagangqiong relative to strong monsoon intensity years. Controlling for the effect of increasing air temperatures in recent years, the monsoon intensity is able to explain ∼50% of the variability in Yala summer mass balance (partial correlation r = 0.498, p = 0.001), but only 15% for Parlung 4 and $\lt$3% at for Mugagangqiong.

Monsoon characteristics (i.e. onset timing) can be associated with variable responses in glacier mass balance at our study sites, but the variability in meteorological conditions during the summer, which are only be partly influenced by the monsoon itself (up to 20% variability in (pre-)monsoon precipitation explained by monsoon timing or intensity—figure S43), play a dominant role in determining surface conditions and resultant mass loss. Figure 6 indicates that air temperature is the strongest control on glacier mass balance for Parlung 4, largely through its impact on the phase of monsoon precipitation (Jouberton et al 2022—figures 6(k)and S49(h)) and surface albedo (figure S49(k)—Zhu et al 2022). Controlling for the effect of air temperature itself, changes in the solid fraction of precipitation can explain 41% of the variability in glacier mass balance at Parlung 4. Meanwhile, precipitation phase changes had a smaller impact on Yala and Mugagangqiong, and mass balances were mostly responsive to the total amount of monsoon precipitation (figures 6, S49 and S50). Although the role of turbulent heat fluxes was more distinct between each site, increases in net shortwave radiation were consistently associated with the most negative mass balance years for all glaciers (figure S50).

Zoom In Zoom Out Reset image size