Challenges in Understanding the Variability of the Cryosphere in the Himalaya and Its Impact on Regional Water Resources

The Himalaya plays a vital role in regulating the freshwater availability for nearly a billion people living in the Indus, Ganga, and Brahmaputra River basins. Due to climate change and constantly evolving human-hydrosphere interactions, including land use/cover changes, groundwater extraction, reservoir or dam construction, water availability has undergone significant change, and is expected to change further in the future. Therefore, understanding the spatiotemporal evolution of the hydrological cycle over the Himalaya and its river basins has been one of the most critical exercises toward ensuring regional water security. However, due to the lack of extensive in-situ measurements, complex hydro-climatic environment, and limited collaborative efforts, large gaps in our understanding exist. Moreover, there are several significant issues with available studies, such as lack of consistent hydro-meteorological datasets, very few attempts at integrating different data types, limited spatiotemporal sampling of hydro-meteorological measurements, lack of open access to in-situ datasets, poorly accounted anthropogenic climate feedbacks, and limited understanding of the hydro-meteorological drivers over the region. These factors result in large uncertainties in our estimates of current and future water availability over the Himalaya, which constraints the development of sustainable water management strategies for its river catchments hampering our preparedness for the current and future changes in hydro-climate. To address these issues, a partnership development workshop entitled “Water sEcurity assessment in rIvers oriGinating from Himalaya (WEIGH),” was conducted between the 07th and 11th September 2020. Based on the intense discussions and deliberations among the participants, the most important and urgent research questions were identified. This white paper synthesizes the current understanding, highlights, and the most significant research gaps and research priorities for studying water availability in the Himalaya.

The Himalaya plays a vital role in regulating the freshwater availability for nearly a billion people living in the Indus, Ganga, and Brahmaputra River basins. Due to climate change and constantly evolving human-hydrosphere interactions, including land use/cover changes, groundwater extraction, reservoir or dam construction, water availability has undergone significant change, and is expected to change further in the future. Therefore, understanding the spatiotemporal evolution of the hydrological cycle over the Himalaya and its river basins has been one of the most critical exercises toward ensuring regional water security. However, due to the lack of extensive in-situ measurements, complex hydro-climatic environment, and limited collaborative efforts, large gaps in our understanding exist. Moreover, there are several significant issues with available studies, such as lack of consistent hydro-meteorological datasets, very few attempts at integrating different data types, limited spatiotemporal sampling of hydro-meteorological measurements, lack of open access to in-situ datasets, poorly accounted anthropogenic climate feedbacks, and limited understanding of the hydrometeorological drivers over the region. These factors result in large uncertainties in our estimates of current and future water availability over the Himalaya, which constraints the development of sustainable water management strategies for its river catchments hampering our preparedness for the current and future changes in hydro-climate.
To address these issues, a partnership development workshop entitled "Water sEcurity assessment in rIvers oriGinating from Himalaya (WEIGH)," was conducted between the 07th and 11th September 2020. Based on the intense discussions and deliberations among the participants, the most important and urgent research questions were identified. This white paper synthesizes the current understanding, highlights, and the most significant research gaps and research priorities for studying water availability in the Himalaya.

WATER IN THE HIMALAYA ALONG DIFFERENT DIMENSIONS
The Himalaya consists of the highest mountains in the world, its orogeny leading to unique geomorphology has given rise to many complex processes such as river meandering, snow-glacierpermafrost interactions, Indian winter and summer monsoons (c.f. Figure 1; Bookhagen and Burbank, 2006;Dimri et al., 2015). The Himalaya is also known as the Third Pole as they contain a large amount of frozen water outside the polar region, which together with the Indian summer and winter monsoons Maharana et al., 2019), ensures that the rivers originating from the Himalaya plausibly never run dry (c.f. Figure 1; Bolch et al., 2019a;Immerzeel et al., 2020). Therefore, these river systems are crucial for supporting the agrarian societies and hydro-power needs of the downstream nations (Dame and Nüsser, 2011;Nüsser et al., 2012;Azam et al., 2021). Precipitation and temperature gradients across and along the Himalaya in conjunction with ridge-valley floor interactions lead to complex manifestation and physical orographic forcings (Bookhagen and Burbank, 2010;Thayyen and Dimri, 2018;Banerjee et al., 2020Banerjee et al., , 2021. The occurrence and frequency of snowfall replenishes the water storage over the Himalaya (Daloz et al., 2020). Downstream into the Indian subcontinent, the Indian summer monsoon precipitation becomes the dominant driver for water mass transport on and below the surface. It stays for a short time in surface water bodies and longer in subsurface water bodies. However, groundwater exploitation over the northwest Indian region in recent decades is changing the residence time of groundwater and affecting water availability (Bookhagen and Burbank, 2010;Singh et al., 2019).
The availability of water also varies both spatially and temporally. Water leaves the Himalayan river basins in the form of water vapor (driven by evapotranspiration) and river discharge. The majority of these river systems, such as the Indus, Ganga, and Brahmaputra, feed multiple nations (Figure 1), and thus the sharing of river discharge data is a politically sensitive issue (Molden et al., 2017). In addition, accurate evapotranspiration modeling requires an extensive network of weather stations, while only a limited number of weather stations are operational, especially near the glacierised regions where topographical challenges are inherent . Hence, remote sensing and gridded/reanalysis datasets becomes one of the preferred tools for estimating hydro-meteorological observables such as, air temperature, precipitation, evapotranspiration, but they come with their own challenges. For example, they may have a relatively coarse spatial resolution: the Tropical Rainfall Measuring Mission (TRMM; refer to Table A1 for abbreviations) precipitation products are at ∼30 km grid (Hofmann-Wellenhof et al., 2007), and the fifth generation European Center for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis ERA5/ERA5-Land data at ∼30/∼10 km grid (Hersbach et al., 2020). Also, several satellite observations are an indirect measure of a physical process and they often require calibration and validation against in-situ observations that are at times not available (Karimi and Bastiaanssen, 2015;Zhang et al., 2016). Further, the length of satellite era is relatively short, and time-series from many of the satellite missions cannot be considered as climate data (Vishwakarma et al., 2021). Satellite products require substantial post-processing to remove atmospheric effects and improve the signal to noise ratio. Due to these challenges, model simulated outputs at high spatiotemporal resolution and longer temporal coverage are also being employed to understand the past, the present, and the future of water availability in the Himalaya (Khanal et al., 2021). For example, the drivers of variability in precipitation have been investigated with the help of global and regional climate model outputs, but no conclusive pattern or trend can be identified (Immerzeel et al., 2010(Immerzeel et al., , 2020Kaser et al., 2010;Lutz et al., 2016;Maharana and Dimri, 2016;Dimri et al., 2018;Huss and Hock, 2018;Maharana et al., 2019). Over the Hindu-Kush Karakoram Himalayan range, Palazzi et al. (2013) have shown decreased summer precipitation, whereas Kulkarni et al. (2013) and Javed et al. (2022) have shown increased winter precipitation. Dimri et al. (2022) explained the elevation dependent warming over the Karakoram region.
Modeling framework-based studies over the Himalaya show different impacts on different river basin systems (Ghimire et al., 2018;Nengker et al., 2018). Atmospheric changes associated with global warming can lead to alteration of the mountain hydrological cycle, in turn the streamflow patterns of the Himalayan rivers (Archer et al., 2010;Lutz et al., 2016;Armstrong et al., 2018). Available basin-wide studies show that glacier melt contributes the largest to total water discharge in the upper Indus basin (41%), whereas rainfall-runoff in the upper Ganga (66%) and Brahmaputra (59%) basins (Lutz et al., 2014). Glacier melt contribution is relatively lower in the upper Ganga (12%) and Bramhaputra (16%) basins. Snowmelt contribution is also lower FIGURE 1 | Map of the Himalayan River Basins-Indus, Ganga and Brahmaputra-feeding from the glaciers. The glacier cover (blue patches) is based on Randolph Glacier Inventory 6.0 (RGI Consortium, 2017). The background is the Stamen terrain map, the Himalayan region subdivision from Bolch et al. (2019a), river basin boundaries are based on the International Center for Integrated Mountain Development (ICIMOD, 2021), and river central lines are from the Global High-Resolution River Centerlines (USGS). RGI's first order regions 14 (South Asia West) and 15 (South Asia East) are also shown (red outlines).
in the upper Ganga (9%) and Bramhaputra (9%) basins, whereas 22% in the upper Indus basin (Lutz et al., 2014). Snowmelt contribution to the total discharge is significantly higher (or highest) in the headwaters of small/medium sub-basins in the Himalayan region, for example, in the Liddar Basin, western Himalaya (60%; Jeelani et al., 2017), Bhagirathi, upper Ganga Basin (60%; Rai et al., 2019), Chandra, a sub-basin of Chenab River (60%; Singh A. T. et al., 2021) etc. Groundwater also contributes significantly to the river discharge, but only noted in the Chandra basin . Nevertheless, the contribution from snow and glacier melt are expected to change in the future due to ongoing climate change (Prasch et al., 2013;Ragettli et al., 2016;Khadka et al., 2020). This is evident in the Baspa basin, a sub-basin of Satluj river, where  applied a physical glacio-hydrological model model and noted an increase in snow melt contribution by 6% between 2000 and 2018. Lutz et al. (2016) concluded that hydrological extremes would be frequent and stronger in a warmer climate, but the uncertainties are large. Based on stable isotope measurements of hydrogen and oxygen together with historical runoff records in the Indus River Basin, Karim and Veizer (2002) have shown a major loss of water to evapotranspiration, which is in line with the predictions in a warming climate (Banerjee and Azam, 2016;Kraaijenbrink et al., 2017). There are observational and modeling challenges over the Himalaya due to the paucity of observations (You et al., 2017).
In general, Himalayan glaciers have been retreating and losing mass during the last decades at a significant rate Azam et al., 2018;Bolch et al., 2019a;King et al., 2019;Maurer et al., 2019;Bhattacharya et al., 2021). It is projected that river runoff in the larger catchments will increase up to 2050s but then decrease sharply outside the monsoon season owing to the decimated Himalayan glacier area and volume (Lutz et al., 2014;Rounce et al., 2020;Azam et al., 2021;Chandel and Ghosh, 2021). Therefore, we can expect to see an increase in extremes of water availability resulting in disasters such as glacial lake outburst flood (GLOF), avalanches, and droughts (Eckstein et al., 2017;Furian et al., 2021;Majeed et al., 2021;Shugar et al., 2021). GLOF is a sudden release of water or flash flood from a moraine-or ice-dammed glacial lake due to dam failure (Shugar et al., 2020;Veh et al., 2020;Sattar et al., 2021). GLOF could be triggered by a number of processes, including ice/debris/rock fall, avalanches, earthquakes, internal piping, etc. (Veh et al., 2020). GLOF mechanisms can be understood by simulating the GLOF chain sequences and flood routing modeling. Understanding these processes is critical for risk management and developing mitigation strategies for the downstream communities .
A majority of the hydrometeorological and glacier-related processes (c.f. Figure 2), e.g., glacier mass balance, runoff, GLOF modeling studies rely on outputs from regional climate models, land surface models, and remote sensing datasets which suffer from high uncertainties over the Himalaya (Lutz et al., 2016; FIGURE 2 | An illustration of simplified hydrological cycle showing the important cryospheric and atmospheric processes, variables and features. The figure is inspired from Azam et al. (2021). Background illustration credit: Nouri Atchabao/Vecteezy. Hock et al., 2019;Gusain et al., 2020;Azam et al., 2021;Nepal et al., 2021). Furthermore, recent studies indicate that human activities, both agricultural and land use change, that are missing from models have significantly impacted the regional water budget (Budakoti et al., 2021;. Therefore, studies based on models can only provide a limited insight. Many attempts have been made to quantify the amount of water entering and leaving these river catchments and factors influencing their variability, but still our knowledge is limited (Lutz et al., 2016;Pritchard, 2019). Therefore, we need to obtain more data-driven inferences by including as many types of observations as possible.
In the following sections, various techniques, observations, and methods that are being or can be employed to understand the water cycle in the Himalaya are discussed.

ESTIMATING GLACIER MASS BALANCE
An accurate assessment of glacier mass balance is vital for understanding the vulnerability of glaciers due to climate change and the rates of glacier denudation (Rabatel et al., 2011). The mass balance of Himalayan glaciers has been measured with different methods, from conventional (Østrem and Stanley, 1969), to glaciological, to remote sensing and geodetic , to various modeling techniques (Bolch et al., 2012;Shea et al., 2015;Azam et al., 2018).

In-situ Measurements
Glaciers respond to even small changes in temperature and precipitation (Oerlemans, 2001). Therefore, field observations of glacier mass balance at annual and seasonal scales are excellent climate records (Zemp et al., 2009;Trewin et al., 2021). Unfortunately, the high altitudes and harsh weather conditions hinder glaciological measurements in the Himalaya; therefore, existing in-situ studies mostly consider the easily accessible, small sized, and less debris-covered glaciers Bolch et al., 2019a).
Glacier mass balance is the sum of all ablation (melt, sublimation, snowdrift, calving, etc.) and accumulation (snowfall, snowdrift, avalanches, etc.) processes (illustrated in Figure 2). These processes are glacier-wide and continuous in time, however, for practical reasons they are observed at point scale. Ablation is measured by installing an ablation stake network over the ablation area of the glaciers and assuming an ice density of 900 kg/m 3 (Kaser et al., 2003;Huss, 2013). In contrast, snow accumulation is measured with the help of snow pits/cores and probing in the accumulation zone and measuring the density of snow along with the depth to the previous year's surface Azam et al., 2016;Mandal et al., 2020;Soheb et al., 2020). These point-scale observations are integrated to estimate the glacier-wide mass balance using glacier hypsometry data obtained from satellite imagery or terrestrial sources (Cuffey and Paterson, 2010). Glacier hypsometry is generally obtained from a Digital Elevation Model (DEM) and the glacier outline (preferably from the recent ablation-season optical satellite images). Since the in-situ observations are important but limited in number to provide a comprehensive overview over the Himalaya, researchers rely on other proxies to estimate snowmelt and glacier runoff.
Due to the difficulties in field measurements in the Himalaya, only about 36 glaciers have glaciological mass balance measurements conducted for 1 year or more ( Table 1). The majority of the measured glaciers are of small to medium size. Of all the 36 measured glaciers, most are located in the western Himalaya (n = 18), followed by the central Himalaya (n = 14), eastern Himalaya (n = 3), and Karakoram (n = 1). About ten glaciers have a series of measurements equal to or longer than 10 years ( Table 1). The longest continuous series spans 17 years for the Chhota Shigri Glacier in the western Himalaya (−0.46 ± 0.40 m w.e. a −1 over 2002-2019; Mandal et al., 2020). The Naimona Nyi Glacier in the northern slope of central Himalaya has the longest record of annual mass balance but discontinuous .
By averaging all available in-situ mass balance records for 36 glaciers (240 annual mass balance data points), the mean mass balance for the period 1975-2020 becomes −0.60 m w.e. a −1 (standard deviation; SD: 0.56 m w.e.). This region-wide average is not much different from the previous estimate of −0.49 m w.e. a −1 (based on mass balance records from 24 glaciers) for the period 1975-2015 . This increased number of glacier measurements certainly improves the representativeness of the region-wide mass balance for the Himalaya. However, considering the topo-climatic heterogeneity across the Himalaya, 36 glaciers is a small sample size of roughly 95,000 total number of glaciers in the region.

Remote Sensing Based Estimates
Various remote sensing-based methods are available to estimate glacier mass balance at different spatiotemporal scale ( Table 2). Among the available methods, the geodetic method is the most widely used, which uses multiple DEMs, from at least two different epochs, to obtain a volume change estimate (Bamber and Rivera, 2007;Gardelle et al., 2012Gardelle et al., , 2013Brun et al., 2017;Hugonnet et al., 2021). These DEMs can be obtained from various remote sensing products like (a) satellite-based optical stereo images, (b) satellite-based radar interferometry, (c) aerial photogrammetry, images from drones, or laser scanning, and (d) Uncrewed Aerial Vehicles (UAV) based data through Structure from Motion techniques (SfM).
The DEMs obtained from different techniques have different accuracy, spatial coverage, spatial resolution and procurement cost. Considering the higher vertical errors in publicly available satellite-based DEMs, they are, in general, most useful for decadal-scale glacier mass balance estimations and less for the annual mass balance Braun, 2016, 2018). With the availability of new generation satellites that have better instruments, very high-resolution DEM data can be generated with reasonable accuracy but only for the past few years. It is to be noted that these products are only available on demand at a price, that too at irregular intervals. A few prominent remote sensing missions that have been used widely for geodetic mass balance studies in the Himalaya are: Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) mission Bisset et al., 2020;Hugonnet et al., 2021), Satellite Pour l'Observation de la Terre (SPOT) , Pléiades (Berthier et al., 2014;Wagnon et al., 2021), WorldView/GeoEye , etc.
A higher temporal resolution of mass balance estimates has been generated by using the laser altimetry-based Ice, Cloud, and Land Elevation Satellite (ICESat), ICESat-2 mission data (Kääb et al., 2012(Kääb et al., , 2015Wang et al., 2021) and CryoSat-2 swath altimetry data (Jakob et al., 2021). However, ICESat based mass balance estimates have been shown to suffer from a bias due to the short satellite time span (Neckel et al., 2014;Brun et al., 2017;Wang et al., 2021).
Another widely used method is a simple empirical approach called Accumulation Area Ratio (AAR) and Equilibrium Line Altitude (ELA) method proposed by Kulkarni (1992). AAR is the ratio of the area of the accumulation zone to the area of the glacier (Cogley et al., 2011), whereas ELA is the line that divides the glacier into ablation and accumulation areas where accumulation equals ablation over a balance year (Braithwaite and Raper, 2009). This method has considerable uncertainty (Pratap et al., 2016;Tawde et al., 2016). This was first applied on the Gara Glacier (from 1977-1978 to 1982-1983) and the Gor-Garang Glacier (from 1976-77 to 1983, and then has been used to estimate the surface mass balance of Himalayan glaciers (e.g., Kulkarni et al., 2004;Mir et al., 2014 andTak andKeshari, 2020). Tawde et al. (2016) developed an improved AAR method in which they used satellite data and a temperature-index model. They implemented a regression between the modeled AAR and fieldbased mass balance to estimate glacier-wide mass balance. This improved AAR method was first demonstrated over the Chhota Shigri Glacier and 12 other glaciers in the Chandra basin (from 2000 to 2009). Following the same methodology, Tawde et al. (2017) estimated the mass balance of 146 glaciers in the Chandra basin (from 1984 to 2012), however these had higher uncertainties.
In another development, Rabatel et al. (2005) estimated annual mass balance using variation between ELA at steady state of a glacier and satellite-derived ELA along with mass balance gradient across the ELA. They demonstrated this method over French alpine glaciers and showed that the mass balance of individual glaciers in a basin can be estimated reliably by using field-based observation of at least one representative glacier in the same basin. In the Himalaya, Chandrasekharan et al. (2018) and Garg et al. (2021) have used this method for estimating mass balance over some glaciers. This approach seems to have good potential for annual mass balance estimation. However, it needs to be tested rigorously for wider application. Dumont et al. (2012) reported that the annual minimum albedo average on the whole glacier at the end of the ablation season is strongly correlated (r 2 > 0.95) with the field-based annual mass balance of the Saint Sorlin Glacier in Saint Sorlin Glacier in France. Based on this correlation, the annual mass balance series of the glacier was reconstructed from  Bhutiyani (1999) + corrected by Zaman and Liu (2015); 2, Kaul (1986); 3, Geological Survey of India (2001); 4, WIHG (submitted to WGMS); 5, Soheb et al. (2020); 6, Sangewar and Sangewar (2007); 7, Shrivastava et al. (2001); 8, Angchuk et al. (2021); 9, NCPOR (submitted to WGMS); 10, Geological Survey of India Report (2017); 11, This study; 12, Raina et al. (1977); 13, Geological Survey of India (1992); 14, Koul and Ganjoo (2010) Laha et al. (2017); 18, Gautam and Mukherjee (1992); 19, Geological Survey of India (1991); 20, Yao et al. (2012); 21, Gurung et al. (2016); 22, Baral et al. (2014); 23, Liu et al. (1996) + Yao et al. (2012; 24, Sunako et al. (2019); 25, Fujita and Ageta (2001); 26, Wagnon et al. (2013) + Sherpa et al. (2017; 27, Raina (2009); 28, Tshering and Fujita (2016); 29, Royal Govt. of Bhutan (submitted to WGMS).  Brun et al. (2015) have successfully tested this method for the Chhota Shigri Glacier in the western Himalaya. However, its efficacy was low when applied over the Mera Glacier in the central Himalaya due to the presence of cloud cover, short fieldbased mass balance data, and availability of fewer pixels for average albedo estimation. A similar approach was attempted by Sirguey et al. (2016) to estimate the annual as well as the seasonal mass balance of the Brewster Glacier in New Zealand, where they observed a strong correlation (r 2 = 0.93) between the average of annual minimum albedo over the whole glacier and annual mass balance, as well as a good correlation (r 2 = 0.86) between cumulative winter albedo and winter mass balance.
Of the techniques discussed above, the geodetic mass balance method based on DEM differencing has been the most widely used over Himalayan glaciers. One initial study for deriving geodetic mass balance in the Himalaya is by Berthier et al. (2007). The authors used SRTM v2 and SPOT5 DEMs for estimating mass balance of glaciers around the Chhota Shigri Glacier located in Lahaul-Spiti region from 1999 to 2004 and obtained mass loss rates of −0.7 to −0.85 m w.e. a −1 . The authors present a method to co-register the two DEMs but did not consider properly the SRTM radar penetration into ice. The authors simply assumed the SRTM DEM which was acquired in February 2000 to be the representative for end of the ablation season in 1999 due to some penetration into snow. However, subsequent studies found that the penetration of the SRTM-C band radar beam can be up to few meters (Gardelle et al., , 2013Kääb et al., 2012). Vincent et al. (2013), therefore, revisited the results but found only a difference of less than 0.1 m w.e. a −1 with respect to Berthier et al. (2007).
In the last two decades several studies have been published investigating the glacier mass balance starting from 2000 to recent years based on DEM differencing in various regions of the Himalaya. One of the earliest studies is Gardelle et al. (2013) who investigated eight different regions from the Pamir, to the Himalaya, to Hengduan Shan for the period 1999 to 2011. The authors considered the penetration based on the  Table 3). The mass balance estimates for HMA were further refined by Shean et al. (2020) who used DEMs derived from ASTER and very high-resolution stereo images such as Worldview-2/3. They report a higher mass loss in the Himalayan region during 2000-2018 as compared to the previous studies, with the highest rates over the eastern Himalayan region (Table 3). Recently Hugonnet et al. (2021), further improved the processing of ASTER DEM and related mass balance estimates and provided not only an estimate for the period 2000 to 2019 for all glaciers on Earth but also for four subperiods therein. They found an increasing mass loss since 2000 with higher overall mass loss as compared to the previous studies (Table 3).
In order to better understand the long-term response of glaciers to climate change, longer and better data coverage is required. Declassified US spy stereo satellite images from the 1960s and 1970s provide an excellent opportunity in this regard. 1960s and 1970s high resolution (7-2 m) Corona KH-4 and KH-4B data were first used by Bolch et al. (2008Bolch et al. ( , 2011 along with ASTER and Cartosat-1 data to investigate the mass balance evolution of the Mt Everest region in the Himalaya. 1970s Hexagon KH-9 metric camera data (8 m), first employed in HMA in the Tien Shan by Surazakov and Aizen (2010) and Pieczonka et al. (2013), was later used to improve the knowledge about glacier mass loss in the Himalaya and other regions of HMA as well. Zhou et al. (2018) applied these images along with the SRTM DEM for estimating the mass balance in some regions of the Tibetan Plateau and the Himalaya from 1975 to 2000. They reported a strong negative mass balance in most of the regions but only small negative values for Lahaul-Spiti (Table 3). This is in line with Mukherjee et al. (2018) who found similar insignificantly negative values for the region using 1971 Corona KH-4B and the SRTM DEM. Maurer et al. (2019) applied ∼1975 Hexagon KH-9 and ASTER DEMs, while King et al. (2019) generated Hexagon KH-9 DEMs along with the 2000 SRTM and the ∼2015 HMA DEM (Shean, 2017) to estimate the mass balance before and after 2000. Both studies found significant increase in mass loss after 2000 ( Table 3). King et al. (2020) used declassified imagery along with aerial images and Cartosat-1 data for up to 7 subperiods during 1962-2018 and reported increasing mass loss from −0.23 ± 0.13 m w.e. a −1 for 1962-1969 to −0.38 ± 0.13 m w.e. a −1 for 2009-2018. Bhattacharya et al. (2021) used various available stereo satellite data (in particular Corona, Hexagon and Pleiades data) to report mass changes for up to six subperiods within 1964 to 2019. The highest increase in mass loss was found for the Langtang region with −0.20 ± 0.09 m w.e. a −1 for 1964-1977 to −0.59 ± 0.14 m w.e. a −1 for 2017-2018.

Glacier Melt and Mass Balance Modeling
Glacier surface mass balance and melt modeling has gained considerable attention in recent years, in part because of the availability of high-resolution gridded climate datasets (e.g., ERA5, ERA5-Land, HAR v2, CORDEX, etc.), and highresolution satellite and DEM datasets . Also, the theoretical and physical understanding of various complex snow and glacier processes, such as influence of debris cover, long-term dynamic change, has improved significantly (Reid and Brock, 2010;Shea et al., 2015;Carenzo et al., 2016). Various surface melt and mass balance models have been implemented in the Himalayan region, such as hydrological models (Bhutiyani, 1999;Immerzeel et al., 2012), temperature-index model (Azam et al., 2014a), enhanced temperature-index model (Litt et al., 2019), albedo model (Brun et al., 2015), surface energy balance (SEB) model (Azam et al., 2014b), distributed SEB model (Arndt et al., 2021;Steiner et al., 2021;Srivastava and Azam, 2022) Rounce et al., 2020]. However, the majority of the modeling studies use temperature-index models , with some modification for better representation (e.g., Pellicciotti et al., 2005). The temperatureindex model is also known as the degree day model (DDM), because in them daily melt depth is computed by multiplying the number of positive degree days (often called cumulative positive degree days; CPDD) by the melt factor (Hock, 2003)This model is simple and can be applied over large areas/basins with limited data input requirements (Lutz et al., 2014;Azam et al., 2021). SEB models account directly for many of the physical processes that affect melt (e.g., surface melt, sublimation), but they require observations of several variables such as air/surface temperature, air pressure, short-wave and long-wave radiations, wind speed and humidity, etc. (Wagnon et al., 1999;Favier et al., 2004;Litt et al., 2019). Due to the scarcity of input data needed for SEB model, they are generally not suitable for modeling glacier mass balance over the Himalaya (Azam et al., 2014b).
Glacier mass balance models have been developed for a single glacier, for the entire basin as well as for the entire Himalayan region depending on the forcing and calibration data availability. For example, recently, Srivastava and Azam (2022) modeled the mass balance of Chhota Shigri (western Himalaya) and Dokriani Bamak glaciers (central Himalaya) for the past seven decades (1950-2020) using a temperature-index model forced with bias-corrected ERA5 data. Khadka et al. (2020) applied OGGM by integrating Glacio-hydrological Degree-day Model (GDM; Gupta et al., 2019) in the Koshi River basin (∼3,000 km 2 glacierised area) in the central Himalaya to estimate the glacier mass balance and the runoff. The model was forced with insitu hydro-meteorological data and historical gridded data from the Climatic Research Unit (CRU). The model performance was good in simulating runoff (r > 0.8) and partitioning various runoff components. Rounce et al. (2020) applied the PyGEM model, forced with geodetic mass balance observations over the entire HMA region. Based on the model projection, they showed that by the end of the century, glaciers in HMA might lose between 29 ± 12% (RCP 2.6) and 67 ± 10% (RCP 8.5) of their total mass relative to 2015.
Glacier mass balance and melt models are extremely helpful to understand the glacier-scale to region-scale glacier mass balance behavior. However, in general, model outputs have large uncertainties because they are sensitive to the meteorological input variables and model parameters (Lutz et al., 2016;Rounce et al., 2020). Despite large uncertainty they have considerably improved our overall understanding of spatiotemporal variability in the Himalayan water resources.

Satellite Gravimetry
The Gravity Recovery and Climate Experiment (GRACE) mission and its successor GRACE Follow On (GRACE-FO, launched in 2018 after the end of GRACE in 2017) are a pair of identical satellites co-orbiting in the same plane (Tapley et al., 2019). The two satellites measure the change in the distance between them via a microwave link, which is a function of the gravity field of the Earth. The GRACE mission provided changes in the gravity field of the Earth at monthly and sub-monthly scale, which can be related to the surface mass changes (Tapley et al., 2019).
The gravity field can be recovered from the satellite mission either as blocks of masses (also called mascons) or as spherical harmonic coefficients (Tapley et al., 2019). The unconstrained spherical harmonic coefficients are noisy and require filtering prior to using GRACE data (Swenson and Wahr, 2006), but filtering introduces leakage and diminishes the signal amplitude (Vishwakarma et al., 2016(Vishwakarma et al., , 2017. Metrics have been designed to understand the efficacy of filtering , and methods have also been developed to account for leakage and the amplitude loss (Vishwakarma et al., 2018). It must be noted here that the mascon solutions are also filtered versions of the GRACE data, where filtering is performed with a regularization scheme that aims to minimize signal leakage with the help of prior information.
One of the most contentious issues with the use of GRACE-FO is its spatial resolution. A wide range of numbers have been advocated for indicating the spatial resolution, for example, Longuevergne et al. (2010) proposed a value of 200,000 km 2 and Vishwakarma et al. (2018) indicated 63,000 km 2 if the application can tolerate an error of 2 cm. It was also shown that the spatial resolution depends on the method used for filtering and leakage correction (Vishwakarma et al., 2018). For the spherical harmonic solutions Devaraju and Sneeuw (2016) devised the method of modulation transfer functions to identify the spatial resolution. However, for the case of mascons, the spatial resolution is believed to be at a spatial scale of 300 km. Despite these issues, it is known that GRACE-FO can observe variations in the gravity even at small spatial scales, if the gravity anomaly is large enough (Sneeuw and Sharifi, 2015). In this context, it has been shown that small spatial scale variations require dedicated modeling of the range acceleration signal obtained from the level 1-b data of the satellite mission (Weigelt, 2017;Ghobadi-Far et al., 2020).
The GRACE derived total water storage (TWS) anomaly includes changes in the soil moisture, canopy water, surface water, snow water, and groundwater (Tapley et al., 2019). Therefore, to arrive at one component, such as either ice mass balance or groundwater change, it is necessary to remove other storage components usually obtained from a model simulation that can have high uncertainty (Tiwari et al., 2009;Chen et al., 2015;Long et al., 2016). Therefore, using GRACE data for estimating mass changes over a glacier, which is much smaller than the spatial resolution of the observation, or in a groundwater aquifer is challenging (Wang Q. et al., 2018;Vishwakarma, 2020). Nevertheless, there have been some efforts in quantifying glacier mass change in the HMA using the point mass modeling approach, which gives an estimate for the whole region but with very high uncertainty. For example, Jacob et al. (2012) estimated the mass loss at a rate of −4 ± 20 Gt/yr from 2003 to 2010, while a recent study by Wouters et al. (2019) reported a trend of −13.5 ± 6 Gt/yr for the complete GRACE time series.
GRACE has been found more useful at studying basin-scale hydrology. While there are spatial and temporal resolution issues with GRACE data, it has proved useful in downstream regions where groundwater is being exploited rapidly . In addition, like any other remote sensing product, there is a need to validate the GRACE data with groundwater monitoring data from government agencies, farmers and stakeholders. In the past, the use of crowd sourced farmer data for groundwater management using the Mywell App has been successful in creating a high-resolution database for groundwater level data at the village scale (Maheshwari et al., 2014;Maheshwari, 2020). Preliminary results indicate that the locals can aid in better management of groundwater resources and can also augment rural database by providing crowd sourcing of data along with GRACE/GRACE-FO data.
Another geodetic technique that has recently emerged as a tool for monitoring water resources near glaciated regions is Global Navigation Satellite System (GNSS). It has been used for monitoring glacier mass change, lake water level, and surface deformation due to hydrological mass change (Dunse et al., 2012;Durand et al., 2019;Elliott and Freymueller, 2020;Koulali et al., 2022).

Global Navigation Satellite Systems
where ϕ is the carrier phase, ρ is the true range between the satellite and the receiver, λ is the wavelength of the carrier, N is the integer number of cycles, c is the velocity of light, δt is the clock error, T is the delay due to the troposphere, and I is that due to the ionosphere, M is the multipath error, S is the combined term of all other systematic errors like antenna phase center error and hardware delays in the satellite and the receiver, and ε is the random error associated with the phase measurement. The term ρ contains the position of the receiver and the satellite, and the terms T and M provide important environmental information.
In GNSS monitoring of glaciers, the position information from GNSS is sensitive to the movement of the glacier and surface elevation changes due to the load changes associated with glacier mass change. The tropospheric delay term T provides information concerning the precipitable water vapor (PWV), which is typically absorbed into weather models for forecasting. The multipath term M is the reflected GNSS signal. Since the GNSS uses a microwave L-band for sending signals, the reflected signal contains information concerning the backscattering properties of the surface. This has been successfully used for retrieving snow height information in glaciers and mountainous areas alike (Larson, 2016). Thus, GNSS can provide a variety of vital information for monitoring snow and glaciers.
The observation strategies vary for the different types of variables that can be retrieved from GNSS positioning. For understanding loading due to glacier mass change, it is typical to have a continuously operating network of GNSS stations in and around the glacier (Jacobsen and Theakstone, 1997).
Such networks can also provide PWV estimates, but to perform GNSS reflectometry with them, the sites with reasonably flat terrain have to be chosen (Larson, 2016). For monitoring glacier movement and shape, real-time kinematic (RTK) GNSS measurements can be conducted in periodic campaigns. DEMs can also be constructed to identify surface elevation changes of the glacier (Jacobsen and Theakstone, 1997). With the advent of low-cost GNSS receivers and antennas, more observations are being retrieved that are used for deriving novel inferences. For example, Koch et al. (2014) and Henkel et al. (2018) planted a low-cost GNSS antenna on a snow-free area in the summer months and estimated the winter snow depth from distorted GNSS signals. The low-cost GNSS receivers perform as good as the geodetic receivers (Odolinski and Teunissen, 2016), and therefore, they can be used for the densification of GNSS networks in the snow and glaciated areas.

GLACIER ICE THICKNESS MODELING
Glacier ice thickness and ice reserves are estimated by insitu measurements, power-law relations, physical models, and artificial neural networks (ANN) (Zemp et al., 2019;Haq et al., 2021;Millan et al., 2022). In-situ measurement of ice thickness such as borehole drilling, radio-echo sounding and ground-penetrating radar (GPR) is often costly, logistically demanding, and unsafe due to the rugged nature and harsh climatic conditions of the glaciated terrain (Saintenoy et al., 2013;Bohleber et al., 2017). The next approach is based on Volume Area Scaling (VAS) (Bahr, 1997). VAS is suitable for obtaining a general ice thickness pattern for most of the glaciers (Bahr et al., 2015), but requires calibration from a DEM and estimates of glacier extents (Möller and Schneider, 2010;Laumann and Nesje, 2017;Banerjee, 2020). It is to be noted that VAS cannot give information regarding the distribution of ice thickness. As a result, glaciologists have begun to rely on ice thickness models that use multi-temporal, multi-spectral and high spatial resolution satellite data to represent varying complexities of the glacial environment (Farinotti et al., 2021).

Types of Glacier Ice-Thickness Models
Recent work by the Ice Thickness Models Intercomparison eXperiment Phase 1 and Phase 2 (ITMIX1 and 2) evaluated all the models developed so far for the mountain glaciers. Several models exist that rely on different principles, viz., (1) using the shallow ice approximation (SIA) and an empirical relation between glacier elevation range and basal shear stress (Haeberli and Hoelzle, 1995) at the local scale; (2) flowlinebased approach considering mass conservation and Glen's ice flow law (Glen, 1955); (3) two-dimensional consideration of the continuity equation (Morlighem et al., 2011). Apart from these, recent developments in data sciences have led to models that that do not rely on the physics of the problem, such as the Artificial Neural Network (ANN) approach (e.g., Clarke et al., 2009;Mey et al., 2015;Haq et al., 2021). Table 4 provides an up-to-date list of ice thickness models that have been used in the Himalaya. Other approaches Artificial neural networks (Clarke et al., 2009;Haq et al., 2021) (No physical laws followed)

Glacier outline and Digital elevation model of the surface (OL + DEM); Surface mass balance (SMB); Surface ice flow velocity (Vel.); Rate of elevation change ( ∂h ∂ t ). #$ Represent type of model approaches.
∧Perfect plasticity assumptions + local scale empirical relation between basal shear stress and elevation range. ∧@

Recent Advancements in Ice Thickness Modeling
The applicability of ice thickness models has evolved significantly due the abundance of satellite remote sensing datasets (e.g., LiDAR, Optical and SAR images) to the development of automated models like the Volume and Topography Automation (VOLTA; James and Carrivick, 2016), the Open Global Glacier Model (OGGM; Maussion et al., 2019), and ice velocity and slope-based methods (Millan et al., 2022). On the other hand, like the Glacier Thickness Database (GlaThiDa; Welty et al., 2020)a global archive for in-situ ice thickness datasets and Global Glacier Thickness Initiative (G2TI; Farinotti et al., 2019)-have opened gateways for exploiting the maximum potential of the ice thickness models. Farinotti et al. (2021) have shown that calibrating the ice thickness models with field observations (even limited) and ensemble-approach is beneficial in terms of improving accuracy and robustness compared to individual model estimates. Millan et al. (2022) have recently modeled the ice thickness for all the glaciers on Earth using high-resolution ice velocity and topographic data based on the SIA approximation. The re-evaluated ice volume estimates by Millan et al. (2022) is the most comprehensive estimate because they rely on recent DEMs and satellite remote sensing image. The next best attempt was by Farinotti et al. (2019) that were based on at least a decade old dataset.

Recent Estimates of Ice Thickness Over the Himalaya
Studies related to ice thickness date to the 1960s (Raina, 2009), and field investigations are still ongoing in sporadic nature at accessible glaciers in the Himalaya. However, only a few glaciers are being surveyed using GPR (e.g., Azam et al., 2012;Vincent et al., 2016;Mishra et al., 2018;Pritchard et al., 2020). GPR survey-based ice thickness measurements are available for only around 21 glaciers and that too at different time intervals . Therefore, VAS and ice thickness models with inputs from glacier inventories (e.g., GLIMS, RGI, etc.) and remote sensing paved the path to study glaciers at a relatively larger scale than with the GPR (e.g., Sattar et al., 2019;Pandit and Ramsankaran, 2020).
The VAS approaches were used to estimate glacier ice volume for the entire Hindu-Kush Himalaya range with the help of the RGI glacier outlines. However, the estimation of glacier ice volume varied from one study to the other due to a different set of parameters being used in different studies (Radić and Hock, 2010;Marzeion et al., 2012;Grinsted, 2013). Bolch et al. (2012), Frey et al. (2014, Linsbauer et al. (2016), andFarinotti et al. (2019) were able to estimate the total glacial volume of the HK region using ice thickness models for the year 2000. Bolch et al. (2012) reviewed VAS models and provided an updated estimate based on a shear-stress model. Similarly, Frey et al. (2014) and Linsbauer   Table 5 provides a summary of the estimated ice volume for the Himalayan region from various methods and sources. Estimates of ice thickness is a crucial input for glacier models that are driven by climate models to project future glacier mass balance. Climate model projections of future glacier mass balance over the Himalaya produce consistent results and suggest that glaciers in the region will lose up to 90% of their mass by 2100 (Barsch and Jakob, 1998;Kraaijenbrink et al., 2017;Bolch et al., 2019a;Shannon et al., 2019). Based on these projections, there is considerable concern among policymakers and infrastructure planners about future water supplies to the wider region and the potential for increased glacier-related hazards impacting populations and infrastructures in the high mountain valleys (Huggel et al., 2020;Immerzeel et al., 2020).

CLIMATE MODELS BASED MASS BALANCE OF KARAKORAM-HIMALAYAN GLACIERS
Various parts of the Earth system, such as oceans, atmosphere, cryosphere, land surface dynamics, and human systems, interact and influence each other. Therefore, Global Climate Models (GCMs) were developed to include all possible sub-systems, their interactions, and response to any change in the Earth system.
GCMs are computationally expensive and run at a coarse spatiotemporal resolutions, which reduces their efficacy at a regional scale. Testing and calibrating GCMs require highresolution present-day climate data along with well-constrained information on past climates of critical regions worldwide, which includes changes in the glacier extent and mass balance. So far, the terrestrial cryosphere is represented in a highly simplified way in majority of the state-of-the-art global and regional climate models (RCMs) as they use static glacier masks, i.e., (i) no changes in ice extent, no feedback to the atmosphere; (ii) no consideration of water volume stored; (iii) either no or simplified runoff generation (see Figure 3; Kotlarski et al., 2010;Kumar et al., 2015Kumar et al., , 2019. Therefore, a more sophisticated approach is necessary to overcome the poor representation of the glacier processes in climate models. RCMs, owing to their more nuanced spatiotemporal resolutions, act as better alternatives for regional studies. Within the framework of the Coordinated Regional Climate Downscaling Experiment (CORDEX), a regional earth system model comprising of REgional atmospheric MOdel (REMO), Max Planck Institute's Ocean Model (MPIOM), and a hydrological discharge (HD) model, coupled using OASIS coupler, has been developed for the CORDEX-SA (South-Asia) domain for climate change studies. This setup is called the Regional Earth System Model (RESM; Kumar et al., 2022). The proper description of cryospheric processes is essential for simulating the complete terrestrial water cycle in climate models, especially for the high mountain regions.
Owing to various socio-economic, topographical, and political constraints, attempts for a holistic study of Hindu-Kush Karakoram Himalayan glaciers have been very sporadic, undermining the fact that the region suffers from an acute FIGURE 3 | A schematic representation of the complex glacial processes involved in mass-balance estimations (top left). The difference in simplified (static glacier masks) and sophisticated (dynamical glacier masks) approaches (top right). Description of the dynamical glacier scheme DYNAMICE and the regional model dynamics of REMO glacier (bottom right). Different model dynamics and parameterisation schemes operate separately for glacierised and non-glacierised regions (bottom left).
shortage of quality observational datasets, both spatial and temporal. In RCM REMO, a dynamical glacier scheme (DGS) was implemented that has the unique ability to simulate the glacier mass balance for even dynamic fraction of each grid box depending on the accumulation and ablation conditions, while accounting for direct physical feedback mechanisms. The scheme represents surface glacier cover on a subgrid-scale and calculates the energy and mass balance of the glacierised part of a grid box (Kotlarski, 2007;Kotlarski et al., 2010;Kumar et al., 2015Kumar et al., , 2019). This coupled model system is referred to as REMO glacier . It has been applied for the Himalaya to show computational effectiveness. Recently, several studies (e.g., Kumar et al., 2015Kumar et al., , 2019Javed et al., 2022) have started to assess the impact of largescale and micro-climate on Himalayan glaciers, focusing on area and mass balance changes at a sub-grid scale.

SNOW COVER, SNOW DEPTH, AND SNOW WATER EQUIVALENT VARIABILITY OVER THE HIMALAYA
In addition to ice-related processes, snow cover has significant spatial and temporal variations over the Himalaya. Snow cover is regulated by seasons, topography, and hydrometeorology (Singh et al., 2014;Gurung et al., 2017;Rathore et al., 2018). Different techniques were proposed to estimate snow cover over different regions. Some of these techniques were adopted for snow cover monitoring over the Himalaya, and they were reported in Sood et al. (2020): normalized difference snow index (NDSI) Rathore et al., 2018), change detection , pan-sharpening , and snow cover mapping using snow depth maps (Gusain et al., 2016). Other than these techniques, snow cover mapping with interferometry coherence analysis using synthetic aperture radar (SAR) data (Kumar and Venkataraman, 2011) and, backscatter ratio-based technique (Thakur et al., 2013) have also been used over different sub-basins in the Himalayan region. The backscatter ratio from wet snow was calculated with the help of reference image from dry snow season. Different threshold values for different landcover values are reported for mapping the wet snow (Thakur et al., 2013). RISAT-1 data (∼50 m resolution) in HV (vertical transmit and horizontal receive) polarization is used for deriving timeseries of snow cover area (SCA) maps over Beas and Bhagirathi river basins during 2013-2014 (Thakur et al., 2017). The accuracy of the SCA was 95% when compared with the maps prepared using optical data. These studies (Kumar and Venkataraman, 2011;Thakur et al., 2013Thakur et al., , 2017 have demonstrated the potential of snow cover mapping using active microwave in Himalayan region. However, majority of the studies have focused mainly on using the optical data in the Indian Himalayan region for mapping snow cover and studying its variability Ahmad et al., 2020). The Himalayan snow cover does not have a uniform pattern or trend. In the upper parts of the Indus, the Ganga and the Brahmaputra basins, the snow cover area varies from 85 to 10% of the total area during ablation season (Singh et al., 2014). There is no trend in the mean snow cover over the Himalaya, however, at a regional scale, there is a statistically significant declining trend over small basins of Jhelum, Kosi, Manas, Gandaki, and Chandra (Gurung et al., 2017;Sahu and Gupta, 2020). For instance, in the lowelevation Ravi Basin, the snow cover decreases in mid-winter from 90 to 55 %, however, in the high-elevation Bhaga Basin, snow cover does not subside until April (Kulkarni et al., 2010). There is a shift in snow cover trends after 2010 and a deceleration in snow/ice cover shrinkage during recent years in the north-west Himalaya . However, due to large inter-annual fluctuation, snow depth depletion is not considered in snow cover variability, especially when the whole Himalayan range is examined (Gurung et al., 2011;Rathore et al., 2018).
The present synthesis shows that procuring a trend in snow cover changes is challenging for the Himalayan region due to the limited number of observations, substantial annual and inter-annual variability, and scarce in-situ snow depth area data . However, overall, it is widely reported that depleting snow cover and its spatiotemporal changes are associated with changes in climatic variables over the Himalaya Desinayak et al., 2022). Understanding variability of snow cover with respect to altitude and air temperature is vital for assessing availability of regional water resources and the impact of climate change in the Himalaya (Immerzeel et al., 2010;Shrestha et al., 2015). Moreover, the snow cover variation patterns, e.g., accumulation and ablation can potentially lead to imbalance in glacier mass balance and can trigger changes in water availability in different river basins (Kaser et al., 2010;Bolch et al., 2012).
Along with snow cover, snow depth (SD) and snow water equivalent (SWE) are also equally important in understanding regional hydrological cycle, and quantifying water availability in the Himalayan region. Annual snow water storage trends observed during 1987-2009 using Scanning Multichannel Microwave Radiometer (SMMR) data present a negative trend over Himalayan region (Smith and Bookhagen, 2018). The major contribution to SWE comes from mid-elevation range, i.e., 4,000-5,000 m a.s.l. in the Himalayan region and strongly affects the meltwater discharge. SWE in this elevation range has been shown to exhibit a strong negative trend (Smith and Bookhagen, 2018) along with shift in seasonality (Panday et al., 2011;Smith et al., 2017). Though exact mechanism behind these changes are not completely understood yet, few studies emphasized that aerosol contamination (Lau et al., 2010), changing precipitation patterns (Lutz et al., 2014), changes in western disturbances (Cannon et al., 2015), and evolving temperature dynamics could be the potential drivers behind SWE changes. Negative SWE trend can potentially impact the downstream water availability in the Himalayan region.
Monitoring of SD and SWE depends on the availability of a dense gauge network, but what we have is rather poor due to operational challenges. Remote sensing datasets, particularly microwave datasets [e.g., Advanced Scanning Microwave Radiometer for Earth observation (AMSR-E), Microwave Emission Model for Layered Snowpacks (MEMLS)] remains one of the most reliable techniques for obtaining SD and SWE in the Himalayan region. Various approaches consisting of non-linear models , physically based models (Liu et al., 2014), data driven approaches (Xiao et al., 2018;Yang et al., 2021) and assimilation-based approaches (Stigter et al., 2017;Kwon et al., 2019) has been applied in different studies for estimating SD and SWE. However, correcting the radar penetration issue is one of the biggest challenges in such remote sensing datasets (Tiwari et al., 2016;Patil et al., 2020;Awasthi et al., 2021).

ROCK GLACIERS IN THE HIMALAYA
Rock glaciers are lobate or tongue-shaped assemblages of icerich angular debris that moves or slowly creeps downslope due to gravity (Owen and England, 1998;Haeberli et al., 2006;Janke and Bolch, 2022). The surface velocity of active rock glaciers varies between 0.1 m to a few meters per year. Rock glaciers contain significant amounts of ice and are climatically more resistant than glaciers due to the thick insulating debris (Jones et al., 2019a;Janke and Bolch, 2022). Even though the origin of the ice is debated they are commonly considered as permafrost landforms (Berthling, 2011). Rock glaciers are gaining more attention because of their importance in the mountain water cycle in particular with the expected reduction of the glacier water resources Janke and Bolch, 2022).
Rock glaciers are also an important component of the debristransport system and their characteristics and occurrence is strongly influenced by the debris supply (Barsch and Jakob, 1998;Janke and Bolch, 2022). Rock glaciers can potentially transition from debris-covered glaciers especially in permafrost environments (Haeberli et al., 2006;Monnier and Kinnard, 2017;Jones et al., 2019b). Characteristics and evolution of debriscovered glaciers and their transition from clean -ice glaciers are also determined by debris supply. Hence, knowledge about the evolution of mountain slopes and geomorphological processes (such as valley-ridge scale interactions including rock slope failure and degradation of lateral moraines) that operate glacial and periglacial systems and how they impact the glaciers and rock glaciers are crucial. However, the ice-debris systems vary as the systems change in response to climate forcing. As a result, viewed from the land system perspective, a debris-ice land system incorporates numerous processes that respond to the climate in different ways over time.
Climate model projections resolve fundamental climatic and terrestrial processes only at broad scales; often around 100 km or so. Even downscaled products such as CORDEX-SA (South Asia) produce projections at around 25 × 25 km grids (e.g., Tangang et al., 2020). As a result, such projections at these spatial scales are challenging in resolving many of the geomorphological processes that operate in glacial systems, and which might affect their response to climate forcing.
For land surface models to enhance climate model projections of glacier mass balance and subsequent water supplies, they will need to more accurately account for paraglacial processes and be able to more accurately project the future evolution of debris-covered glaciers and rock glaciers to assess which glaciers will undergo the transition to debris-covered glaciers and rock glaciers, and which glaciers will not. The implications of such processes for future water supplies are likely profound. Recent research  has identified ∼25,000 rock glaciers across the Himalaya. Most are located within the central Himalaya (∼40%, n = 10,060), with ∼30% situated in the eastern Himalaya and ∼29% in the western Himalaya. The estimated areal coverage of rock glaciers across the area is 3,747 km² (i.e., intact and relict), representing ∼16% of that covered by glaciers (22,829 km²). Water volumes are estimated as 51.80 ± 10.36 km 3 , which translates to a rock glacier: glacier WVEQ (water volume equivalent) ratio in the central Himalaya of 1:17, with lower ratios in other regions. However, in some regions (e.g., in western Nepal), rock glacier WVEQ ratios reduce to between 1:3 and 1:5, suggesting strong regional contrasts in their significance. Also, rock glaciers are being mapped regionally, e.g., Bolch et al. (2022) identified 370 rock glaciers which cover an area of more than 10% of the total glacier area in the Poiqu basin in the Central Himalaya in Tibet. Pandey (2019) estimated 516 rock glaciers covering ∼350 km 2 area in Himachal Pradesh (northern India). Other forms of buried ice (i.e., ice-rich permafrost, ice-cored moraines, can also contain a significant amount of ice (Bolch et al., 2019b). No estimates have been made of the amount of ice in other non-glacial landforms, but it has been argued that this is likely to be substantial and should be accounted for better understanding the water resources in the Himalaya Janke and Bolch, 2022).

PERMAFROST IN THE HIMALAYA AND SURROUNDING AREAS
Permafrost is commonly defined as rock or soil with a temperature below 0 • C for more than two consecutive years. Permafrost ground often contains ice lenses or ice located in pores of rock and soil. Hence, permafrost impacts hydrology but also slope stability. Thawing permafrost can destabilize mountain slopes leading to hazards and impacting geomorphological processes (Krautblatter et al., 2013;Gruber et al., 2017;Bolch et al., 2019a).
In the Himalaya, most of the existing studies about the permafrost are from remote sensing methods. In-situ studies confirming the permafrost presence in the Himalaya are still rare (Gruber et al., 2017;Wani et al., 2020). A surface temperature model suggests that the total permafrost area in the HKH (including the Hindu-Kush and Tibetan Plateau) is 16 times higher than the HKH glacierised area (Gruber, 2012;Gruber et al., 2017). Another study suggested that rock glaciers can be used as a proxy to infer the existence of permafrost (Schmid et al., 2015).
Cryospheric components -including snow, glacier and permafrost-release meltwater from seasonal to multi-century frozen storage and control the spatiotemporal changes in the river runoff. Therefore, in a warming climate, permafrost thawing is expected to have similar impacts on river runoff like snow and glaciers (Hewitt, 2014;Biskaborn et al., 2019;Wang C. et al., 2019). However, unlike glaciers, detailed studies on ice rich permafrost are not yet available in the Himalayan region.
As a result, permafrost is not yet included in the hydrological modeling framework, and the impacts of permafrost thaw on the hydrological cycle remain unknown in the Himalayan region .
In some other parts of the Northern Hemisphere including the neighboring Tibetan Plateau, studies suggested permafrost degradation under recent climate change has resulted in hydrological changes, destabilization of rock glaciers, and frequent landslides (Zhang and Wu, 2012;Kargel et al., 2016;Rogger et al., 2017). These studies underscore the importance of including permafrost while designing frameworks for monitoring and assessing impact of climate change over catchments in the Himalaya . Understanding the volumetric and seasonal shifts in river runoff due to contribution of permafrost thaw would require a precise assessment of permafrost distribution, volume, and freeze/thaw processes . To fill these gaps, techniques should be developed that assess permafrost mass changes at large-scale by integrating time-varying changes of gravity with the GRACE satellites, digital elevation model differencing from multispectral images, and InSAR; while catchment-scale in-situ measurements can be used to validate satellite measurements . These in-situ data may include those from the automatic weather stations (Wani et al., 2020), geophysical surveys measuring electrical resistivity, data from ground penetrating radar, and temperature measurements from boreholes (Sjöberg et al., 2015;Wang C. et al., 2019).

DATA INTEGRATION AND ASSIMILATION
Combining different datasets can help in improving the data quality and also in determining an unknown variable when the others are known. It has been shown that the effective resolution of a dataset can be improved by assimilating information at better spatial resolution (Reichle et al., 2001;Houborg et al., 2012;Peng et al., 2017;Miro and Famiglietti, 2018). Various data assimilation techniques have been devised and successfully implemented, such as improving operational weather forecast, predicting ocean dynamics, and modeling soil moisture content (Jackson et al., 1981;Peng et al., 2017) in streamflow estimates Ramsankaran, 2017, 2018). Recently several studies have assimilated hydro-geodetic data and hydrological models to estimate, calibrate, or validate hydrological flux variables. For example, modeling river runoff with the help of hydrogeodetic approaches (Tourian et al., 2013;Sneeuw et al., 2014), estimating catchment-scale water budget using a Kalman filter framework (Pan and Wood, 2006;Lorenz et al., 2015), and calibration or/and validation of hydrological model outputs using GRACE datasets Eicker et al., 2014). The ensemble Kalman filter approach has been used effectively to assimilate GRACE TWSC into a Land Surface Model (LSM) to improve model performance (Zaitchik et al., 2008;Houborg et al., 2012;Eicker et al., 2014). Several non-parametric methods have also been proposed to improve spatiotemporal knowledge of hydrological variables. For example, Sun (2013) predicted groundwater level changes by incorporating GRACE with hydrometeorological variables in an ANN framework. Long et al. (2014) demonstrated the efficacy of ANN to predict TWSC from precipitation, soil moisture, and temperature, and Seyoum and Milewski (2017) used ANN to produce high-resolution TWS estimates. Recently Vishwakarma et al. (2021) used multivariate regression to improve the resolution of GRACE satellite products to half degree grids. Such higher spatial resolution datasets are required to better attribute spatiotemporal changes in water mass redistribution.
Data integration can also be carried out by using budget equations that rely on the conservation of mass or energy to define relations between various physical observables. These equations can be used to estimate one variable when others are known, and such a dataset is a derived dataset. For example, the water budget equation writes precipitation in a catchment as a sum of storage change, evapotranspiration, and runoff. If any three of these variables are known, then the fourth dataset can be estimated (Sneeuw et al., 2014;Lehmann et al., 2022). Such approaches have been used to estimate evapotranspiration, which, when compared to evapotranspiration from models, can provide an assessment of the accuracy of the model. It is to be noted that the efficacy of such an approach is limited by the accuracy of the variable with the highest uncertainty.

MAJOR GAPS
The Himalaya has the largest glacierised area outside poles, and therefore, it requires a massive collaborative effort to measure a representative number of glaciers in the field. Working on Himalayan glaciers is challenging, and only about 36 glaciers have been studied so far for mass balance (Table 1), which constitute about 1% of the total glacierised area in the Himalaya (Bolch et al., 2012;Azam et al., 2018). Therefore, the biggest gap is the lack of glacier in-situ mass balance data from the region. Likewise, our current knowledge of the volume of water stored and ice thickness distributions in individual glaciers is also poor (c.f. Figure 4). Hence, one of the fundamental research questions in front of us is "How are climate change and glacier volume responses affecting the glacier movements of Himalayan glaciers?" So far, this is still under investigation and has been recently highlighted by Azam et al. (2021). In Figure 4 we show a schematic that dissects our current and required knowledge in the spatial and temporal domains. We are yet to obtain a comprehensive understanding of frozen water at various spatiotemporal scales.
Another gap is the lack of in-situ meteorological datasets that are required to model glacier mass balance accurately. Temperature data from several reanalysis, gridded and model datasets (ERA5, ERA5-Land, APHRODITE, CRU, HAR, among others) are the best available options (Kumar et al., 2015Kanda et al., 2020). However, precipitation values from these sources are either over or underestimated due to the unpredictability of precipitation in high altitude regions (e.g., Immerzeel et al., 2015;Wortmann et al., 2018). Also, information from the inaccessible areas and avalanche contribution to glacier mass balance is typically ignored (Laha et al., 2017).
The downstream countries would be affected the most by changes in the Himalayan water storage. There are several hydro-power projects being developed that would be severely affected by changes in water availability. Recently several waterrelated disasters were reported that were triggered by sudden release of large water volumes, e.g., the Chamoli disaster in 2021 (Shugar et al., 2021) and the Kedarnath disaster in 2013 (Dobhal et al., 2013b;Allen et al., 2016). Although the communities in the downstream countries recognize the need for improving the quantity/quality of the observation data, but still there is a lack of cooperation between surface water and groundwater management agencies (including trans-boundary, inter-and intra-government agencies). In such scenarios, data from remote sensing platforms (e.g., Altimeters, GRACE) can be used to understand the state of water use. However, due to poor spatiotemporal resolution of remote sensing data (c.f. Figure 4), there is a need for exploring data-integration or assimilation methods that may predict hydrometeorological variables of interest at various scales. It is also recommended to develop a framework for upscaling the use of crowd sourcing-based data and remote sensing data, along with observation data from government agencies to formulate management plans.

CHALLENGES AND FUTURE RESEARCH DIRECTIONS
We have made significant advances in our understanding of water resources variability in the Himalaya and downstream during the past few decades. These advances resulted from an improved understanding of the physical processes and dynamics and the advances in data collection and automation. However, we still need to overcome some of the challenges to direct future research.
Installing more weather stations and collecting longterm continuous in-situ data at different elevation over the Himalayan regions for a comprehensive spatiotemporal coverage is needed. Instrumentation and data collection schemes should be developed to monitor changes in glacier mass, snow hydrological processes and permafrost storage. In-situ meteorological observations are sub-daily, while glaciological measurements are usually seasonal or annual, but both are at point scale. Remote sensing can help us obtain spatially consistent, sub-monthly or monthly time-series data. Hence, there is an eminent need to develop frameworks that can assimilate/integrate both in-situ and multi-satellite remote sensing observations to create high resolution time-series estimating various hydro-glaciological processes.
Models are an invaluable tool for predicting future water availability and understanding the role of individual drivers. Current climate models lack in representing feedback processes involving human interferences, snow and glacier processes, and the corresponding basin runoff responses. For example, at the regional scale, coupled glacier-climate model setups such as REMO glacier have already shown good potential in replicating the observed glacier mass-balance variability. Furthermore, estimation of uncertainties in model outputs must be improved. Ensemble modeling, for example, has been shown to provide a realistic assessment of uncertainty range.
Geodetic observations can provide critical data for monitoring snow, glaciers, surface, and sub-surface water storage. The FIGURE 4 | Spatiotemporal timeline of research investigation and understanding of the cryosphere and regional water resources in the Himalayan region. Relative data availability from various techniques/methods are also shown.
low-cost GNSS receivers can be deployed as they provide a wide variety of information on snow and glaciers. Particularly, a dense network of GNSS stations that can also perform GNSS reflectometry will be of immense value for snow and glacier monitoring in the Himalaya. The development of micro-electromechanical-system (MEMS) based gravimeters (Carbone et al., 2020) should make in-situ gravity measurements cheaper. It must be noted that terrestrial gravimetry is the most underutilized of all the geodetic techniques in glacier monitoring. Dedicated modeling of GRACE-FO data for retrieving Himalayan glacier mass change must be pursued to extract data for nearly two-decades. With the success of GRACE-FO, a number of future gravimetric missions are being proposed to ensure the continuous availability of temporal gravity data. Thus, investment in dedicated and integrated geodetic observations will open new vistas in the Himalayan glacier research.
Tracking changes in glacier mass balance, ice volume, permafrost storage, and snow hydrology over the Himalaya is challenging and we can minimize the above-mentioned gaps by capacity building and collaborations. Indian subcontinent hosts a large number of universities and institutions. Development of curriculum in this area of expertise in many of these institutions will increase the chances of minimizing the gaps. Also, capacity building from the secondary education level will prepare the upcoming generation to contribute significantly to this field. Collaboration and easy data sharing between institutions like India Meteorological Department (IMD), Indian Space Research Organization (ISRO), Defense Geoinformatics Research Establishment (Formerly Snow Avalanche and Study Establishment; SASE), Defense Research and Development Organization (DRDO), National Center for Polar and Oceanic Research (NCPOR), and other regional intergovernmental institutions such as the International Center for Integrated Mountain Development (ICIMOD) will also help immensely.

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
BDV, RR, and JB procured the funding for the workshop and organized it. Everyone was involved in discussions, contributed text relevant to their expertise and helped in revising the document. BDV, RR, MFA, and AM collated these contributions and continuously improved the manuscript. AM and SS prepared the