Positive mass budgets of high-altitude and debris-covered fragmented tributary glaciers in Gangotri Glacier System, Himalaya

Glacier-wide mass balances (MBs) of the Gangotri, Chaturangi, Raktavaran, Meru, and Gangotri Glacier System are reconstructed with a temperature-index (T-index) model using bias-corrected ERA5 data at a daily temporal resolution over 1979–2020. The model output is calibrated against available geodetic MB for Gangotri Glacier System and validated with satellite-derived snow line altitudes (SLAs) for Gangotri Glacier. Gangotri and Meru glaciers show mean mass wastage of –0.88 ± 0.31 m w. e. a‒ˡ (meter water equivalent per year) and ‒0.17 ± 0.29 m w. e. a‒ˡ, respectively whereas the mass budgets of fragmented tributary Chaturangi and Raktavaran glaciers are positive with the mean values of 0.49 ± 0.17 m w. e. a‒ˡ and 0.62 ± 0.15 m w. e. a‒ˡ, respectively over 1979–2020. Gangotri Glacier’s tongue is covered by thick debris having several supra-glacial lakes and ice cliffs (considered as melting hotspots); therefore, despite the presence of thick debris, we assume the melting over this area as of a clean glacier. The whole Gangotri Glacier System shows a moderate wastage of ‒0.27 ± 0.25 m w. e. a‒ˡ. The positive MBs of the Raktavaran and Chaturangi glaciers are due to their high area-elevation distribution and heavily debris-covered tongues. The positive MBs on these fragmented tributary glaciers are due to non-climatic topographic reasons and should not be misunderstood as climate change deniers or compared with Karakoram Anomaly. Modelled MBs are most sensitive to the threshold temperature for melt. The altitudinal MB sensitivities to all model parameters become negligible above 6,200 m a.s.l.

(including the Karakoram Range), the Himalaya implies the largest concentration of glaciers outside the polar regions (Nuimura et al., 2015). Emerging evidence of continuous, accelerated, and spatially heterogeneous rates of glacier mass loss in the Himalaya since the mid-nineteenth century raised the scientific community's attention to the regional differences in climatic conditions and the paucity of in situ observations in the Himalaya (Brun et al., 2017;Sakai and Fujita, 2017;Azam et al., 2018;Bolch et al., 2019;Maurer et al., 2019;Shean et al., 2020). In the west of the Himalaya, Karakoram glaciers have demonstrated an exceptional worldwide case of balanced or positive mass budgets, termed the "Karakoram Anomaly" (Hewitt, 2005;Kääb et al., 2012;Gardelle et al., 2013).
In situ glacier-wide mass balance (MB) monitoring of glaciers provides an understanding of the quick response of the glaciers to local meteorological conditions (Oerlemans, 2001). Several in situ glaciological studies are being conducted over the different basins of the Himalaya (Dobhal et al., 2013;Mandal et al., 2020;Wagnon et al., 2021;Stumm et al., 2021) but these observations are still sparser compared to the other glacierized regions of the world (Zemp et al., 2015;Azam et al., 2018). Due to the harsh climatic conditions, low oxygen level, and steep terrain of the Himalaya, the in situ MBs have only been observed for a short period and on a few small glaciers . However, as satellite missions and remote sensing methods progressed, it has become possible to work at a regional scale using geodetic method Shean et al., 2020), but the uncertainty associated with the sensors and inability to estimate the interannual and seasonal variations of glacier MB limit its applicability to understand the glacierclimate relationship .
Due to the scarcity of in situ MB observations in the Himalaya, modelling approaches have been used to understand the MB variability with climatic parameters . These models range from simple temperature-index (T-index) to complex surface energy balance (SEB) models (Hock and Holmgren, 2005;Pellicciotti et al., 2005;Srivastava and Azam, 2022a). The SEB model applications, required to understand the impact of different climate variables on glacier health, are limited in the Himalaya because of the need of extensive input data (Azam et al., 2014a;Mandal et al., 2022;Oulkar et al., 2022). On the other hand, the T-index models require minimal data, often only the temperature (Hock, 2003). The good performance of T-index models compared to SEB models is attributed to the fact that many components of energy balance, such as long wave radiation, sensible heat flux, are strongly correlated with temperature (Ohmura, 2001). Therefore, despite the simple computation of the T-index models, they have been applied widely to the Himalayan glaciers to reconstruct the long-term glacier MBs (Azam et al., 2014b;Shea et al., 2015;Litt et al., 2019;Pratap et al., 2019;Srivastava and Azam, 2022b).
Glacier MB reconstructions using the simple T-index models have been carried out in different parts of the Himalaya, especially on the glaciers where some in situ data is available . A short-term study on the Chandra Basin (8 glaciers from whole basin, area ranging from 2 to 78 km 2 ) showed the mean MB of -0.71 ± 0.34 m w. e. aˡ over -2009(Tawde et al., 2016. Another study modelled the mean MB of -0.68 m w. e. aˡ over 1985-2014 on four selected glaciers (Naradu, Shaune Garang, Gor-Garang and Gara glaciers) in the Baspa Basin (Gaddam et al., 2017). The MB of Naimona'nyi Glacier (14.4 km 2 ) was reconstructed as -0.40 ± 0.17 m w. e. aˡ over 1974-2014 (Zhao et al., 2016). Patsio Glacier (Bhaga Basin having area of 2.5 km 2 ) showed the modelled mean mass loss of -0.10 ± 0.10 m w. e. aˡ over the period 1993-2018 (Kumar et al., 2021). The longest reconstructed MB series is available from Chhota Shigri and Dokriani Bamak glaciers with moderate mass wastages of -0.12 ± 0.28 m w. e. aˡ and -0.09 ± 0.35 m w. e. aˡ over the last 7 decades   (Srivastava et al., 2022). These recent applications suggest that the T-index models can reconstruct the mass balances, especially in datascarce region like the Himalaya.
Further, the in situ observations, as well as the model applications, are often available from small glaciers, and the large glaciers have not been investigated in the Himalaya. The presence of ice cliffs and supra-glacial lakes on the debris-covered surface of the large Himalayan glaciers is another challenge that cannot be included in T-Index model due to extensive data requirement to compute melt over these surfaces.
For the present study, we have selected Gangotri Glacier System (Gangotri, Chaturangi, Raktavaran and Meru glaciers) and reconstructed its annual and seasonal MBs since 1979 using a simple T-Index model combined with an accumulation model. The Gangotri Glacier System is selected because it has one of the biggest glaciers (Gangotri Glacier) in India and has already been studied using the geodetic approaches (Bhattacharya et al., 2016;Bhushan et al., 2017). The major objectives of this study are 1) to test T-index model applicability for a large, highly-debris covered glacier with ice cliffs and supra-glacial lakes, 2) to model the longterm annual and seasonal MB series for all glaciers of the Gangotri Glacier System and 3) to understand the modelled MB sensitivity to input parameters and climate data.

Study area and climatic conditions 2.1 Study area
The Gangotri Glacier System (30.72°-31.02°N, 78.99°-79.29°E) is the largest glacier system in the Bhagirathi Basin, located in the Garhwal range of the central Himalaya in the Uttarakhand state of India. It is having four glaciers namely: Gangotri, Chaturangi, Raktavaran and Meru (Figure 1; GSI, 2009). The system originates from the Chaukhamba massif~7,000 m a.s.l, flows from SE to NW within the granitic terrain, and terminates at 4,000 m a.s.l. with a total glacierized area of~252 km 2 ( Figure 1). Gangotri Glacier System shapefile was extracted from the GAMDAM Glacier Inventory, in which the glacier outlines were delineated manually using high resolution Google earth images, DEM, and Landsat ETM data from 1999 to 2003. (Nuimura et al., 2015). Gangotri Glacier System also has its religious importance, and its snout is called Gaumukh (mouth of a cow). The Bhagirathi River flows from the Gaumukh to Devprayag, where it joins the Alaknanda River to form the Ganga River.
The Gangotri Glacier System is surrounded by mystic peaks (ranging from 6,000-7,000 m a.s.l.). The main trunk of Gangotri Glacier is~32 km long and 1-3 km wide, with Meru Glacier on the left and Chaturangi and Raktavaran glaciers on the right (Figure 1). Historical evidences indicate that Meru, Chaturangi and Raktavaran glaciers were tributary glaciers of the main Gangotri Glacier and got fragmented in the past (GSI, 2004). The topographical details of all glaciers in the Gangotri Glacier System are given in Table 1. All four glaciers of the Gangotri Glacier System are heavily debris-covered ( Figure 1) and have around 20% debris cover area, found increasing with an increasing rate over recent decades (Bhattacharya et al., 2016). Gangotri Glacier has the highest debris cover (24% of its total area) and has numerous ice cliffs and glacial lakes on its lower ablation area (inset of Figure 1) that are expected to act as melting hotspots hence amplify the local melting (Brun et al., 2018;King et al., 2019;Shugar et al., 2020). Debris covers, ice cliffs and supraglacial lakes were manually delineated on the Google earth image of 6 October 2019, near the end of the melting season for clear visibility of all these features. The climatic condition of the Gangotri Glacier System is influenced regionally by the Indian summer Monsoon (ISM) and Indian winter Monsoon (IWM) (Dimri et al., 2016;Kotlia et al., 2018).

Data, bias correction and climatic conditions
Daily reanalysis temperature and precipitation data available at 0.25°X 0.25°resolution from ERA5 (https://cds.climate. copernicus.eu) were used to determine the annual and seasonal MBs of all glaciers in the Gangotri Glacier System since 1979. The ERA5 reanalysis temperature and precipitation data were downloaded at the nearest grid point (~8 km) to the Bhojbasa Base Camp (3,800 m a.s.l.) (Figure 1).
The mean monthly temperature data from May 2006 to April 2007 from an Automatic Weather Station (AWS) at Bhojbasa Base Camp (Agrawal et al., 2018) was used for bias correction of raw ERA5 temperature data. Raw ERA5 mean monthly temperature exhibited a high correlation with AWS temperature (Figure 2A). In situ precipitation data was available only for the summer months (May-October) over -2003(Singh et al., 2006 and utilised for bias correction of the raw ERA5 precipitation data. ERA5 precipitation data also showed high correlation but high overestimation ( Figure 2B) with the available mean monthly precipitation data. The multiplicative monthly bias-correction factors, developed from the linear regression equation between the mean monthly ERA5 and in situ temperature and precipitation data, were applied to bias-correct the daily data. The multiplicative factors were 1.12 and 0.42 for the temperature and precipitation data, respectively (Figures 2A,B).
The long-term in situ meteorological data is not available hence in this study we used the bias-corrected ERA5 temperature and precipitation for the period 1979-2020 to understand the local temperature and precipitation distribution over the year. The bias-corrected ERA5 annual mean temperature at the Bhojbasa Base Camp was 1.89°C, with the highest and lowest annual mean temperatures of 3.26°C and 0.62°C for 2016 and 1997, respectively ( Figure 3A). The mean temperatures in the summer (May-October) and winter (November-April) months were 6.89°C and -3.21°C, respectively. July was the hottest month with a mean temperature of 10.1°C and January was the coldest with a minimum temperature of -6.74°C over 1979-2020 ( Figure 3B). The bias-corrected ERA5 mean annual precipitation was 492 mm, with 54% occurring during the summer season and 46% during the winter season, indicating that Gangotri Glacier System receives nearly equal precipitations in both the seasons ( Figure 3B).
Gangotri Glacier System is located on the leeward side of the mountain range. The mean annual precipitation at Bhojbasa Base Camp (492 mm) is only 30% of that of Dokriani Bamak Glacier Base Camp (1,616 mm), a glacier that is located on the orographic front in the same range around 30 km SW of Gangotri Glacier (Azam and Srivastava, 2020). Because highaltitude orographic forcing influences the amount of  Frontiers in Earth Science frontiersin.org 04 precipitation (Bhaskaran et al., 1996;Dimri et al., 2013), the ISM reach over the Gangotri Glacier System is quite poor. Dokriani Bamak Glacier receives its 77% of annual precipitation from the summer months and is thought to be a summer-accumulation type glacier (Srivastava and Azam, 2022a), while Gangotri Glacier System receives an almost equal amount of summer and winter precipitations ( Table 1) that makes it difficult to say whether it is summer-or winter-accumulation type glacier system. This comparison clearly demonstrates that different types of glaciers can co-exist in the same range due to the complex topography that gives rise to the local microclimate.

Mass balance model
The annual and seasonal MBs of all glaciers in the Gangotri Glacier system are reconstructed by applying a T-index model together with an accumulation model that has been applied in several previous studies for MB reconstructions in the Himalaya (Azam et al., 2014b;Srivastava and Azam, 2022a). The model is forced with the bias-corrected ERA5 daily precipitation and temperature data since 1979. The model runs at a daily time scale and estimates the daily accumulation and ablation for each altitudinal range of 50-m for a full hydrological year, from November 1 to October 31.
The daily precipitation and temperature are extrapolated to each 50-m altitudinal range by applying the altitudinal precipitation gradient (P g ) and temperature lapse rates (LRs), respectively.
Where P is the calculated precipitation at different altitudinal bands, P i is the precipitation at base camp, ΔH is the altitudinal difference between Bhojbasa Base Camp and the corresponding altitudinal band (m) and P g is the precipitation gradient in % precipitation change per 1,000 m altitude. P g has been used to extrapolate the precipitation at a daily scale because our model runs at a daily scale.
Where T is the calculated temperature at different altitudinal bands, T i is the temperature at Bhojbasa Base Camp, ΔH is the altitudinal difference between base camp and corresponding altitudinal band (m) and LR is the monthly lapse rate. The daily accumulation (A) is estimated by using the threshold temperature for snow/rain (T P ) and the daily ablation (M) is estimated by using the threshold temperature for melt (T M ). Accumulation of snow is taken place when T < T P , and if T > T P all the snow is melted.
The daily accumulation A (mm w. e. d −1 ) at each altitudinal range is computed by: Where P and T are daily precipitation (mm) and temperature (°C) respectively extrapolated at each altitudinal range and T P is the threshold temperature (°C) for snow-rain. Ablation (M) is taken place when T > T M , otherwise it is computed to be zero. During melting, first, the accumulated snow is melted out and then ice is started melting depending on the surface condition whether it is debris cover or clean ice.
The daily ablation M (mm w. e. d −1 ) at each altitudinal range is computed by: where, DDF S/I/D denotes the degree-day factor (mm d −1°C−1 ) for snow (S), ice (I) and debris-covered ice (D) surfaces, and T is extrapolated daily air temperature (°C) at each altitudinal range.

FIGURE 3 (A) Bias-corrected annual mean temperature (blue dots) and annual precipitation sums (green bars) over 1979-2020 and (B)
bias-corrected mean monthly precipitation and temperature patterns at Bhojbasa Base Camp (black dots represent the temperature, blue-green bars represent the winter precipitation, and the grey bars represent the summer precipitation) and the pie chart inset shows the percent of seasonal precipitation contribution.
Frontiers in Earth Science frontiersin.org 05 In addition, DDF I was also employed for ablation on the thick debris-covered surface of Gangotri Glacier which has numerous ice cliffs and supra-glacial lakes (Section 5.2).
Daily mean altitudinal MB (b i ) at each 50-m altitudinal range is calculated using daily M and daily A (from Eqs. 3, 4): Glacier-wide MB (B a ) is calculated as: where b i is the altitudinal MB at each altitude, a i is the area for each altitude, and A is the total glacier area.

Model parameters
The temperature is the most important variable controlling the distribution of snow/rain on the surface and the melting of snow and ice surface in T-index glacier MB modelling (Shea et al., 2015). The in situ DDF for snow, ice, and debris-cover are not available for Gangotri Glacier System and were taken as 6.1, 7.7 and 4.8 mm d −1°C−1 , respectively, which were calculated (Azam and Srivastava, 2020) based on the previous in situ observations from different studies on the Dokriani Bamak Glacier located~30 km west of Gangotri Glacier System (Singh et al., 2000;Pratap et al., 2015). The T P is taken as 0.7°C from (Jennings et al., 2018), which corresponds to 70-80% relative humidity ranges for Gangotri Glacier System. The temperature is extrapolated at different altitudinal ranges using the monthly LRs developed on the Dokriani Bamak Glacier catchment (Azam and Srivastava, 2020).
Monthly LRs vary with the highest mean monthly LR (6.94°C km −1 ) in May and the lowest mean monthly LR (5.42°C km −1 ) in August. T M and P g are used to calibrate the model (Section 3.3). All the model parameters are listed in Table 2.

Model calibration
The precipitation distribution in the complex Himalayan terrain is poorly known because of very limited in situ observations (Azam et al., 2021), hence the suitable selection of P g for precipitation distribution on the glaciers is always challenging in the mountainous region Bolch et al., 2019). Another critical parameter is T M , for which the MB models are very sensitive (Engelhardt et al., 2017;Azam and Srivastava, 2020). Further, studies suggest that sometimes melting does not even happen at an air temperature of >0°C (Kuhn, 1987;Hock, 2003) however, some T-index model studies suggest that the calibrated T M can be negative at the daily time step (van den Broeke et al., 2010;Engelhardt et al., 2017;Azam and Srivastava, 2020). So, P g and T M are chosen as the calibrating parameters and varied to fit the model output with the observed geodetic MB data. The T M and P g were altered between -5 and 5°C, and 0% km −1 and 100% km −1 , respectively. The Monte Carlo simulations (10,000 runs) were performed for the model calibration while minimizing the root mean square errors (RMSE) between the modelled and available geodetic MB (-0.29 ± 0.19 m w. e. aˡ) over 2006-2014 for Gangotri Glacier System (Bhattacharya et al., 2016). The geodetic MBs for different glaciers in Gangotri Glacier System are not available separately hence the geodetic MB available for the whole Gangotri Glacier System has been used for model calibration. T M of 0°C and P g of 46% km −1 showed the least RMSE of 0.08, corresponding to a mass wastage of -0.27 ± 0.25 m w. e. aˡ over 2006-2014 for the Gangotri Glacier System. Same calibrated parameters were applied separately on Gangotri Glacier and other fragmented tributary glaciers of the Gangotri Glacier System to simulate the corresponding annual and seasonal MBs.

Model validation
The model is validated against the manually delineated snow line altitudes (SLAs) derived from the Landsat satellite images with the modelled SLAs. The modelled SLAs are independently derived with the calibrated parameters at the 5-m altitudinal range. The uncertainty in modelled SLAs was assumed to be equal to 10 m as the model ran at a 5-m altitudinal range. Total of 19 SLAs are derived from Landsat satellite images on cloud-free days from 1989 to 2019. The modelled and satellite-derived SLAs showed a good agreement but high RMSE (r = 0.69 and RMSE = 208 m, Figure 4). The possible reasons for mismatch could be due to an error in the manual delineation of SLAs caused by the presence of shadows in the satellite images as well as the coarser resolution of the images. Additionally, the assumption of a single SLA for the entire glacier in the model estimation could be the cause of the larger RMSE, as SLA may have several altitudes in different parts of the glacier due to various aspects and local topography (Racoviteanu et al., 2019).
The uncertainty in SLA is estimated using the similar method as followed in (Racoviteanu et al., 2019), in which the accuracy of the digital elevation model (DEM) and buffer size are used for the partition of surface and SLA extraction, respectively and the uncertainty in the snow and ice area estimates are also considered. We used Cartosat-1 DEM for surface partition and Landsat satellite images for snow and ice area estimation (Racoviteanu et al., 2019). The SLA uncertainty as ± 18.03 m was calculated with the formula: Where ε DEM is the vertical error of the Cartosat-1 DEM as ± 9.6 m (Singh et al., 2013) and ε SNOWLINE is half the pixel size of the Landsat image as ± 15 m (Racoviteanu et al., 2019).

Uncertainty estimation
Each model parameter is altered one by one, within a 10% range of its calibrated value (except for parameters that are constrained using in situ data, and uncertainty range is known and used) to estimate the parametric uncertainty on the model outputs, while the other parameters are maintained constant (Anslow et al., 2008;Heynen et al., 2013;Ragettli et al., 2013). The plausible ranges of model parameters are given in Table 2. The overall model uncertainties for the modelled MBs were estimated using error propagation law and found as ± 0.31, ±0.17, ±0.15, ±0.29, ±0.25 m w. e. aˡ for Gangotri Glacier, Chaturangi Glacier, Raktavaran Glacier, Meru Glacier and Gangotri Glacier System, respectively.
The areal extent and hypsometry of the study area have been kept fixed for modelling the MBs in the present study. The areal shrinkage on Gangotri Glacier between 1962 and 2006 was only 6% of its area in 1962 which accounts for a very low rate of 0.14% aˡ (Negi et al., 2012) so the uncertainty due to the fixed-area assumption can be assumed negligible. In line, (Bhattacharya et al., 2016), also estimated a very limited area loss of 0.33% (0.01% aˡ) between 1965 and 2015 for Gangotri Glacier. The changes in surface elevation due to the temperature change also induce some MB uncertainty, which is negligible compared to the overall uncertainty in MB  hence it was ignored.

Modelled annual and seasonal glacierwide mass balance
The annual and seasonal MBs were computed separately for Gangotri, Chaturangi, Raktavaran and Meru glaciers, and for the whole Gangotri Glacier System over 1979-2020 ( Figure 5). Frontiers in Earth Science frontiersin.org 07 summer MB was estimated between 1st May to 31st October while the winter MB was estimated between 1st November and 30th April for each hydrological year.
The Gangotri Glacier continuously lost mass for all years over 1979-2020, while the Meru Glacier experienced a few balanced/slightly positive MB years (1982, 1989, 1992, and 2005) (Figures 5A,D). Conversely, the Chaturangi and Raktavaran glaciers showed positive MBs continuously for all years ( Figures 5B,C). As a result, the Gangotri Glacier System showed a continuous but moderate mass loss from 1979 to 2020, except for 1989 when the MB was slightly positive ( Figure 5E). The exceptional positive MBs on Chaturangi and Raktavaran glaciers, despite a general wastage in the Himalaya, have been discussed in Section 5.1.

Equilibrium line altitude and accumulation area ratio
The mean modelled ELAs were found to be 5358, 5350, 5347, 5358, and 5356 m a.s.l. for Gangotri, Chaturangi, Raktavaran, Meru glaciers and Gangotri Glacier System, respectively over 1979-2020 ( Figure 5). The nearly equal ELAs and their similar interannual variability on all glaciers ( Figure 5) were due to the model calibration process where the whole Gangotri Glacier System is calibrated against the geodetic MB (Section 3.3) that resulted in the same distribution of the temperature and precipitation fields on all glaciers. However, the mean annual AARs, which depend on the ratio of ablation and accumulation areas, were quite different, with the values of 45%, 81%, 86%, 54%, and 61% corresponding to an annual mean MBs of -0.88, 0.49, 0.62, -0.17, and -0.27 m w. e. a −1 for Gangotri, Chaturangi, Raktavaran, Meru glaciers, and Gangotri Glacier System, respectively.
Surprisingly, the modelled ELAs on Chaturangi and Raktavaran glaciers are slightly below the debris cover (Figure 7). These glaciers are expected to receive comparatively more solar radiation because of their southwest aspect, resulting in more melting hence higher ELAs than Gangotri and Meru glaciers. As our model provided uniform precipitation and temperature fields hence the accumulation and ablation on all the glaciers, the modelled ELAs on Chaturangi and Raktavaran glaciers are expected to be lower than the actual ELAs or SLAs at the end of ablation season. We delineated the SLAs (see details in Section 3.4) for Chaturangi and Raktavaran glaciers and found that the SLAs are higher than both the debris cover as well as the modelled ELAs. For Chaturangi and Raktavaran glaciers, the mean SLA (delineated for 2013, 2015, 2019 and 2020) and the mean modelled ELA (for same years) differed by 101 m and 112 m, respectively. Modelled ELAs and SLAs on Gangotri Glacier (Section 3.4) as well as some previous studies, comparing modelled ELAs with field-based ELAs, have shown more than 200 m difference among modelled ELAs, delineated SLAs and in situ ELAs (Wang et al., 2017;Azam and Srivastava., 2020). Though, the difference in our modelled ELAs and delineated-SLAs is only around 100 m, yet it is clear that the modelled ELAs (below debris cover) are underestimated; therefore, it is advised to be cautious when this data is used.
The relationships between the modelled MBs, ELAs, and AARs for all glaciers were established individually using 41 years of modelled values ( Figure 6). The steady-state ELAs (ELA 0 ) were also estimated as 5177, 5529, 5838, 5374, and 5267 m a.s.l. for Gangotri, Chaturangi, Raktavaran, Meru glaciers, and the Gangotri Glacier System, respectively using regression fits for zero MB over 1979-2020. The ELA 0 of 5267 m a.s.l. for the Gangotri Glacier System was in good agreement with the ELA 0 of 5280 m a.s.l. for the Dokriani Bamak Glacier (Azam and Srivastava, 2020). Similarly, the steady-state AARs (AAR 0 ) were 60, 65, 70, 58, and 66% for Gangotri, Chaturangi, Raktavaran, Meru and Gangotri Glacier System, respectively. The AAR 0 values for our study area were well-matched to the AAR 0 ranges (43%-73%) compiled by Pratap et al. (2016) for various Himalayan glaciers.

Positive MBs of fragmented tributary Chaturangi and Raktavaran glaciers
The modelled MBs of Chaturangi and Raktavaran glaciers were positive with mean values of 0.49 ± 0.17 and 0.62 ± 0.15 m w. e. aˡ, respectively over 1979-2020 ( Figure 5). Further, their summer mean MBs were also positive with occasional negative summer MBs ( Figure 5). This seems surprising when the Himalayan glaciers have been losing mass over the last 5-6 decades Bolch et al., 2019). Chaturangi and Raktavaran are quite big (75 km 2 and 30 km 2 , respectively) and high-altitude (4,450-6,900 m a.s.l., and 4,600-6,600 m a.s.l, Frontiers in Earth Science frontiersin.org 09 respectively) glaciers. The area-elevation distribution is also very high on both the glaciers hence most of their glacierized area is above ELA (Chaturangi = 81% and Raktavaran = 86% accumulation area, respectively) -giving rise to a very small ablation area exposed to positive summer temperatures ( Figures  7B,C). Further, ablation areas are covered by thick debris (~15% area covering these glaciers from snouts to slightly above ELA) that protects these glaciers from higher melting. This exceptionally high area-elevation distribution coupled with debris-covered tongues resulted in very limited wastage, easily compensated by the large summer accumulation ( Figures 7B, 5). Hence, Chaturangi and Raktavaran glaciers showed positive summer mean MBs that further produced positive annual mean MBs. Gangotri and Meru glaciers have significant ablation areas below ELA (45% and 54%, respectively; Figures 7A,D) hence their modelled mean annual MBs were negative. The whole Gangotri Glacier System also has a significant ablation area below ELA (39%, Figure 7E) that produced negative summer MBs resulting in the moderate annual mass loss (Section 4.1).
The modelled cumulative MBs for the Chaturangi and Raktavaran glaciers were 20.2 ± 1.12 and 25.5 ± 0.94 m w. e, respectively over 1979-2020. These continued cumulative positive MBs must have resulted in terminus advancement through higher ice flux transfer to the terminus, and as a result, both the glaciers must have advanced and rereconnected to the main Gangotri Glacier trunk, but these are still fragmented from the Gangotri Glacier (Figure 1). Previous literature suggests that Chaturangi, Raktavaran and Meru glaciers were connected to the main Gangotri Glacier in the past and got fragmented from the main glacier (GSI, 2004). Chaturangi Glacier was an independent glacier during the 1960s, then it advanced and attached to the Gangotri Glacier in 1971 (Srivastava, 2012), and then again detached sometime around the early 1990s (Bisht et al., 2020), while Raktavaran seems to be an independent glacier even during the last Little Ice Age (GSI, 2004). The information about the Meru Glacier detachment is not available in the literature. Chaturangi Glacier is having retreat and advance periods however our model-not involving glacier dynamics-cannot provide insights on this. During the field surveys in 2010 (Pottakkal et al., 2014), it was observed that the Chaturangi Glacier snout lies in a steep and deep gorge bounded by unstable steep valley walls. In such a situation, the continuous rockfalls and gushing water probably provide some mechanical ablation to the snout. However, detailed research about Chaturangi Glacier's snout advance/retreat (Srivastava, 2012) is needed using some sophisticated models including the glacier dynamics and the local topography in the modelling scheme.
We hypothesize that the cumulative positive MBs on Chaturangi and Raktavaran glaciers lead to the glacier advancement and a possible reconnection to the main Gangotri Glacier trunk. Therefore, to determine the impact of the existing unglacierized area between these glaciers and the Gangotri Glacier trunk on modelled MBs, we re-ran the model to reconstruct the MBs assuming that both glaciers are linked with the Gangotri Glacier trunk. This assumption gives an additional ablation area of 0.41 km 2 on Chaturangi Glacier and 0.59 km 2 on Raktavaran Glacier. With these additional ablation areas, the modelled mean MBs of 0.47 and 0.55 m w. e. aˡ were obtained on Chaturangi and Raktavaran glaciers, respectively over 1979-2020. The minor differences of 0.02 m w. e. aˡ and

Frontiers in Earth Science
frontiersin.org 0.07 m w. e. aˡ on Chaturangi and Raktavaran glaciers, respectively show that the advancement has almost no effect on the modelled MBs on these glaciers.

Assumption of DDF for thick debris having ice cliffs and supra-glacial lakes and suitability of T-index model
Supra-glacial debris surface is ubiquitous on the Himalayan glaciers that strongly impacts glaciers' response (Scherler et al., 2011;Banerjee and Shankar, 2013). While the thin debris cover of a few millimeters accelerates the melting, the thick debris cover of more than a few centimeters can protect the glacier from higher melting (Shah et al., 2019). It has also been observed that the thick debris-covered glaciers may have ice cliffs and supra-glacial lakes that result in enhanced, localized melting; hence ice cliffs and supra-glacial lakes are considered as the melting hotspots Miles et al., 2016). Further, several studies focusing on large-scale glacierized regions have revealed that both debris-covered and clean glaciers experience similar mass wastage (Kääb et al., 2012;Brun et al., 2019), a phenomenon that has been coined as the "debris-cover anomaly" (Pellicciotti et al., 2015;Vincent et al., 2016). This is due to the presence of hotspots on the debris-covered area of some large glaciers that give rise to higher localized melting Miles et al., 2016) that compensates for the lower melting on the debris-covered region; therefore, near-similar mass wastage on debris-covered and clean glaciers at regional scale (Kääb et al., 2012).
In Gangotri Glacier System, ice cliffs and supra-glacial lakes are evident only on the lower ablation area (4,000-4,750 m a.s.l.) of the Gangotri Glacier (Figure 1), while the lower ablation area on other glaciers (Meru, Chaturangi and Raktavaran) has debris cover with some occasional hotspots. The estimated area for ice cliffs and supra-glacial lakes on the Gangotri Glacier was 1.09% (with the maximum contribution of 3.9% from the altitude 4,200-4,250 m a.s.l) and 0.5% (with the maximum contribution of 1.6% from the altitude 4,350-4,400 m a.s.l) of total debris cover area (Figure 9), but these have been found to enhance the melting significantly Miles et al., 2016;Brun et al., 2018). Therefore, we used DDF D to compute the MBs of the debris-covered area on all glaciers except the Gangotri Glacier, where we applied DDF I to compute the MB over the ablation area with prominent ice cliffs and supra-glacial lakes,

FIGURE 7
The 50-m hypsometry for four selected glaciers (A-D) and the Gangotri Glacier System (E). Light grey bars represent the glacier area, and the dark grey bars represent the debris-covered area on each glacier. The outer pie charts show the percentage of debris-covered and clean-ice glacier surface area, and the inner pie charts show the areas below and above the equilibrium line altitudes (ELAs). The black dashed line is the mean ELA of the Gangotri Glacier System.

Frontiers in Earth Science
frontiersin.org 11 assuming that the accelerated melt due to hotspots would compensate for the lower melt over the debris-covered area on Gangotri Glacier tongue (Kääb et al., 2012). DDF I is applied over 4,000-4,750 m a.s.l. because this range has hotspots while the higher range 4,800-5350 m a.s.l. has only debris cover (hence DDF D was applied).
To check our assumption of DDF I application for melt generation over the lower debris-covered area on the Gangotri Glacier, we extracted the glacier thinning rates (dh/dt) from Shean et al. (2020) for all four glaciers of our study area as well as for some other glaciers (clean and debris cover) (Figure 8). The thinning patterns on Meru, Chaturangi and Raktavaran glaciers show that the thinning is subdued on these glaciers due to the presence of debris cover over the ablation area ( Figure 8B). In agreement, Bhagirathi Kharak and Satopanth glaciers in the same region, which too originate from Chowkhamba Massif (flowing SE, opposite to Gangotri Glacier), also showed subdued thinning over the debris-covered ablation area ( Figure 8B). On the other hand, Gangotri Glacier witnessed continued thinning over the lower debris-covered area ( Figure 8A). We selected a random clean glacier (inset in Figure 8A) from the same region (60 km NW from Gangotri Glacier System) and a large clean glacier of Durung Drung (68 km 2 ) from the Zanskar Range and compared their thinning patterns with that of Gangotri Glacier. The higher thinning over lower elevations on Gangotri Glacier is analogous to the thinning patterns of a randomly selected clean glacier and Durung Drung Glacier ( Figure 8A). This reveals that the lower debris-covered ablation area (4,000-4,750 m a.s.l.) of Gangotri Glacier is wasting its mass like clean glaciers and supports the application of DDF I for computation of melt over the debriscovered area with prominent melting hotspots. Another geodetic study also suggested that the lower ablation area of Gangotri Glacier showed considerable thickness loss despite thick debris cover on its surface, while other glaciers in the system showed subdued thickness losses over the debris-covered area (Bhattacharya et al., 2016).
Though the comparison of thinning patterns ( Figure 8A) supports the use of DDF I for melt computation over the lower debris-covered area on Gangotri Glacier, it is an indirect support and can be questioned. Some studies directly compared the modelled altitudinal MBs with the in situ observations (Azam et al., 2014b;Azam et al., 2019), but such a direct comparison is not possible with geodetic thickness changes as later is the net result of glacier wastage and emergence velocity (Cuffey and Paterson, 2010), that is not accounted in our model. However, our analysis is qualitative, yet it supports the application of T-index model for the MB reconstruction on a large debriscovered glacier with supra-glacial lakes and ice cliffs.
To investigate the impact of our DDF I application for melt generation over the debris-covered area of Gangotri Glacier, we also ran the model with DDF D for melt computation over this area. The mean MB with DDF D was computed as -0.60 m w. e. aˡ against the mean MB of -0.88 m w. e. aˡ with DDF I for Gangotri Glacier over 1979-2020. The additional ablation contribution, according to our DDF I assumption, is 32%, which is the combined effect of both ice cliffs and supra-glacial lakes. Our result is in good agreement with the Brun et al. (2018) study, which found a 23% ablation contribution from ice cliffs on Changri Nup Glacier, central Himalaya in Nepal. Figure 9 shows the modelled mean altitudinal MB distribution for Gangotri Glacier System over 1979-2020 along with the altitude-wise percentage areal coverage of ice cliffs as well as the supraglacial lakes on the Gangotri Glacier trunk. The altitudinal MB values vary from -8.22 to 1.21 m w. e. for Gangotri, -2.80 to 1.20 m w. e. for Chaturangi, -2.19 to 1.13 m w. e. for Raktavaran, and -3.38 to 1.11 m w. e. for Meru glaciers, and -8.22 to 1.21 m w. e. for Gangotri Glacier System. The highest altitudinal mass loss was observed over the lower ablation area of the main trunk of the Gangotri Glacier (Figure 9). Conversely, the fragmented tributary (Chaturangi, Raktavaran and Meru) glaciers showed relatively less altitudinal mass loss over their lower ablation areas (Figure 9). This is because the lower ablation areas of other glaciers are debris-covered and located at relatively higher

Frontiers in Earth Science
frontiersin.org altitudes that protect from higher melting on these glaciers (Section 5.1). We have also estimated MB gradients using the 50-m modelled altitudinal MBs by fitting linear regressions individually for each zone based on the glacier surface morphology (debris or clean ice) and location of ELA (as a demarcation for ablation and accumulation MB gradients) (Figure 9). The MB gradients were estimated as 0.74, 0.36, and 0.06 m w. e. per 100 m for the debris-covered ablation area having hotspots, debris-covered ablation area and accumulation area, respectively for Gangotri Glacier ( Figure 9A). Despite considering debris-covered surface with hotspots as clean ice surface (detail in Section 5.2), the Gangotri Glacier MB gradient (0.74 m w. e. per 100 m) was slightly lower than the Dokriani Bamak Glacier MB gradient for the clean ablation [0.91 m w. e. per 100 m, (Azam and Srivastava, 2020)]. The debris cover ablation area of the Chaturangi and Raktavaran glaciers showed very low MB gradients of 0.30 and 0.29 m w. e. per 100 m, respectively because of the substantial debris cover and high altitudes that collectively give rise to less negative altitudinal MBs (Figures 9B,C). For Meru Glacier, MB gradients were calculated for debris-covered ablation area, clean ablation area and accumulation area, and found as 0.32, 0.39, and 0.08 m w. e. per 100 m, respectively ( Figure 9D).
For the entire Gangotri Glacier System, the resulted MB gradient for the ablation surface was estimated to be 0.61 m w. e. per 100 m because of its high altitudinal range and high debris cover, which is slightly lower than the MB gradient value for the debris surface of nearby Dokriani Bamak Glacier (0.70 m w. e. per 100 m, Azam and Srivastava, 2020) ( Figure 9E). Very low MB gradients in accumulation zones of all glaciers were found because of little or no melting there, and it is attributable solely to snow accumulation for all the glaciers, which is similar on all glaciers due to calibration process (Section 3.3).

Modelled mass balance sensitivity
The sensitivity of the modelled MB to each parameter was determined by rerunning the model with a uniform change in each model parameter (Oerlemans et al., 1998;Braithwaite and Zhang, 2000;Azam et al., 2014b), while remaining all other model parameters unaltered. These sensitivities were estimated by (Oerlemans et al., 1998): Where B a is the average annual and seasonal MB for the period 1979-2020 and µ H and µ L are the highest and lowest values of parameters, respectively. Table 2 provides the computed sensitivities for each model parameter. The sensitivity tests were carried out independently for each glacier. The sensitivity of MB to LR was determined by altering the monthly LR values with the standard deviation of monthly LRs. For the MB sensitivity tests, T P , T M and DDFs were modulated by 1°C, 1°C, and 1 mm d −1°C−1 , respectively, and P g was varied by 10%. In this analysis, sensitivity tests were also carried out for input precipitation and temperature, with variations of 10% and 1°C, respectively.
The MBs are most sensitive to T M with sensitivities of 0.81, 0.25, 0.19, 0.55, and 0.56 m w. e. aˡ for Gangotri, Chaturangi, Raktavaran, Meru glaciers, and Gangotri Glacier System, respectively ( Table 2). The similar sensitivities were also found on the other Himalayan glaciers: MB was most sensitive to T M with the sensitivities of 0.44 m w. e. aˡ for the Chhota Shigri Glacier and 0.77 m w. e. aˡ for the Dokriani Bamak Glacier Azam and Srivastava, 2020). Followed by T M , the modelled MB was most sensitive to the LRs with sensitivities of 0.29, 0.17, 0.14, 0.28, and 0.23 m w. e. aˡ for Gangotri, Chaturangi, Raktavaran, Meru glaciers, and Gangotri Glacier System, respectively. That is in line with Heynen and others, (2013), which stated that the model is quite sensitive to the LR for the large altitudinal-ranged glaciers. The model was found to be less sensitive to other model parameters (Table 2).
Figures 10A-E shows the relative sensitivity of modelled MB for each model parameter compared to T M sensitivity, kept as 100%. For the Gangotri Glacier and the whole Gangotri Glacier System, the modelled MBs were least sensitive to DDF D due to the assumption of a debris surface with hotspots as clean ice (Section 5.2), whereas for the Chaturangi and Raktavaran glaciers, the modelled MBs were least sensitive to DDF I due to the large coverage of debris cover on the ablation surfaces ( Figure 10).
We also plotted the altitudinal MB sensitivities to all the model parameters ( Figures 10F-J). The MB sensitivity to the most sensitive parameter, T M decreases with altitude, and becomes almost insensitive above 6,200 m a.s.l. for all the glaciers. This is because at lower altitudes the air temperature is usually higher than T M and there are always some positive temperatures up to 6,200 m a.s.l. Model sensitivity to T M changes abruptly around 4,800-5000 m a.s.l. on all glaciers because of sudden changes in surface conditions (debris-cover to clean ice or vice versa) ( Figure 10I). All glaciers' sensitivity to LR increases slightly with altitude in the ablation area and decreases with altitude in the accumulation area. The sensitivity to LR varies insignificantly in the debris-covered ablation area (Figures 10G-I) due to the debris cover's protection (hence less melting and sensitivity); however, at higher altitudes debriscovered area percentage gradually decreases (Figure 7) and thus sensitivity to LR gradually increases (Figures 10F,J). Above this transition zone, where debris-covered and cleanice areas coexist at the same elevation band, the LR sensitivity gradually decreases with air temperature and becomes zero at around 6,200 m a.s.l.
We assumed the lower debris-covered area with hotspots on Gangotri Glacier to be similar to clean ice (Section 5.2); therefore, the lower debris-covered area on Gangotri Glacier (as well as Gangotri Glacier System) shows high sensitivity to DDF I instead of DDF D , that show higher sensitivities on rest of the debris-covered glaciers ( Figure 10). The altitudinal MB sensitivity to DDF S and T P on all glaciers is negligible. The altitudinal sensitivity analysis clearly indicates that practically all the model parameters become insensitive at higher altitudes (6,200 m a.s.l.), and there is no effect on glacier health.
We also ran the model with 10% and 1°C variation in input precipitation and temperature data, respectively to understand the relative control of precipitation and temperature on the studied glaciers. The MB sensitivities were 0.07, 0.09, 0.09, 0.07, and 0.08 m w. e. aˡ to a 10% increase in precipitation, while the sensitivities were computed to be -0.54, -0.23, -0.18, -0.43, and -0.40 m w. e. aˡ to 1°C increase in temperature for Gangotri, Chaturanga, Raktavaran, Meru glaciers, and Gangotri Glacier System, respectively. We ran the model numerous times to estimate the additional precipitation required to compensate for a 1°C increase in temperature. The results showed that 83%, 32%, 24%, 66%, and 59% increase in precipitation is required to compensate 1°C rise in temperature on the Gangotri, Chaturangi, Raktavaran, Meru glaciers, and Gangotri Glacier System, respectively. The amount of precipitation increase by 59% to compensate for a 1°C temperature rise on the Gangotri Glacier System is in good agreement with the nearby Dokriani Bamak Glacier, where a 49% increase in precipitation was suggested to compensate for a 1°C temperature rise (Azam and Srivastava, 2020). Our results for Chaturangi and Raktavaran glaciers showed good agreement with Braithwaite et al. (2003), who reported a 30-40% increase in precipitation to compensate for the effect of a 1°C temperature rise based on 61 glaciers and ice caps from different regions of the world. The significant precipitation requirements for the Gangotri and Meru glaciers were due to a considerable number of days with temperatures above the T M , leading to more melting.
Frontiers in Earth Science frontiersin.org

Comparison with other studies
A few studies estimated the MBs of the Gangotri Glacier or Gangotri Glacier System using geodetic and modelling approaches. Our modelled mean MB on the Gangotri Glacier System was -0.27 ± 0.25 m w. e. aˡ over 1979-2020 while Bhattacharya et al. (2016) estimated a mean mass wastage of -0.19 ± 0.12 m w. e. aˡ over 1968-2014 using the geodetic approach. Our model was calibrated using mean MB data from the same study over 2006-2014 (Section 3.3). The slightly higher mass wastage in our study may be attributed to different time periods; nevertheless, the modelled mean MB is in good agreement considering the uncertainty bounds in MBs in both the studies. Another geodetic study estimated a mean mass wastage of -0.55 ± 0.42 m w. e. aˡ on Gangotri Glacier over 1999-2014 (Bhushan et al., 2017). Over the same period, our model showed a higher mean mass wastage of -0.89 ± 0.31 m w. e. aˡ on Gangotri Glacier. However, our modelled MBs on Gangotri Glacier are in better agreement with (Agrawal et al., 2018) who estimated the mean MB of -0.98 m w. e. aˡ on Gangotri Glacier using a simplified energy balance model over 1985-2005 while our study estimated a mean wastage of -0.84 ± 0.30 m w. e. aˡ over the same period.
The same study also estimated the mean MB as -0.92 m w. e. aˡ using an ice-flow velocity method over 2001-2014 that is in good agreement with the mean MB of -0.87 ± 0.30 m w. e. aˡ estimated in our study over the same period.
Previous studies estimated a much higher wastage when considering only Gangotri Glacier (Bhushan et al., 2017;Agrawal et al., 2018) compared to Gangotri Glacier System (Bhattacharya et al., 2016) and the same has been observed in our study, but no attention was paid to these significant differences. Our study, focusing on the overall Gangotri Glacier System as well as all individual glaciers explains these large differences. The higher mass wastage on Gangotri Glacier (-0.88 ± 0.31 m w. e. aˡ) is compensated by the slight mass loss on Meru Glacier (−0.17 ± 0.29 m w.e. a −1 ) and positive mass budgets on fragmented high-altitude tributary Chaturangi and Raktavaran glaciers (0.49 ± 0.17 and 0.62 ± 0.15 m w. e. aˡ, respectively) hence the whole Gangotri Glacier System has been experiencing a moderate mass loss of -0.27 ± 0.25 m w. e. aˡ over 1979-2020. This mass wastage is similar to a nearby debriscovered (13% debris surface) Dokriani Bamak Glacier that witnessed a mean mass wastage of -0.25 ± 0.37 m w. e. aˡ over 1979-2018 (Azam and Srivastava, 2020). Further, the modelled mass wastage on Gangotri Glacier System is also in good agreement with the mean geodetic mass wastage of

FIGURE 10
Panels (A-E): Bubble plots of relative MB sensitivities for all parameters corresponding to T M as 100% for Gangotri, Chaturangi, Raktavaran, Meru and Gangotri Glacier System, respectively. Panels (F-J): Altitudinal MB sensitivities for all the parameters for Gangotri, Chaturangi, Raktavaran, Meru and Gangotri Glacier System, respectively.

Frontiers in Earth Science
frontiersin.org 15 -0.37 m w. e. aˡ in the whole Himalaya over 1975.
The Himalayan glaciers have been in general wastage conditions, while the glaciers in the Karakoram Range (west to the Himalaya) are in a near-balanced state at least since 1970s Bolch et al., 2019). The near-balanced condition of Karakoram glaciers is linked with the increasing snowfalls (hence more accumulation) due to increased local irrigation Farinotti et al., 2020). With increasing global warming, the Karakoram glaciers also started thinning post-2010, and they are expected to lose their mass in the coming decades (Hugonnet et al., 2021). The positive MBs on Chaturangi and Raktavaran glaciers should not be misunderstood as climate change deniers or compared with Karakoram glaciers. The positive MBs on these fragmented tributary glaciers are due to non-climatic topographic reasons: very high area-altitude distribution and the presence of thick debris cover.

Conclusion
The MBs of Gangotri Glacier, its fragmented tributary glaciers (Chaturangi, Raktavaran and Meru glaciers) and the Gangotri Glacier System were reconstructed using T-index model over 1979-2020. Most of the model parameters were taken from the nearby Dokriani Bamak Glacier except T M and P g , which were derived by calibrating the modelled MB of the Gangotri Glacier System with the available geodetic MB from Bhattacharya et al. (2016) over 2006-2014. The model was validated using the modelled SLAs against the manually delineated snow line altitudes (SLAs) derived from the Landsat satellite images.
The fragmented tributary Chaturangi and Raktavaran glaciers showed positive mean summer and annual MBs due to the high area-elevation distribution and debris-covered tongues. Literature suggests retreat and advance periods on Chaturangi Glacier. Further research, using some sophisticated models including the glacier dynamics and the local topography, is needed to better understand the dynamics of Chaturangi Glacier.
Because of the existence of melting hotspots (ice cliffs and supra-glacial lakes), the debris-covered tongue of the Gangotri Glacier was assumed as a clean-ice surface, and DDF I was employed to compute the ablation. Our model computed an additional ablation contribution of 32% due to ice cliffs and supra-glacial lakes. The altitudinal MB gradients of our study region were determined to be slightly less than those of other Himalayan glaciers due to the debris-covered ablation surface and high altitudes. The modelled MBs are most sensitive to T M followed by LR and other parameters (that have quite low sensitivities). The sensitivity analysis showed that the change in annual MBs due to 1°C rise in temperature can be compensated by increasing the precipitation by 83%, 32%, 24%, 66%, and 59% for Gangotri, Chaturangi, Raktavaran, Meru glaciers and Gangotri Glacier System, respectively.
Our study for the first time showed that there are highaltitude, heavily debris-covered fragmented tributary glaciers that might be facing balanced or positive mass budgets in the Himalaya. Further, we also showed that a simple T-index model, with some limitations, can be applied to simulate the MBs of a debris-covered glaciers having ice cliffs and supraglacial lakes. We stress that the positive MBs on fragmented tributary Chaturangi and Raktavaran glaciers should not be considered as climate change deniers as these MBs are due to topographic settings, and with increasing global warming these glaciers are also expected to waste their ice mass.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions
MFA conceived the study. MAH, MFA, and SS developed the model. MAH did the simulations and developed the analysis and figures with the help of MFA and SS. All authors contributed to writing the manuscript.