Modelling rock glacier velocity and ice content, Khumbu and Lhotse Valleys, Nepal

Rock glaciers contain significant amounts of ground ice and serve as important freshwater resources as mountain 10 glaciers melt in response to climate warming. However, current knowledge about ice content in rock glaciers has been acquired mainly from in situ investigations in limited study areas, which hinders a comprehensive understanding of ice storage in rock glaciers situated in remote mountains over local to regional scales. In this study, we develop an empirical rheological model to infer ice content of rock glaciers using readily available input data, including rock glacier planar shape, surface slope angle, active layer thickness, and surface creep rate. The model is calibrated and validated using observational data from the Chilean 15 Andes and Swiss Alps. We apply the model to infer the ice content of five rock glaciers in Khumbu and Lhotse Valleys, northeastern Nepal. The velocity constraints applied to the model are derived from Interferometric Synthetic Aperture Radar (InSAR) measurements. The inferred volumetric ice fraction in Khumbu and Lhotse Valleys ranges from 71% to 75.3% and water volume equivalents lie between 1.40 to 5.92 million m3 for individual landforms. Considering previous mapping results and extrapolating from our findings to the entire Nepalese Himalaya, the total amount of water stored in rock glaciers could 20 be in the magnitude of 10 billion m3, equivalent to a ratio of 1:17 between rock glacier and glacier reservoirs. Due to the accessibility of the input parameters of the model developed in this study, it is promising to apply the approach to permafrost regions where previous information about ice content of rock glaciers is lacking, and ultimately to estimate the water storage potential of the remotely located rock glaciers.

research, was first noted by Corte (1976); despite this, research on the role of rock glaciers in maintaining hydrological stores 30 in mountain catchments remains limited.
In regions such as the Himalaya, recent research has argued that rock glaciers might represent the end member of an evolutionary process where some glaciers transition to debris-covered glaciers, a proportion of which will then undergo further transition to rock glaciers Knight et al., 2019). This process would be triggered by the paraglacial response of high mountain slopes as glaciers undergo downwasting, which produces rock slope failures and mountainside collapses, 35 increasing the flux of rock debris to glacier surfaces. Depending on the debris cover thickness, this would be expected to limit ice melting and increase the resilience of glaciers to climate change (e.g., Reznichenko et al., 2010).
Recent work (Jones et al., 2021) was the first to show that around 25,000 rock glaciers exist in the Himalayas, covering 3747 km² and containing 51.80 ± 10.36 km³ of water volume equivalent. The comparative importance of rock glacier ice content versus that in glaciers in the region was 1:25, ranging from 1:42 to 1:17 in the Eastern and Central Himalaya and falling to 1:9 40 in Nepal. Importantly, we expect these existing ratios to reduce significantly as glaciers melt and undergo transitions to rock glaciers; yet the rates of transition from glacier to rock glacier are not understood. We also expect rock glaciers to provide water supplies long after glaciers have melted; in other high arid mountains, such as the Andes, ice-cored rock glaciers have persisted in valleys long after glacier recession (Azócar and Brenning, 2010;Monnier and Kinnard, 2015a). However, there lacks modelling studies to test these postulations and to assess the likelihood of glacier-rock glacier transition and the 45 hydrological implications of this process. A significant gap in our understanding of the likely future hydrological role of rock glaciers in arid mountains is the absence of quantitative information concerning their ice content. Currently, estimates of ice content in rock glaciers have focused on empirical information from drilling cores and boreholes (Hausmann et al., 2007;Kinnard, 2013, 2015a, b;Fukui et al., 2007;Arenson et al., 2002;Berthling et al., 2000;Croce and Milana, 2002;Florentine et al., 2014;Fukui et al., 2008;50 Guglielmin et al., 2004;Guglielmin et al., 2018;Haeberli et al., 1998;Haeberli et al., 1999;Krainer et al., 2015;Leopold et al., 2011;Steig et al., 1998), and from geophysical surveys (e.g., for reviews see: Hauck, 2013;Kneisel et al., 2008;Scott et al., 1990). However, these approaches are costly, time-consuming and extremely difficult to apply to rock glaciers at high altitudes and in remote mountains. It is therefore desirable to develop alternative approaches to understanding the likely ice content of rock glaciers, especially for regional scale estimates. 55 Ice content is one factor controlling the movement of rock glaciers by influencing the driving force and the rheological properties of materials which constitute the permafrost core (Arenson and Springman, 2005a;Cicoira et al., 2020), thus it is feasible to infer ice content using rheological modelling and observed kinematic data. Here we adapt an empirical rheological model by integrating rheological properties of rock glaciers derived from laboratory experiments (Arenson and Springman, 2005a), and parameterise the rheological model based on the structure and composition data of Las Liebres rock glacier 60 (Monnier and Kinnard, 2015b;Monnier and Kinnard, 2016). We then apply the model to simulate surface velocities of four rock glaciers with known ice content in the Swiss Alps and evaluate the modelling results to determine a suitable parameterisation scheme. Finally, we present results from modelling the kinematic responses of the coherently moving part of https://doi.org/10.5194/tc-2021-110 Preprint. Discussion started: 14 July 2021 c Author(s) 2021. CC BY 4.0 License. five rock glaciers in the study region of north-eastern Nepal and assess the modelled movement as a proxy for ice content by using remote sensing-derived downslope velocities as constraints. 65

Study area
Our study area comprises the Khumbu and Lhotse valleys in north-eastern Nepal (Fig. 1a). The glaciers draining Everest and Lhotse (e.g., Khumbu and Lhotse glaciers) are the highest in the world and have well defined debris-covered snouts. The tributary valleys contain a variety of rock glaciers and composite landforms where glaciers are transitioning to rock glaciers 70 Knight et al., 2019). There are five rock glaciers in the study area, namely Kala-Patthar, Kongma, Lingten, Nuptse, and Tobuche (Fig. 1b). The five rock glaciers examined in this study are situated at 4900-5090 m a.s.l., near the altitudinal boundary of discontinuous permafrost in the region. Previous seismic refraction surveys conducted on active rock glaciers indicate that the lower limit of permafrost occurrence in this region to be ≈5000-5300 m a.s.l. (Jakob, 1992), which is consistent with an earlier estimate of 4900 m a.s.l. based on ground temperature measurements (Fujii and Higuchi, 1976). 75 Meteorological data provided by the Pyramid Observatory Laboratory near Lobuche village on the western side of the Khumbu Glacier (5050 m a.s.l.) reveal that the dominating climate of this area is the South Asian Summer Monsoon. For the period of 1994-2013, recorded accumulated annual precipitation is 449 mm yr -1 , with 90% of the precipitation concentrated during June-September (Salerno et al., 2015). The mean annual air temperature is -2.4 °C (Salerno et al., 2015).

Methods
Our methods are divided into two parts and detailed in the subsections below. First, we derived surface kinematics of rock 85 glaciers using Interferometric Synthetic Aperture Radar (InSAR; Sect. 3.1). Second, we developed a rheological model for estimating ice content of rock glaciers by using surface velocity as a constraint. Sect. 3.2 describes the model design, calibration, and validation procedures, as well as application of the model.

Deriving surface kinematics with InSAR 90
Nineteen L-band ALOS PALSAR images and twenty-one ALOS-2 PALSAR-2 images acquired during 2006-2010 and 2015-2020, respectively, were used to form more than fifty interferograms to measure the surface displacements of the landforms in the study area (Table 1). We selected SAR data to achieve high interferometric coherence by following the criteria such as: (1) short temporal spans (less than 92 days for ALOS pairs and 70 days for ALOS-2 pairs); (2) short perpendicular baselines (smaller than 800 m for ALOS pairs and 400 m for ALOS-2 pairs). We estimated and removed the topographic phase with the 95 1-arcsec digital elevation models (DEM) produced by the Shuttle Radar Topography Mission (SRTM) (spatial resolution ~30 m). Multi-looking operation and adaptive Goldstein filter (8×8 pixels) were applied to the interferometric processing, which was implemented by the open-source software ISCE version 2.4.2 (available at https://github.com/isce-framework/isce2). The interferograms were unwrapped using the SNAPHU software (Chen and Zebker, 2002). We randomly selected three pixels located at flat and stable ground near each ice-debris landform and averaged their phase values to re-reference the unwrapped 100 phases measured within the landforms. By doing so, atmospheric artefacts including the water vapour delay and ionospheric effects can be effectively removed because these are spatially long-wavelength features and can be assumed as constant within the range of our study objects (Hanssen, 2001).
We then derived the surface velocities along the SAR satellite line-of-sight (LOS) direction from the unwrapped interferograms and projected the LOS velocities to the downslope direction of the landforms. Uncertainties were quantified by considering 105 the error propagation of the InSAR measurements and associated geometry parameters (Hu et al., 2021).
After that, to ensure high data quality, we selected the InSAR observations meeting the following criteria as valid results for further analyses: (1) the pixels showing low coherence (<0.3) are masked out before velocity statistics, and more than 40% of pixels remain after the masking procedure; (2) the mean velocity of the landform is larger than 5 cm yr -1 .
Next, we defined and outlined the coherently moving part of the landform by considering the time series of downslope velocity 110 of each pixel acquired during all the observational periods. If the InSAR-measured velocity is higher than 5 cm yr -1 in more than half of the periods at a given pixel, it was included into the coherently moving part of the landform.
Finally, we analysed the velocity values of all pixels within the coherently moving part of the landform and selected the mean, median, and maximum values for each observation to characterise the surface kinematics of the landforms.

Estimating ice content from a surface-velocity-constrained model
This subsection describes the process of model development. The workflow is illustrated in Fig. 2.

Model design and assumptions
Active rock glaciers creep as a result of two internal processes: plastic deformation of the ice-rich permafrost core, and deformation at the shear horizon at depth (e.g., Arenson et al., 2002;Berthling, 2011;Cicoira et al., 2019b;Haeberli, 2000;Kenner et al., 2019). Many previous modelling studies depict the deformation mechanism of rock glaciers based on Glen's flow law (e.g., Arenson and Springman, 2005a;Cicoira et al., 2020;Whalley and Azizi, 1994), which essentially relates strain 125 rate () with effective shear stress ( ) and describes the rheology of ice flow (Glen, 1955): where and are parameters reflecting variations in environmental conditions (mainly including temperature and pressure), material properties (such as composition, structure, and texture), and operating creep mechanisms (e.g., diffusion and dislocation). 130 In this study, we primarily adopted a creep model of ice-debris mixture, proposed by Moore (2014), based on Glen's flow law: where is a strain enhancement factor; is a parameter reflecting the strength of the ice-debris mixture, associated with the volumetric debris content ( $ ). When $ is less than a critical volumetric debris content ( $% ), the strength of the mixture is governed by interparticle friction, and the value of equals one. Theoretically, $% is around 0.52 (Moore, 2014). "# is a 135 threshold stress imparted by the frictional strength between debris particles, also depending upon the volumetric debris content ( $ ). Assuming that "# ≪ , $ < $% , and = 1, Eq. 2 can be reduced to the following form (Monnier and Kinnard, 2016): where is the effective viscosity and is equal to 3 Following a common setup in glaciology (Cuffey and Paterson, 2010), we consider each rock glacier as a slab with uniform width and thickness and a semi-elliptical cross-section, resting on a bed of constant slope. It consists of two layers: an active layer and a permafrost core. The active layer is a mixture of debris and air, and the permafrost core consists of ice, water, debris and air. Both layers are assumed as homogeneous. Movement of rock glaciers is caused by the steady creep of the permafrost core in the plane parallel to the bed slope. The active layer moves passively along with the inner core, which has 145 been validated by observations (Arenson et al., 2002;Haeberli, 2000). From Eq. 3 and the structure and geometry illustrated in Fig. 3, we have: where $, $-is the velocity derivative relative to the depth in the permafrost core. At a given depth , the driving stress is imparted, taking into account the loading of the above material and the effect of 150 frictional drag occurring between the lateral margins and surrounding bedrocks, which is represented by a shape factor . (Cuffey and Paterson, 2010): where is the slope angle; is the gravitational acceleration; /0 and %123 are the densities of the active layer and the permafrost core, respectively; ℎ /0 is the active layer thickness. 155 The shape factor is expressed as (Oerlemans, 2001): where and are the width and thickness of the rock glacier, respectively.
The integration of the velocity profile (Eq. 4 and 5) is expressed as: where ? is the surface velocity as illustrated in Fig. 3. When is set as the thickness of the ice core (ℎ %123 ) and basal sliding is assumed to be absent, ? is then expressed as: The densities of the active layer ( /0 ) and the permafrost core ( %123 ) are given as: where $,/0 and /,/0 are the volumetric contents of debris and air in the active layer, respectively. The volumetric contents of the components in the inner core, namely debris, air, ice and water, are expressed as $,%123 , /,%123 , F,%123 , and G,%123 , respectively. $ , / , F , and G are the densities of debris, air, ice, and water, respectively. 170 For the flow law exponent ( ), we first used an empirical average value as assumed in numerous glaciological studies: We also adopted a linear relationship between and the volumetric ice content ( F,%123 ) based on laboratory experiments undertaken on borehole samples from two rock glaciers (Arenson and Springman, 2005a):

Model calibration
Combining Eq. 9-11 with Eq. 12 or 13, we formulated several models depicting the relationship between the surface velocity and properties of rock glaciers, including their composition, structure, and geometry. We then calibrated the models to determine the curve of best fit between the effective viscosity ( ) and the volumetric ice content ( F,%123 ) using observational 185 data of Las Liebres rock glacier in Central Chilean Andes (collected by Monnier and Kinnard, 2015b). This dataset includes information of structure (ℎ %123 and ℎ /0 ), geometry ( and . ), and composition ( $,%123 , /,%123 , F,%123 , and G,%123 ), all of which were derived from Ground Penetrating Radar (GPR) measurements. Surface velocities ( ? ) were provided by a Differential Global Positioning System (DGPS) along the central creep line at 14 locations on Las Liebres rock glacier, (detailed in Kinnard, 2015b &2016). 190 First, we adopted the exponential -F,%123 relation estimated by Monnier & Kinnard (2016) with the same dataset and a constant creep parameter (Eq. 12) (Fig. 4a). Then by integrating the relationship between and ice content (Eq. 13), we applied both a 2 nd -degree polynomial regression model and an exponential regression model to determine the -F,%123 relationship (Fig. 4b, c). The polynomial regression model is used to capture the subtle increase in effective viscosity when the ice fraction becomes larger at the ice-rich end. This trend is also depicted in the laboratory experiment conducted by 195 Arenson and Springman (2005a) as a parabolic relationship between the minimum axial creep strain rate and the volumetric ice content. Finally, we obtained three candidate parameterisation schemes expressed as: For simplicity, the parameterisation scheme proposed in Monnier & Kinnard (2016) is labelled as Scheme 1 (Eq. 14). The parameterisation scheme considering the empirical relation between and F,%123 (Eq. 13) and parameterisation scheme derived from the polynomial and exponential relationship between and F,%123 are marked as Scheme 2 and Scheme 3 (Eq. 15 and 16), respectively.

Model validation 210
The three parameterisation schemes (Eq. 14-16) were validated using observational data of four rock glaciers in the Swiss Alps, namely Murtèl-Corvatsch, Gruben, Muragl, and Schafberg (Cicoira et al., 2019a;Arenson et al., 2002;Hoelzle et al., 1998;Barsch et al., 1979). We simulated the surface velocity ( ? ) of each rock glacier by varying volumetric ice content ( F,%123 ) of the permafrost core and inferred its ice fraction by comparing the modelled velocity and the measured velocity from Terrestrial Geodetic Surveys (PERMOS, 2019). We then referred to the previously estimated ice content of the selected 215 rock glaciers to validate our predicted results.
To derive the input parameters, we first outlined the boundaries of the four rock glaciers, from which their shapes and areal extents can be extracted. An empirical relationship established by Brenning (2005b) was then applied to calculate the rock glacier thickness ( ) from its areal extent ( 29 ): where the area ( 29 ) is in km 2 . The width of each glacier was quantified as the width of its minimum envelop rectangle. We took the mean value of the active layer thickness obtained from borehole measurements in the PERMOS network as the input parameter ℎ /0 for each rock glacier. The surface slope ( ) was calculated based on the SRTM DEM with a spatial resolution of 30 m. Table 2 lists the values of the above parameters. The permafrost core thickness (ℎ %123 ) can be obtained by subtracting ℎ /0 from the total thickness calculated using Eq. 17. 225 We assumed the volumetric ice content ( F,%123 ) of the permafrost core to be between 40% and 100%, considering the prerequisites of the modified ice-debris mixture flow law (Eq. 3) that the debris fraction ( $,%123 ) should be less than the threshold ( $% ) (Sect. 3.2.1). We varied the ice content ( F,%123 ) by 1% in each step to model the corresponding surface velocities ( ? ). We fixed the air content in the permafrost core as 7.5%, which is a mean value of the air fraction in ice-rich permafrost samples (Arenson and Springman, 2005b). At near 0 °C, the volumetric content of water ( G,%123 ) displays a 230 positive correlation with the debris fraction ( $,%123 ) (Monnier and Kinnard, 2016). Thus, we calculated the $,%123 -G,%123 correlation based on the data published in Monnier and Kinnard (2015b) and assumed the constitution of the selected rock glaciers followed the same linear relationship (Fig. 4d). The debris density ( $ ) was given as 2450 kg/m 3 (Monnier and Kinnard, 2016). The density of air ( / ) is determined by the elevation of each rock glacier. The ice density ( F ) is 916 kg/m 3 and the water density ( G ) is 1000 kg/m 3 . 235

Sensitivity analysis
To explore how uncertainties of the input parameters contribute to the final output of the developed approach, we tested the response of the model to varying input parameters by performing a series of synthetic sensitivity experiments. For these experiments, we simulated surface velocities of the rock glacier with varying ice fractions and inferred the current ice content from the velocity constraint. The parameters explored here are detailed in Table 3. A reference scenario is set up with the parameters of Murtèl-Corvatsch rock glacier and labelled as Sc-1.0. We designed eight scenarios extending from Sc-1.0, naming each scenario after a multiplication factor which indicates the ratio between the parameter in each scenario and that in the reference scenario; with the exception of two parameters, namely debris density (ρ Q ) and debris fraction in the active layer (θ Q,RS ), where we changed the upper or lower boundary of the value range to be consistent with the usual value range in reality. 245 We performed the sensitivity experiments by varying one parameter at a time while keeping the other variables constant. Table 3. Parameters of the sensitivity experiments. Scn-1.0 is the reference scenario that adopts the parameters of Murtèl-Corvatsch rock glacier. The other scenarios are designed by multiplying the reference value of each variable with the corresponding factor in their scenario labels.

Model application 250
We applied the validated model with the optimal parameterisation scheme to infer ice content of the coherently moving part of five rock glaciers in the Khumbu and Lhotse Valleys. The geometric and structural data used as input parameters are detailed in Table 4. Area, width, and slope angle are quantified using the same method as described in Sect. 3.2.3. Active layer thickness was determined as the mean value over the extent of each rock glacier during 2006-2017 from the European Space Agency Permafrost Climate Change Initiative Product (ESA CCI) (Obu et al., 2020). The same empirical relation for calculating rock 255 glacier thickness as used in the validation procedure was adopted here to obtain the thickness parameter. The surface velocity constraint is the range of InSAR-derived downslope velocity during the observed period (Sect. 3.2.2); except for Tobuche RG where the abnormal value in 2015 is removed from the range (see Sect. 4.1 for details). Finally, we calculated the water volume equivalent to estimate the amount of water stored in rock glaciers by considering the inferred ice content, areal extent, and permafrost core thickness. 260

Results
In this section we first summarise the surface kinematic characteristics of rock glaciers in Khumbu and Lhotse Valleys measured by InSAR. Then we show the results of model validation and sensitivity experiments. Finally, we present the inferred ice contents and estimated water storage of rock glaciers in the study area.

InSAR-derived surface kinematics of rock glaciers 270
We used InSAR to derive the downslope surface velocities of five rock glaciers situated in the study region. Surface kinematics of debris-covered glaciers were also quantified and presented in Fig. S1 and S2 in the supplementary materials. Figure 5 shows the time series of the InSAR-derived surface velocities of the rock glacier coherently moving parts. We observe that the median and mean velocities of each landform have similar values, and both are capable of characterising the kinematic status of the landforms. By selecting the mean velocity as the representative value, most rock glaciers, except for Tobuche, 275 moved at a nearly stable rate, ranging from 5 cm yr -1 to 30 cm yr -1 during the observational period, with the largest standard deviation being 3.4 cm yr -1 for Lingten. The maximum velocity represents the local extreme of downslope motion and was as high as 112.1±12.4 cm yr -1 for Lingten during 2019/07/15-2019/08/26. Tobuche displayed similar stable behaviour before 2010 but had accelerated by more than four times from 14.9±0.2 cm yr -1 to 81.4±2.4 cm yr -1 since 2010. The maximum velocity reached was 181.0±57.4 cm yr -1 for the period 2015/03/18-2015/03/22. However, the associated uncertainties during this 280 period were high: the relative uncertainties of mean, median, and maximum velocity were 2.9%, 38.2%, and 31.7%, respectively. Therefore, the acceleration of Tobuche cannot be confidently revealed by our data. The extents of coherently moving parts of the five rock glaciers are presented in Fig. 6, with the average velocities derived from the interferograms obtained during the past several years.

Model validation
With the input parameters and model setup as detailed in Sect. 3.2.3, we simulated the surface velocities ( ? ) of each rock glacier using Schemes 1-3. Uncertainties, generated through the statistical analysis used to establish the model (as shown in https://doi.org/10.5194/tc-2021-110 Preprint. Discussion started: 14 July 2021 c Author(s) 2021. CC BY 4.0 License. Fig. 4), have all been considered in the simulation. We used the annual mean surface velocities, calculated from the Terrestrial 295 Ground Survey data (PERMOS, 2019), as the constraint for inferring the ice content.
For each rock glacier, an inferred ice content range is derived based on the velocity constraint and modelled ? -F,%123 relationship. The median of the range is selected as the inferred ice content and compared with the reference ice content, which is taken as the average value of the estimated ice content based on previous field measurements (Cicoira et al., 2019a;Arenson et al., 2002;Hoelzle et al., 1998;Barsch et al., 1979). 300 Comparing the reference and inference ice content from the three schemes, Scheme 2 is the optimal one for the following two reasons: (1) the reference ice content is within the range inferred from Scheme 2 (Fig. 8); (2) Scheme 2 gives the smallest average bias (8.4%) compared with Scheme 1 (12.9%) and Scheme 3 (13.3%) (Table 5). However, the above bias is not statistically useful for correcting the modelling results due to the limited amount of validation data.

Model sensitivity
The results of sensitivity experiments are shown, normalised to the corresponding values of the reference scenario (Scn-1.0), in Fig. 10. We observe that the inference result remains stable in response to most varying parameters, with a bias of less than 5%, relative to the reference scenario (Scn-1.0). Surface slope angle influences the result most: in the extreme scenario (Scn-325 0.2), the inferred ice content can be altered by 15%. In non-extreme cases (e.g., Scn-0.8, Scn-0.6), the influences of varying slope angles can be well constrained within the 5% range. In general, the model is insensitive to the uncertainties of any single input parameter.  Figure 11 and Table 6 present the inference ice contents of rock glaciers in the study area. The inferred average ice fractions of the landforms are in the range of 71.0-75.3%; the water volume equivalent ranges from 1.40 to 5.92 million m 3 for individual 335 landforms. The maximum range of ice fraction is estimated to be 57.5-92.0%; the corresponding water volume equivalent ranges from 1.13 to 7.24 million m 3 . Nuptse stores the most ice by volume due to its largest dimensions. The total amount of water stored in each of the five rock glaciers lies between 10.61 and 16.54 million m 3 , with an average value of 13.57 million m 3 .

Discussion
The following subsections firstly investigate the potential water storage of rock glaciers over the Nepalese Himalaya by extrapolating the estimated water storage in the Khumbu and Lhotse Valleys (Sect. 5.1). We then discuss the validity of the 350 approach developed in this study for inferring ice content by analysing the general assumptions and design of the surfacevelocity-constrained model (Sect. 5.2), based upon which we discuss the application of the method to large-scale regions and further improvements (Sect. 5.3).

Potential water storage in rock glaciers in the Nepalese Himalaya
The inferred average ice content of the five rock glaciers in the study area lies within a narrow range (71.0-75.3%), mainly 355 due to their similar observed downslope velocities (5-30 cm yr -1 ), used as modelling constraints ( Fig. 5; Fig. 6; Sect. 4.1). In general, rock glaciers typically creep at a rate ranging from decimetre to several meters per year (Delaloye and Echelard, 2020), thus the average ice content of the five rock glaciers may not be able to represent the motion of all rock glaciers situated in the entire mountain range.
A previous study has compiled an inventory including 4226 intact rock glaciers over the Nepalese Himalaya, and a first-order 360 estimate has indicated that these landforms contain between 16.72 and 25.08 billion m 3 of water (Jones et al., 2018b). By extrapolating from the estimated water storage in the five rock glaciers found in this study, the total amount of water stored in all the intact rock glaciers ranges from 8.97 to 13.98 billion m 3 over the entire Nepalese Himalaya, which is in the same magnitude predicted by the previous research (Jones et al., 2018b). In the Nepalese Himalaya, the ratio between the amount of water stored in rock glaciers (11.47 billion m 3 ) and in glaciers (197.63 billion m 3 ) is 1:17. By using the minimum and 365 maximum inference values, the estimated ratio ranges between 1:21 and 1:13. Our modelling-based results are lower than earlier estimates (1:9), yet reveal higher hydrological importance than across the Himalayas (1:24) (Jones et al., 2018b;Jones et al., 2021).

Validity of general assumptions and model design
Our aim was to develop a rheological model that allows for inferring ice content of rock glaciers in areas where in situ 370 investigations are scarce. To achieve this objective, certain simplifications have been applied to the model setup. Here we discuss the validity of the method through the following six aspects of the model assumptions and design, including (1) assumption of a steady-state rock glacier creep, (2) homogeneous warm permafrost hypothesis, (3) neglect of shear horizon in the model design, (4) accuracy of rock glacier thickness derivation, (5) identification of the coherently moving part of rock glaciers, and (6) generalisation of statistical relationships derived from Las Liebres rock glacier. 375

Steady-state creep of rock glaciers
By using the adapted form of Glen's flow law (Eq. 2), we primarily assume the rock glacier movement to be steady-state creep driven by viscoelastic deformation of the ice-debris mixture (Moore, 2014). This premise indicates that our method is applicable to rock glaciers currently moving at a relatively stable rate. Recent research has reported abrupt and significant acceleration of rock glaciers triggered by abnormal surface warming events (Delaloye et al., 2013). These destabilised rock 380 glaciers are beyond the applicability of our method. In this study, we quantified surface kinematics of rock glaciers over multiple years to quantify the stability of the rock glacier motion. The seasonal variations in creep rate are neglected and sudden acceleration events are excluded in the velocity range used as the model constraint (Sect. 4.1).
We also exclude the component of basal sliding processes in our model design (Fig. 3), which operates at the base of some rock glaciers, as observed in the Tien Shan (Harrison unpubl.). This is not surprising; many clean glaciers undergo basal sliding 385 in certain situations (e.g., Vivian and Bocquet, 1973) although, a debris-rich layer with high water content at the base undergoes enhanced deformation and some sliding (e.g., Boulton and Jones, 1979;Echelmeyer and Wang, 1987). This process also occurs in rock glaciers with a high ice content and is accompanied by the disruption of sediment and vegetation at the front of such features. The kinematics of these rock glaciers cannot be appropriately simulated by our current approach.

Homogeneous warm permafrost 390
Ground temperature is one of the factors controlling the creep parameter ( ) (Eq. 1), as described by the Arrhenius relation (Mellor and Testa, 1969). As ground temperature changes with depth, primarily due to heat diffusion, creep parameter ( ) varies along the vertical profile of the rock glacier. Previous studies implemented different relationships between creep parameter and temperature, and integrated a heat diffusion model (proposed by Carslaw and Jaeger, 1959) to consider this effect (Azizi and Whalley, 1996;Arenson and Springman, 2005a;Kääb et al., 2007;Ladanyi, 2003). 395 In our model design, we use the effective viscosity ( ) to absorb the intricate effects of strain enhancement factor ( ), threshold stress ( "# ), and most importantly, the creep parameter ( ) (Eq. 3), which reduces the number of input parameters and allows for developing an empirical relationship between effective viscosity and ice content based on an existing observational dataset and laboratory findings (Eq. 13). The data and relationship between variables we used for model calibration are derived from observations of rock glaciers situated in a warm permafrost environment (>-3°C) (Monnier and Kinnard, 2016;Arenson and 400 Springman, 2005a) Measurements of ground temperature in the study area are scarce in general. However, we infer that these rock glaciers develop in a warm permafrost environment for the following reasons: (1) the landforms are located near or below the altitudinal limit of permafrost distribution in Nepal (Sect. 2) (Fujii and Higuchi, 1976;Jakob, 1992), indicating that the local environment is at the critical condition of permafrost occurrence; (2) based on empirical relationships between mean annual ground 405 temperature (MAGT), mean annual air temperature, latitude, and altitude, the estimated MAGT is >0.5°C, which suggests that permafrost in this area is in a warm and unstable condition (Nan et al., 2002;Zhao and Sheng, 2015).

Neglect of shear horizon
The shear horizon is discovered from borehole investigations and is defined as the thin layer situated at more than ten meters deep where the majority of internal deformation takes place (Arenson et al., 2002;Buchli et al., 2018;Haeberli et al., 1998). 410 Field observations and numerical modelling suggest that unfrozen water within the shear horizon plays an important role in controlling the seasonal variations in rock glacier creep (Buchli et al., 2018;Cicoira et al., 2019b;Kenner et al., 2019). This short-term feature of rock glacier kinematics is insignificant to modelling the relationship between ice content and multi-annual average movement velocity in our study.
Previous studies have considered the enhanced deformation occurring in the shear horizon additionally, but it requires detailed 415 knowledge of the internal structure, i.e., the depth of shear horizon (Frehner et al., 2015;Ladanyi, 2003). To tackle the issue of data insufficiency of internal rock glacier structure in this study area -as with most permafrost areas -we neglect the distinct rheology in the shear horizon and assume a constant effective viscosity instead. This simplification has also been adopted in other research aiming at studying rock glacier dynamics over a large-scale extent (Cicoira et al., 2020).

Accuracy of rock glacier thickness derivation 420
The accuracy of rock glacier thickness is discussed here because it influences the surface kinematics most significantly. As shown in Eq. 8, the surface velocity is proportional to the thickness to the power of + 1, resulting from the vertical integration of Eq. 7. We adopt the empirical relationship between rock glacier area and thickness (Eq. 17) based on field observations in the Andes (Brenning, 2005a). The derived thicknesses of the four rock glaciers in the Swiss Alps, used for validation, are consistent with previous in situ measurements. However, another rock glacier, namely Ritigraben, situated in the same region, 425 does not follow this empirical relationship and has a bias as large as ten meters compared with the field estimates. We also test a linear relationship between surface slope angle and thickness, recently established by Cicoira et al. (2020). The accuracies of the results turn out to be at the same level, with four out of five rock glaciers having good estimation results (Table 7). Thus, the uncertainty introduced by thickness derivation when applied to rock glaciers without known information of structure is unavoidable. 430

Identification of the coherently moving part of rock glaciers
Taking advantage of the multi-temporal observations and continuous spatial coverage of the InSAR measurements, we define 435 the coherently moving part of each rock glacier in our study area ( Fig. 6; Sect. 3.1) and infer the ice content of the coherently moving part using the developed modelling approach. We introduce this concept because it corresponds with the general model setup (Fig. 3). Moreover, with the assistance of displacement maps generated by InSAR, the defined boundary of coherently moving rock glacier avoids the ambiguities involved when delineating rock glaciers solely based on geomorphologic features from a highly dynamic environment where complex glacial, periglacial, and paraglacial processes take place (Jones et al.,440 2019; Delaloye and Echelard, 2020). In addition, defining the coherently moving part associates with the uncertainty in the rock glacier area, which is unlikely to affect the modelling result significantly due to the insensitive response of our model to this input parameter (Fig. 10). Finally, employing this kinematics-based definition, may also contribute to the currently inevitable uncertainties in thickness estimation (Sect. 5.2.4), as Eq. 17 (proposed by Brenning, 2005b) uses rock glacier area as a variable, without considering the landform activity. 445

Generalisation of statistical relationships derived from Las Liebres rock glacier
The model developed in this study essentially relies on statistical relationships, especially the effective viscosity-ice content relationship, derived from geophysical surveys conducted on Las Liebres rock glacier in the Andes (Monnier and Kinnard, 2015b). Applying the calibrated model to rock glaciers in other areas is primarily based on the assumption that a common flow law governs the rheology of rock glaciers developed in warm permafrost environment, irrespective of locality. In addition to 450 the general hypothesis, we tackled this issue in two ways. First, we validated the model using samples from a different region, i.e., the Swiss Alps (Sect. 3.2.3). Secondly, the uncertainties introduced by the statistical analysis have all been quantified in the validation and inference procedures, which leads to a wide range of the inference ice contents (Sect. 3.2.3 and 4.4).

Application of the model and outlook
By adopting the assumptions and setup as discussed in Sect. 5.2, we firstly present an approach for estimating ice content of 455 rock glaciers with simple input parameters, including the planar shape (area and width), surface slope angle, active layer thickness, and surface velocity, all of which can be obtained directly or derived from remote sensing techniques and products. Therefore, this established method for estimating the amount of ice stored in numerous rock glaciers with well-quantified uncertainties is ready to be applied to many remote alpine environments.
To improve the performance of the approach, more data obtained from field and geophysical investigations, especially detailed 460 data of rock glacier composition, can be integrated in the future to calibrate and validate the empirical rheological model. More reliable methods for estimating rock glacier thickness will also improve the accuracy of the modelling results.

Conclusions
We develop an empirical rheological model for inferring ice content of rock glaciers and apply it to estimate the water storage of rock glaciers situated in the Khumbu and Lhotse Valleys using surface-velocity-constraints derived from InSAR 465 measurements. The main findings are summarised as follows: (1) An empirical rheological model is presented in this study for estimating ice content of rock glaciers using five input parameters, namely rock glacier area, width, surface slope angle, active layer thickness, and surface velocity, all of which can be obtained from readily available remote sensing products or emerging datasets.
(3) The inferred average ice contents of rock glaciers in Khumbu and Lhotse Valleys are in the range of 71.0-75.3%; the water volume equivalent ranges from 1.40 to 5.92 million m 3 for individual landforms. Nuptse RG stores the most ice due to its largest dimensions among the five studied rock glaciers.
(4) The inference range of ice content of the landforms lies between 57.5 and 92.0%. Total amount of water stored in the five 475 rock glaciers in Khumbu and Lhotse Valleys ranges from 10.61 to 16.54 million m 3 , with an average value of 13.57 million m 3 .
(5) Considering previous estimates and extrapolating from our inference results in Khumbu and Lhotse Valleys, the total amount of water stored in rock glaciers over the Nepalese Himalaya is in the magnitude of 10 billion m 3 , and the ratio between water storage in rock glaciers and glaciers ranges from 1:13 to 1:21, averaging at 1:17. 480 This study demonstrates the effectiveness of inferring ice content of rock glaciers by using a surface-velocity-constrained model. The estimated ice content and water storage in the study area highlights the hydrological significance of rock glaciers in the Nepalese Himalaya. The approach developed here can be readily applied to other alpine permafrost regions where rock glaciers are widespread for a preliminary water resource evaluation.

Code and data availability 485
The source code of ISCE is available at https://github.com/isce-framework/isce2. The ALOS PALSAR and ALOS-2 PALSAR-2 data are copyrighted and provided by the Japan Aerospace Exploration Agency through the EO-RA2 project ER2A2N081.
Data for the rock glacier kinematics in the Swiss Alps are available at http://www.permos.ch/data.html. The ESA CCI permafrost data are available at http://catalogue.ceda.ac.uk/uuid/1f88068e86304b0fbd34456115b6606f. The code of the modelling approach for estimating ice content will be provided by Yan Hu upon request. 490

Author contribution
YH developed the code, performed the data analysis and interpretation, visualised the results, and wrote the majority of the manuscript. SH conceptualised the research goal, supervised the study, and wrote Sect. 1 of the draft. LL advised YH and actively helped the investigation process. JLW helped formulate the initial framework of the method and collect research data.
All the authors contributed to the reviewing and editing of the manuscript. 495