The influence of water percolation through crevasses on the thermal regime of Himalayan mountain glaciers

In cold and arid climate, small glaciers with cold accumulation zone are often thought to be entirely co d based. However, Ground Penetrating Radar (GPR) measurement s on Rikha Samba Glacier in the Nepal Himalaya reve al a large amount of temperate ice that seems to be influenced by the presence of crevassed areas. We used a coup led thermomechanical model forced by a firn model accounting for firn heating to interpret the observed thermal regime. We show that the addition of water percolation and refreezing in crevassed areas using a simple energy conservative approach is able to 20 explain the observations. Model experiment shows th at both steady and transient thermal regimes are si gnificantly affected by latent heat release in crevassed areas. It makes half of the glacier base temperate, resulting in t he dynamics mainly controlled by basal friction instead of ice deforma tion. Timescale of thermal regime change in respons e to warming climate is also greatly diminished with a potential switch from cold to temperate basal ice in 50-60 years in the upper part of the glacier while it would take 100-150 years without t he crevasses effect. This study highlights the cruc ial role of water 25 percolation through the crevasses on the thermal re gime of glaciers and validates a simple method to t ake it into account in glacier thermo-mechanical models. 30 https://doi.org/10.5194/tc-2019-172 Preprint. Discussion started: 14 August 2019 c © Author(s) 2019. CC BY 4.0 License.


Technical corrections
(grammatical suggestions annotated as "(gr)") 16: (gr) "In cold and arid climates" Done 18: How does GPR reveal temperate ice? Is it more correct to say that GPR suggests/implies temperate ice due to bed reflectivity measurements? We modified the sentence: "However, scattering in Ground Penetrating Radar (GPR) measurements on Rikha Samba Glacier in the Nepal Himalaya suggests a large amount of temperate ice that seems to be influenced by the presence of crevassed areas" (lines 17-18). Can any additional information be provided about the type of antenna and transmitter used? How wide is the frequency band over which energy is produced? Is this a frequency domain or impulse-type system? Description of the antenna and transmitter have been added to the text also clarifying that it is impulse-type system with a frequency band of .
73: Could you clarify whether all reflectors were picked, or only the one thought to represent the glacier base? The reflectors were picked from the assumed glacier bed when possible as well as from the interface between cold and temperate ice identified from the signal scattering to quantify the thickness of cold and temperate ice. Clarification added in the text (lines 75-76).
82: What's the spatial resolution of this imagery and of the DEM produced? This is now specified in the revised manuscript (lines 85).

86: Perhaps the model variogram parameters could be shared in an appendix?
The model variogram sill and range has been added in the text (lines 90). 95: (gr) "greater than what could be supported by ice deformation alone and thus implies basal sliding" The sentence has been corrected 103: I'm not sure that I understand "are calibrated ... by linear regression method" well enough that I could reproduce this independently given the data. Could you describe in greater detail (or if this is done elsewhere, provide a reference)? We now provide more detail here (lines 109-111). eq1: do we constrain R as less than to equal to M, or is it possible for refreezing to exceed melting (either locally or over the entire domain)? R is locally less than to equal to M (it is defined by Eq. 6).
125: I don't think the description for rs/m is correct. As written, one would expect that it could never exceed 1, however in eq4 it's clear that it can (i.e. in the accumulation zone). Perhaps it would be better described as "the ratio of annual accumulation to melt"? Ok, we modified the sentence. eq3: Doesn't the second case imply that if the annual accumulation is half the melt, we model the radiative melt factor as if the surface we ice-covered half the time? It's not obvious to me that this should be true. Yes, this the principle of this parametrized approach. Since the parameters and are calibrated to fit the mass balance data, their values may compensate the uncertainty linked to the approximation of / made here. This approach has been used in Gilbert et al. (2016) where it provided a robust estimate of the surface mass balance despite of the crude approximation made to calculate / . eq5: Since this is an annual precipitation rate summed over 365 days, don't we need to divide by 365 somewhere? Yes, this has been modified in the revised manuscript. I also find the units given for dP/dz surprising -I would have expected it to be m w.e. / (a m) to match Pref. The reference precipitation Pref is multiplied by a factor depending of elevation which is expressed by (1+(z -zref)*dP/dz). To keep homogenous units (z -zref)*dP/dz should be without unit so dP/dz is in m -1 . In the previous version of the manuscript we made a mistake in this factor, this has been corrected (line 137).
144: Is this parameterized in any way? Coefficients or parameterizations used in the flow law should be described. This is not parametrized, the Stokes equation are fully solved. It has been clarified. 153: (gr) "well-constrained," Done eq8: It'd be nice to drop the parentheses around T to disambiguate it from function application (as it's used on the LHS and in one of the integration limits). Also, I'm not sure we've defined $T_m(P)$ yet. Done. Tm is defined in the sentence following the equation.

171: can this assumption (crevasses go to the bedrock) be justified?
Our assumption is not that crevasses go to the bedrock but that liquid water is able to percolate through ice in the crevassed areas via cracks opening under water pressure. Indeed, it has been shown that meltwater is able to penetrate up to 400m of cold ice in Greenland (Lüthi et al., 2015). This is also theoretically confirmed by Van Der Veen (2007) where the author shows that crevasses subjected to inflow of water will continue to propagate downward until the bedrock with the propagation speed controlled primarily by the rate of water injection. Those references have been added in the revised manuscript (line 178).
197: Does the grid move with the surface or is it fixed? The grid moves with the surface and the variables are interpolated from the old grid to the new grid at each time step. This is now specified in the revised manuscript.
209: Can we simplify eq12 by expressing df in aˆ-1 (to avoid the scaling parameter)? Ok done.

215: In equation 14, should F be Fref?
No, F is the firn thickness at a given time whereas Fref is the initial firn thickness.
216: The notation might be a bit muddled here -the function is parameterized by zf (undefined, I think?), but that doesn't appear on the RHS. A diagram might help.
We re-organize the description of the computation of the density profile which is now clearer. We also checked that all the variables and parameters are well defined.
225: I'm pretty sure the units don't work out here -c(t) is accumulation per day, but df.F is per second The equations are now homogeneous. 226: Again, I think a more structured way of describing the algorithm, like a simplified code listing, would be helpful; the text feels too ambiguous. We re-organize the text in a way that the procedure we have done here is not ambiguous anymore.
234: "shifted our temperature forcing" -how and how much? (If this is described later,I had a hard time making the connection.) This is described in the result section 4.2: "Balanced conditions for the 2014 geometry are reached for a climate that is 0.7°C colder than the 1980-2016 climate with an Equilibrium Line Altitude (ELA) of 5770 m a.s.l. (19805770 m a.s.l. ( -2016Fujita and Nuimura, 2011). This provides a mean surface mass balance and a melting rate to force the steady state glacier simulation". We now provide more details in this section 3.4.2. (lines 250-251).
236: I'm not quite sure at first what "reported to the bedrock topography" means; are you altering the bedrock by an amount equal to the free surface change? Later in the paper I gather that you also changed the bedrock significantly in places where there is ground truth, which I think deserves justification. Yes, this is what we have done, it is now clarified. Our bedrock correction avoids major flow divergence in the velocity field that would result in strong vertical ice advection that lead to unrealistic steady state temperature field. Also radar measurements have their own uncertainties and cannot be considered as ground truth especially in temperate area where the bedrock reflection can be pretty weak and undetermined. We therefore favor a bedrock topography that satisfy mass conservation with the prescribed surface mass balance. Also Figure 6c shows that bedrock correction where radar measurement are existent rarely exceed 20 m apart of a really few exception (5 points of the 2010 measurements). We added a comment in the revised manuscript about this point (lines 306-307).
246: "reasonable accordance with the observation" -what exactly does this mean, and what are the criteria for "reasonable"? And do we even expect the observations to be similar to the steady state? We mean here that, for the purpose of the study, which is a thermal regime study, our method allows a good approximation of the glacier dynamic and geometry. We delete this sentence here since the result section provides the quantitative comparison needed to evaluate how good the accordance is. This comment in our manuscript was unnecessary here. 261: How do we know that scattering isn't due to the crevassing itself? See explanation in the General Comments. We added a short discussion in the revised manuscript (lines 275-279) 287: s/inexistent/nonexistent Done 325: (gr) not a complete sentence The sentence has been corrected.
339: I presume this means "no snow or firn over the surveyed part of the glacier," as there does seem to be a large part of the glacier above the ELA where we might expect firn? The mean ELA over the period 1980-2016 is at about 5880 m (Figure 4b) (contrary to the steady state ELA at 5760 m (Figure 4a)). Over the last ten years the ELA was even above 5900 m which is in accordance with the fact that no firn were present over the surveyed part of the glacier in 2015.
340: What are the uncertainties in the surface DEM? Are they meaningful (i.e., they're used as an input to create the modelled bedrock topography, IIUC)? Surface DEM was made from high resolution Pleiades images resulting in uncertainty bellow 1m. The surface DEM was acquired on November 7, 2014, 6 months before the radar measurement. With a mean thinning rate of 0.5 m w.eq. yr -1 , the uncertainty introduced by the temporal lag between the two measurements (radar and surface DEM) should be inferior to 50 cm making the uncertainty on the surface topography negligible in comparison to the one coming from ice thickness estimation. The bedrock topography uncertainties are therefore not significantly affected by the surface DEM uncertainties. We added a sentence in the revised manuscript (line 370).
346: It's interesting that these uncertainties are much smaller than the differences between the measured is thickness and the modelled ice thickness, and it's hard to believe that 20 m of horizontal uncertainty account for the rest. Any idea where the remaining difference could come from? Is it all from the assumed friction parameter? These uncertainties are coming from the theoretical vertical resolution of the radar measurements. However, the ice/rock interface is manually picked on the radargram which can introduce extra uncertainty when the reflector is weak (especially in temperate area, see figure 2). And yes, as you pointed out, the modeled thickness is strongly linked to the friction parameter which is not known but also by the surface mass balance which is modeled and also introduces uncertainties in the reconstruction. We add a sentence in the revised manuscript (line 367).
356: s/delicate/difficult Done 359: I might be misinterpreting this statement, but I'm not sure I agree; the advection of heat should depend on whether motion is at the surface or distributed throughout the thickness, shouldn't it? Yes, this sentence is not clear, the way that motion is distributed is important. We meant here that if the modeled surface velocities are in accordance with the measurements, it is likely that the advection processes are well represented because the way that motion is distributed is solved through the stokes equation that should accurately represent ice deformation in our 3D setup. We removed this sentence of the manuscript which led to confusion.

Referee #2 (Martin Lüthi)
The manuscript The influence of water percolation through crevasses on the thermal regime of Himalayan mountain glaciers by Gilbert and others is a very nice study on the poly-thermal regime of mountain glaciers. The paper describes a good modeling study that can explain the observed interesting polythermal structure.
Nevertheless, a lot of small changes are needed to render the manuscript ready for publication. Obviously, we as non-native English "users" have a big disadvantage, and our formulations might be not idiomatic. Despite this, care should be taken to at least consistently use "the", (the apparently random) capitalization rules, proper words ("localization" instead of "position" or "area", English -if somewhat similar -is not some kind of French), etc.
Once the many below concerns have been addressed (and many more which I did not all mention, since this is more time consuming than just rewriting the text), the paper will be ready for publication.

General comments
Overall, this is a very nice paper on an important topic. The modeling study is nicely crafted, comes to meaningful results, and elucidates some very relevant processes in poly-thermal glaciers. What I really liked is • good and comprehensive Introduction, • very good and thorough approach to the problem at hand, • very nice modeling study to explain noteworthy features in a unique data set, • good, meaningful Discussion of the important processes and their general significance.
This being said, the paper unfortunately is presented with many small warts that should be cured before it is ready for publication.
It took me unusually long to write this review, especially since I commented on many small things that should have been improved by the proof-readers before submitting the manuscript. Using colons (:) before equations, citing equations without parentheses around equation numbers, using × for multiplication in formulas etc. are things I have never seen in any of the usual journals. Please adapt to the conventions of the journal. Also, the paper would certainly profit from streamlining by a native speaker (I only indicated a few (many!) small issues).
We are grateful for the careful job you have done in pointing out all those small issues. We have taken into account all of them and the manuscript has also been streamlined by a professional English proofreader.
In the Methods section care should be taken to clearly describe the algorithms used, especially these rather ad-hoc rules that describe meltwater percolation. Especially section 3.3 is very opaque, and would profit from a flow diagram or a formula describing the iterations (I think). Also the algorithm in Section 3.4.2 is not easily understood, and details on the procedure are missing in some steps (e.g. Step 2). Some formulas (e.g. Eq. 14) appear to contain errors (outlined below).
With the help of reviewer 1 as well, we improved the clarity of those sections by adding more details and restructured how the different steps are ordered.
Some figures, and many figure captions should be adapted and improved. We modified and improved all the figure according to your comments.
My feeling is, that the main results could be explained with less figures, and some of them could be moved to an appendix / supporting material section.
We decided to keep all the figures as 12 figures still fits the format requirements is The Cryosphere and we prefer to have all the information in a single document.
The bibliography appears complete with DOIs (although given as URL). ISBN/ISSN should be given for the cited books.

Specific comments
41 "the small number of boreholes gives ..." Done 43 "extrapolated". The whole sentence is somewhat awkward, better reformulate and split in two. Done 45 "Scattering of electromagnetic waves ..." Done 46 start new sentence after "GPR data". Done 47 also Ryser et al. (2013) nicely showed the relation between ice temperature and scattering. We added this reference.
50 "rare" instead of "rearer" Done 74 or168 mμs−1 Done 84 better "imagery" Done 86 "the Kriging algoritm" (give reference). Done. We added a reference. 92 and 94: "stake" (singular) Done 95 "support" seems wrong here, better use "can be explained". And then, one wonders which ice flow parameters and stress regime have been used to arrive at this conclusion. We made the change in the revised manuscript. This statement comes from the cited reference (Fujita et al., 2001) where the authors used measured ice thickness and slope to estimate surface velocities using realistic flow rate factor in the shallow ice approximation (slab). Comparison with observed surface velocities that are much higher led to the conclusion that sliding occurs under Rikha Samba Glacier. We added details in the revised version of the manuscript (lines 99-101).
96 Since ice temperature is the main focus of this study, one is curious about the measurement process. Were the holes drilled mechanically or thermally? Was temperature measured once, or logged? Type of sensors and measurement equipment would also merit description. 100 "as input data" (leave away "an") Done 104 leave away "method" Done 106 "aims at identifying" Done 106 "observed": better say "deduced from" Done 107 no comma after "which" Done 112 leave away the colon (here, and before all equations). TC does not use this style. Done. We removed it before all equations 115 only give units once (they have to be the same anyway). Ok done 117 this contradicts line 111. Only say you used and improved model there.

Ok done
118 How does short-wave radiation affect ice dynamics (leave away)? We removed this sentence. It was not necessary here. 125 from where do you know the frad values? Those parameters are used as tuning parameters and are calibrated to match the measurements. This is detailed in the result section and in Table 2.
142 what value do all these factors have? Are they taken from literature, or determined from local measurements? This is also detailed in Table 2. We added a sentence to refer to the result section about this. 151 Using Q for a source term is unfortunate notation, since often Q denotes fluxes. Indeed, it is called a flux on line 174. Using Q for the volumetric flux coming from latent heat release (often called "heating rate") is a common notation used in most literature. We kept the notation but now refers to heating rate instead of source term to be consistent with the literature refereeing to this term.
152 σ and ε should also be written bold in the text Done 153 "constrained" Done 155 replace "defined" by "written in terms of" Done 159 A dot is missing in the number 3·10 3 The number seems to be correct at this line.
162 Here it would be interesting to say what the vertical resolution in meters is. We now specify the vertical resolution here (~10 m).
162 What kinds of elements have been used? Geometry and approximation order should be mentioned here. The elements we used are now described. "The mesh is constructed from 2D triangular linear element extruded in the vertical direction. It gives a mesh made of triangular prism unstructured in the horizontal direction and structured in the vertical direction." (line 169-170).
171 Why and how is water percolating to the bed? Ice (even temperate ice) is pretty impermeable.

"neglected"? Not clear what you want to say with this. Neglected from what?
We meant that surface meltwater coming from outside of the crevassed area through surface runoff is neglected in the amount of water available for refreezing in crevassed area. The revised manuscript has been clarified (line 180).
174 Why a heat flux, and in which direction? This should be a source term in Equation (7). We now call it heating rate, yes it is a source term.
175 How is the amount of refreezing water determined? This is described in the following paragraph. The annual surface melting is distributed vertically as available latent heat from refreezing as long as the ice temperature is cold.

Top to bottom of what? Of the glacier or of the layer?
Of the glacier. We added it in the revised manuscript. Yes, we now call it the heating rate.
185 This whole description of the algorithm is somewhat opaque. Please consider a flow diagram, or a clear description of what happens. This was also pointed out by reviewer 1. We now clearly specify the different steps to make the description clearer. We also moved this part in the section 3.4 about steady state strategy. (line 239-244).
186 "steady state": what are you doing here? Do you mean you do a fixed-point iteration until convergence? Are the steps iterative steps (for the solution at one grid point), or time steps? I doubt that you calculate a steady state of the whole glacier here. We calculate a steady state for the whole glacier here. We clarified the procedure in the revised manuscript (line 239-241).
191 "distribute": I think you use the 1-d model at every grid point, independent of all other grid points? If so, please say that. Yes, we now specify this (line 200-201).

"for the 3D model" Done
194 "It provides a high ..." Done 200 how sensitive is the temperature profile to the choice of the Gaussian standard deviation? The temperature profile is pretty sensitive to this choice since higher standard deviation lead to greater positive degree-day and therefore more melting. However, this standard deviation is well constrained by the fact we impose than the mean temperature seasonal cycle conserves the amount of positive degree day and therefore the amount of modeled melting. We added a sentence about this (line 210-211).
203 "reaches" Done 209 why not just using a commensurate valued with time unit in years−1, then tyr could be left away. Yes, we changed the units here. 217 "Combining Equations (13) and (14) gives" with parentheses around equations, "Equations" capitalized, and the sentence without colons (like any other article in TC). This has been corrected.
222 The time step should be named ∆t (not the infinitesimal dt), and already named in the text. Done 227 why is "Topography" capitalized. Even if English has no meaningful rules, titles should be capitalized consistently. We tried to be now consistent with capitalized words in the revised manuscript.
229 "known" instead of "resolved" Done 233 "the" instead of "our" Done 234 leave away parentheses around β. Done 235 Was this "best" determined with an optimization? This best was manually estimated to minimize the difference between measurements and model. Indeed, we assumed a uniform coefficient in this first step. This is now specified.
235 What is the meaning of these units? This is the value of the uniform coefficient mentioned in this first step. We moved the value sooner in the sentence to make this clear.
236 What is done during the "reporting"? This is probably not the right expression, and something happens to update the bed topography. This now clarified as suggested by the reviewer 1 as well (line 254).

"permits" instead of "allows"
This sentence has been removed (see reviewer 1 comments).
249 "were performed" (also the past/present tense should be used consistently) Done 256 I'm not sure whether a comparison to other glaciers would rather belong to the Discussion. We moved this comparison to the discussion (line XX). Done 271 "provides" sounds wrong, why not "the model is in good agreement with ..." Yes, we corrected this.
287 "in areas without radar measurements" Done 303 weird sentence, please correct Done 308 "in equilibrium" Done 310 "the climate change" is pretty meaningless. Probably you mean "surface warming" or similar? Also line 317. Yes, we modified this. 312 please call this consistently "ERA-interim". The time series should also be shown ,maybe in Figure 4c. Ok. We now added the annual temperature from ERA-interim in Figure 4. 320 "crevasses" Done 321 "linearly increasing temperature", the trend is not increasing. Done 325 Do you mean polythermal glaciers in general. Then the "the" should be discarded. Yes, this sentence has been modified according to reviewer 1 comment.

this sentence is incomplete.
This sentence has been corrected.

weird sentence
We improved this sentence 336 weird sentence We improved this sentence 337 "values". Better rephrase Ok we changed "mean values" by "average" 338 "calibrated on data" is not proper English (IMHO) We corrected the sentence.

A GPS (Garmin) should be accurate to 5-10 meters, so what is the problem here?
Even if the horizontal accuracy of the GPS is within 5-10 m, the accuracy of the positioning decreases when moving on the steep slopes in the mountainous terrain due to e.g. a weaker satellite signal, unevenly distributed satellite coverage or the GPR antenna not being completely aligned because of the rough terrain.
344 Considering how much effort it is doing these measurements, why don't you use some cheap real-time corrected GPS, such as the Emlid Reach? In this paper, we used rental equipment with the GPS provided with the GPR system. Other, more accurate methods will be definitely worth using when possible in future.
352 "the friction ..." Done 355 complicated sentence, rephrase Done 355 Why is this "mass flux conservation" special? It is part of the solution of the Stokes equations.

Yes, the mass conservation is always satisfied when solving the Stokes equation but can lead to strong surface elevation change if the flux divergence significantly differs from surface mass balance. We meant here that geometry is in equilibrium with surface mass balance. We clarified this in the revised manuscript (lines 378-379).
362 "coming from": better "derived from surface melt" Done 368 "the Greenland ..." Done 370 "the thermal regime" Done 372 not really "observed", but "inferred for" Done 384 "position" instead of "localization" Done 385 "warming" instead of "climate change" (which is a generic catch phrase without any particular meaning for this mountain area -there could also be local cooling) Done 392 "combined" Done 393 "reveals" Done Done 395 "facilitates/permits/enables" instead of "allows" Done 396 "affects" Done 397 "the thermal..." Done 398 "surface temperature increase" or "surface warming" (again, climate change is unspecific) Done Eq 5 Why are you using % and (1/100) here? Its easier to understand if you just use the ratios. Yes we changed this.
Eq 12 write equations without "×" (also in many others) We removed the "x" in all the equations.
Eq 14 the integration variable should be dz, not dzf. Maybe the upper integration limit should be zf. This should be written carefully! We corrected the equation.
Eq 16 Do not use × in any of these formulas, they become unreadable. We removed the "x" in all the equations.    (589) What does "steady" mean in this context of an oscillating forcing? Is this a limiting cycle (stationary periodic response at depth)? Yes we modified the caption here.   Tab 2 "Precipitation Lapse rate": consistent capitalization! Also the units are weird, why not justm−1(although I think this should bea−1, so this is completely wrong). Done. The notation dP/dz is not adequate since we use a multiplicative factor, this why the unit is km-1. We changed the notation to avoid confusion.
Tab 2 Are the radiative rates per square meter? So the units are wrong. frad is multiplied by potential solar radiation (W m -2 ) to obtain melting rate (m w.eq. d -1 ). The unit (m w.eq. d -1 W -1 m 2 ) should be correct here.
Abstract. In cold and arid climates, small glaciers with cold accumulation zones are often thought to be entirely cold based.
However, scattering in Ground Penetrating Radar (GPR) measurements on the Rikha Samba Glacier in the Nepal Himalaya revealHimalayas suggests a large amount of temperate ice that seems to be influenced by the presence of crevassed areas. We used a coupled thermo-mechanical model forced by a firn model accounting for firn heating to interpret the observed thermal 20 regime. WeUsing a simple energy conservative approach, we show that the addition of water percolation and refreezing in crevassed areas using a simple energy conservative approach is able to explain theexplains these observations. Model experiment showsexperiments show that both steady and transient thermal regimes are significantly affected by latent heat release in crevassed areas. It makes half of the glacier base temperate, resulting in the dynamicsice dynamic mainly controlled by basal friction instead of ice deformation. TimescaleThe timescale of thermal regime change, in response to atmospheric 25 warming climate, is also greatly diminished, with a potential switch from cold to temperate basal ice in 50-60 years in the upper part of the glacier whileinstead of the 100-150 years that it would take 100-150 years without the effect of the crevasses effect. This study highlights the crucial role of water percolation through the crevasses on the thermal regime of glaciers and validates a simple method to take it into account for it in glacier thermo-mechanical models.

Introduction
The thermal regime of a mountain glacier controls its hydrology, flow rheology, and basal conditions affecting glacier 35 dynamics and, consequently,, which in turn affects its behavior in response to climate change. It influences erosion rates (Bennett and Glasser, 2009), potential glacier hazards (Faillettaz et al., 2011;Gilbert et al., 2015)), and water resources in the glaciated catchments (Miller et al., 2012). It is thus essential to understand the processes causing and maintaining temperate basal conditions, as well as the mechanisms leading to changes in the thermal regime of glaciers.
Very little is known about the thermal regime of the Himalayan glaciers due to the harsh conditions and logistical difficulties 40 makingwhich make direct observations challenging in the remote, high altitudes areas. Borehole temperature measurements, such as carried out on Khumbu, Yala and Gyabrag Glaciers in the Himalayas (Miles et al., 2018;Mae, 1976;Watanabe et al., 1984;Liu et al., 2009), provide direct observations of the glaciers' thermal condition. However, a restrictedthe small number of boreholes gives only very limited information about the spatial distribution of the ice temperatures within the glacier and need in that case to be extrapolate from numerical modeling to estimate the . An estimation of the thermal structure of the 45 glacier in such case requires the use of numerical modeling in order to extrapolate the boreholes measurements (Wang et al., 2018, Zhang et al., 2013. Scattering of the electromagnetic waves signal in glacier ice is commonly interpretedregarded as diagnostic ofdiagnosing temperate ice in ground penetrating radar (GPR) data, and continuous. Continuous GPR profiles can thus provide information about the spatial distribution of thermal ice zones within a glacier (e.g. Ryser et al., 2013;Wilson et al., 2013;Gusmeroli et 50 al., 2012;Irvine-Fynn et al., 2006;Pettersson et al., 2007). Wilson et al. (2013) showed that the interface between cold and temperate ice matching with the localization of temperatures reaching pressure melting point in the boreholes could be identified with a 10 MHz GPR on two sub-Arctic polythermal glaciers. In the Himalayas, such GPR data are rarer whilerare although Sugiyama et al. (2013) showed with GPR data that the Yala Glacier in Nepal is polythermal, which was in agreement with two previous two borehole measurements in the ablation and accumulation areas of the glacier (Watanabe et al., 1984). 55 In this study, we revealdetermine the polythermal structure and ice thickness of a high altitude glacier in the Nepal Himalaya using GPR. We combine GPR data from 2010 and 2015 with other field data to determine ice thickness and to estimate the amount of temperate versus cold ice in the glacier. Measurements are interpreted using a 3D thermo-mechanical model for which we have developed new methods in order to;: (i) determine the thermal surface boundary condition, and (ii) take into account water percolation and refreezing in the crevassed areas. The model is forced by a surface mass balance model 60 calibrated with the field measurements, andwhich is then run to determine the steady state and transient thermal regimes of the glacier. We compare our modeling results with the GPR data to concludedraw conclusions about the processes defining the thermal regime of the glacier, and to provide recommendations on how to take them into account forin further modeling studies.

Study area 65
The Rikha Samba Glacier is a south-east orientated, medium sized glacier (5.5 km 2 ) located in the Hidden Valley in western Nepal (Fig. 1a). The glacier iswas about 5.5 km long with an elevation ranging from 5420 to 6440 m a.s.l. in 2014. The glacier melt water contributes to the Kali Gandaki Basin of the Ganges River. The Hidden vValley falls under a rain shadow climate, receiving the least precipitation in Nepal with an annual precipitation of 370 mm (Fujita and Nuimura, 2011). The annual mean temperature measured with an automatic weather station (AWS) in the vicinity of the glacier (Fig. 1a) was -5 °C in 70 2014. The glacier was first visited in 1974, and it hasd been losing mass at a mean rate of -0.5 m w.e. yr -1 between 1974 and 2010 (Fujita et al., 1997;Fujita and Nuimura, 2011).

Ground penetrating radar (GPR)
We used a Malå GeoScience ProEx impulse-type ground penetrating radar (GPR) with a 30 MHz Rough Terrain Antenna (RTA) with a frequency band from 15 to 45 MHz to measure the ice thickness and thermal regime of the Rikha Samba Glacier 75 in 2015 (Fig. 1b). The continuousRTA is an unshielded radar antenna in which the transmitter and receiver elements are mounted in a row in a flexible tube making it suitable for rough terrain such as glacier surface.
Continuous profiles were filtered and some of them migrated, and while the GPR reflectors from the assumed bedrock surface were manually picked from the data. Picked two-way travel times of the radar signal were converted to ice thickness using a wave velocity of 1.68 × 10 8 168 m sµs -1 (Robin, 1975). Strong scattering of the radar signal within the ice was interpreted as 80 temperate ice whereas ice without internal reflectors was classified as cold ice., and the interface between these zones was also manually picked in the data. In addition, we have also used point data collected in 2010 with an impulse GPR transmitter (Ohio State University) with a set of half-wavelength 5 MHz dipole antennas (Fig. 1b). These data were used only for ice thickness.
The time increment of five years between the GPR measurements was corrected by projecting the 2010 data to 2015. This was done by assuming that the glacier was thinning withat the same rate between 2010 and 2015 as the long-term thinning rate for 85 100-m elevation bands obtained for a 12-year period between 1998 and 2010 (Fujita and Nuimura, 2011).

Glacier geometry and crevasse localization
A digital elevation model (DEM) of 30m resolution was generated for the Rikha Samba Glacier from Pleiades tri-stereo satellite imagery (0.5m resolution) for November 7, 2014. Crevassed areas on the glacier were visually identified from the imagery and refined using Google Earth imagery (WorldView, September 21, 2012) in which the crevasses were more visible. 90 The ice volume and bedrock topography were initially estimated by defining ice thicknesses as zero at the margins of the glacier and interpolating the ice thickness data measured with the GPR. For interpolation, we assumed a spherical semivariogram (fitted to data with a sill of 3743 m 2 and a range of 784 m) and applied krigingthe Kriging algorithm. (Matheron, 4 1963). This method is widely used to interpolate ice thicknesses measured with a GPR to estimate the volumes of mountain glaciers (e.g. Fischer, 2009). Since the GPR data do not cover the entire glacier, it results in high uncertainties in the 95 interpolated bedrock topography in some part of the glacier. The initial bedrock topography is thus corrected using the ice flow model (Sect. 3.4.2).

Glacier mass balance, surface velocities, and ice temperature
We constrained the surface mass balance model from stakes network measurements in 2012 and 2013 (Gurung et al., 2016) and from the total volume change estimated by a geodetic survey (Fujita and Nuimura, 2011) over the period 1974-1994 0.57 m w.e. a -1 ) and 1998-2010 (-0.48 m w.e. a -1 ). Stakes displacement monitored during the 1998-1999 period shows horizontal surface velocities between 9 and 24 m a -1 , which is (Fujita et al., 2001). By comparing those velocities with the deformation velocities computed from shallow ice approximation and realistic flow rate factors, Fujita et al. (2001) concluded that observed velocities are greater than what could support the be explained by ice deformation of ice andalone, which thus revealing the existence ofimplies basal sliding (Fujita et al., 2001). Ice. The ice temperature at 10 m depth was measured with 105 a thermistor chain manufactured by Stum Bohr AG Switzerland in the lower ablation area (5600 m a.s.l., Fig. 1b) for 2014-2015. AirThe temperature was logged over one-year period in six-hour intervals. The air temperature and precipitation were observed with an AWS in the vicinity of the glacier at 5320 m a.s.l. (Fig. 1a).

Meteorological parameters
We used the ERA-Interim reanalysis data (Dee et al., 2011) at a daily timescale over the period 1980-2016 as an input data 110 for the mass balance model. We only used air temperature data and assumed a constant precipitation rate in time (no precipitation seasonality) to avoid complexity in the simulation. Temperature and precipitation awere then distributed on the glacier according to altitudinal gradients to reproduce the observed mass balances. Bias in the ERA-Interim air temperature are calibratedwas corrected using a linear function determined from the comparison with the local AWS data (Fig. 1b) over the period 2011-2015 by linear regression method.. 115

Modeling methods
The modeling study aims to identifyat identifying which physical processes lead to the thermal structure observed bydeduced from the GPR measurements. First, we focus on steady state simulation for which, ice flow and thermal regime are in equilibrium with constant surface boundary conditions (surface temperature and mass balance). We then use the steady state simulation as initial condition offor the transient model experiments. 120

Surface mass balance model
Mass balance is modelled using a degree-day method following Gilbert et al. (2016). Net) in which we include the influence of the spatial variability of the shortwave radiation. The net annual surface mass balance is determined by: where is the net annual surface mass balance (m w.e. a -1 ), is the annual snow accumulation (m w. e. a -1 ) , is the rate of refreezing (m w. e. a -1 ) and is the annual melting (m w. e. a -1 ).. Meltwater is computed from the sum of two components following Pellicciotti et al. (2005): In this study, we updated a degree-day model described in Gilbert et al. (2016) by including the influence of the spatial variability of the shortwave radiation to constrain both ice dynamics and thermal regime of the glacier. Meltwater is computed 130 from the sum of two components (Pellicciotti et al., 2005): where is the daily melt (m w.e. d -1 ), is the air temperature (K), ℎ is a temperature threshold for melting (K), is the melting factor (m w.e. K -1 d -1 ), is the potential solar radiation (W m -2 ) and is the radiative melting factor (m w.e. W -135 1 m 2 d -1 ). Following a similar approach as used in Gilbert et al. (2016) 140 The annual ratio / is computed assuming that: where, is the dailyannual precipitation rate (m w. e. a -1 ) at the elevation (m a.s.l.), ⁄ is the altitudinal precipitation lapse rate (% factor (m -1 ), is the elevation (m a.s.l.) and is a temperature threshold that distinguishes between snow and rain (K). Assuming that refreezing in the previous year creates impermeable layers and thus occurs only to 150 a depth equal to the annual accumulation rate, we write: where is a refreezing factor. The determination of the parameter values is described in section 4.2.

Thermo-mechanical model 155
The ice flow model is based onsolves the Stokes equations for incompressible flow adopting Glen's flow law for viscous isotropic ice (Cuffey and Paterson, 2010) and is coupled to an energy conservation equation using the enthalpy formulation where is the firn/ice density (kg m -3 ), is the time (s), is the enthalpy (J kg -1 ), is the ice velocity vector, is the enthalpy diffusivity (kg m -1 s -1 ), is the source termheating rate coming from meltwater refreezing (W m -3 ), and (̇) is the strain heating (W m -3 ) with and respectively the stress and strain rate tensors. A basal heat flux of 4.0 × ·10 -2 W m -2 is assumed for the basal boundary condition. This value is not well -constrained ranging from 2.0 × ·10 -2 W m -2 (observed at Rongbuk 165 Glacier in the Everest region (Zhang et al., 2013)) to 8.0 × ·10 -2 W m -2 as predicted by a large scale model . The enthalpy is defined fromwritten in terms of ice temperature (K) and water content ω: with where is the heat capacity of ice (2.05 ·10 3 J K -1 kg -1 ), 0 is the reference temperature for enthalpy (set to 200 K), is the 170 melting point temperature (K), is the latent heat of fusion (3.34 ·10 5 J kg -1 ) and is the ice pressure (Pa).
Changes in the glacier surface elevation are computed by solving a free surface equation (Gilbert et al., 2014a). We adopt a linear friction law as a basal boundary condition for the Stokes equation with where is the basal shear stress (MPa), is the sliding velocity (m a -1 ) and is the friction coefficient (MPa a m -1 ).
The model is solved using the finite-element software Elmer/Ice ( where is the basal shear stress (MPa), is the sliding velocity (m a -1 ) and is the friction coefficient (MPa a m -1 ).

Modeling crevasse influence via water percolation
In order to determine the areas where the crevasses are likely to form on the glacier, we compute the maximal principal Cauchy 185 stress (MPa) at the glacier surface from the stress tensor. We compare it with a threshold value ℎ (MPa) to identify where damage production occurs (Pralong and Funk, 2005;Krug et al., 2014) and define the crevassed areas where > ℎ .
In Crevasses are commonly considered to be the crevassed areas, weplace where surface meltwater can enter the body of the glacier (Fountain and Walder, 1998). The development of a hydrologically driven crack network initiated by water-filled crevasses seems to be able to transport liquid water down to the bedrock even through cold ice (Boon and Sharp, 2003;Van 190 Der Veen, 2007). In our model, we therefore make an assumption of free vertical percolation of the meltwater down to the 8 bedrock, in which local surface meltwater in the crevassed area is the only source of liquid water percolating into the crevasses.
This means that we leave out any water coming from the surface runoff and draining to outside of the crevassed area is neglected.and draining into the crevassed area. Assuming that water refreezes in the first cold layer, we compute a latent heat volumetric fluxheating rate (see eq. 7) from the available annual meltwater and the ice temperature of the current iteration. 195 At each vertical layer i of the model, is computed from the amount of refreezing water each year (kg m -2 ).) with where is the thickness of the layer (m) and is the timestep (s). The amount of refreezing water is distributed from top to bottom of the glacier with the condition that the enthalpy for a certain layer has first reached the value corresponding to 200 theenthalpy of fusion of waterHf before the water can access the next layer downwards. Starting from the surface melt 1 = , the amount of liquid water available for refreezing in the next layer +1 is computed following: Using the estimated flux, a new steady state enthalpy field is computed and can be updated from the new temperature 205 field. The procedure is repeated until reaching a steady state. WithIn adopting this approach, we assume that all the energy used to melt ice at the surface in crevassed area is released into the deeper ice body. It can be seen as an energy conservation approach rather than the modeling of water routing through crevasses.

Enthalpy surface boundary condition including firn/snow influence 210
For this study, we develop a new method in order to determine surface boundary conditions for enthalpy. We use the 1D semiparameterized approach developed in Gilbert et al. (2014b) and distribute it over the entire glacier. The method takes into account water percolation and refreezing in both firn and seasonal snow to determine the adequate surface boundary condition offor the 3D model (Gilbert et al., 2012). The 1D model is solved on a one-dimensional 10-m-depth vertical profile at each surface node of the 3D model. independently of all other surface nodes. It allows a higher temporal and vertical spatial 215 resolution in order to explicitly solve percolation and refreezing processes. in the snow and firn.
Starting from an initial uniform temperature profile, firn/ice temperature is solved at daily time steps along the vertical profile with a 0.06 m resolution. The vertical grid is updated as the surface elevation evolves and all the variables are interpolated from the old grid to the new one at each time step. The 1D model is forced by air temperature and by the surface mass balance model that provides snow accumulation and surface melting for the corresponding surface node. To compute the steady state 220 condition, the model is driven by a mean annual cycle of air temperature which is determined at the daily time scale from the meteorological data. A Gaussian random noise is added to the computed mean annual temperature cycle to plausibly represent the daily temperature variability. The standard deviation of the Gaussian function is adjusted to match the number of positive degree-days in our mean annual cycle to the mean one in the data. The 1D model run during in order to conserve the amount of modeled melting. It provides a strong constraint on this parameter making our results independent of this choice. The 1D 225 model runs for several years with the same cycle until the 10m-depth temperature (approximate limit of the thermally active layer) reaches a mean annual equilibrium value that will be used as the surface boundary condition of the thermo-mechanical model.
At each surface node, the initial density profile is calculated from the steady state firn thickness FrefF (m w.e.) which is computed from the steady state mass balance assuming: as a function of time as follow: 230 where ( ) is the daily net surface accumulation (m w.e. yr -1 ), is a firn densification rate parameter (syr -1 ) and tyr (s) is one year in second. Δ the timestep (yr). The initial value of is set to which is computed from the steady state mass balance assuming 235 The density is then calculated assuming a linear evolution of density with depth between the surface density 0 (kg m -3 ) and the ice density (kg m -3 ) at the firn/ice interface. It gives: where z is the vertical coordinate (m a.s.l.), is the elevation of the surface (m a.s.l.), and is the elevation of the firn/ice transition (m a.s.l.). From mass conservation, zice has to satisfy: where is the water density (1000 kg m -3 ), and is the firn thickness (m).). Combining the equation 13 and Equations (14) and (15) gives: In order to take into account the snow seasonal variability due to snow/rain threshold (Table 1), the density profile is updated 250 at each time step by computing the evolution of : where ( ) is the daily net surface accumulation (m w.e. d -1 ). The density profile and the surface elevation are updated only if > by adding the corresponding amount of snow at density 0 . The initial value of is set to . 255 The surface enthalpy field obtained using the distributed 1D model is used as a surface Dirichlet condition to solve the enthalpy equation in the 3D model. The 3D thermo-mechanical steady state of the glacier is then computed by following an iterative procedure: 1. Compute the steady state enthalpy field from Eq. 7 with Qlat = ⃗ = () = 0 forced by the surface enthalpy field given by the distributed 1D model 260 2. Compute Qlat from Eq. 10 and Eq. 11 3. Solve the ice velocity and deformational heat from the Stokes equation.
4. Compute the steady state enthalpy field from Eq. 7 with Qlat determined in step 2 Steps 2, 3 and 4 are repeated until a steady state is reached.

Bedrock Ttopography and basal sliding condition 265
The main challenge in determining the glacier thermo-mechanical equilibrium is that: (i) the bedrock topography is not resolvedknown everywhere underneath the glacier, and (ii) the glacier is sliding, which means that a friction coefficient has to be quantified. In order to resolve these issues, we usedadopted the following approach.steps: Step 1: Starting from the measured surface topography in 2014, we run the coupled thermo-mechanical model with the interpolated bedrock topography and a uniform basal friction coefficient during a 10-year period forced by the steady state 270 surface mass balance and enthalpy. In order to obtain athe steady state mass balance, we shifted ourthe temperature forcing (1980-2016 corrected reanalysis) of 0.7°C to obtain balanced mean conditions during the simulation period. The resulting 11 mean surface mass balance is used as a steady mass balance associated with the 2014 glacier geometry. Here, we assume a uniform friction coefficient ( (10 -2 MPa a m -1 ) that allows the best match with the measured surface horizontal velocities (10 -2 MPa a m -1 ). ; 275 Step 2: The computed changes in the free surface in Step 1 are reported to the bedrock topography. by altering the bedrock by an amount equal to the free surface change (Gilbert et al., 2018); Step 3: After a few iterations between Steps 1 and 2, we obtain a corrected bedrock topography where major flux divergences are avoided. Using this bedrock topography and the measured surface topography we inverted for the friction coefficient by constraining the surface velocity onwith emergence velocities, which are taken as opposite to surface mass balance. This is 280 done by using a controlled inverse method to minimize a cost function defined from the misfit with measured surface velocities and a regularization term (Gillet-Chaulet et al., 2012;Gagliardini et al., 2013). Following Gilbert et al. (2016Gilbert et al. ( , 2018 we define the cost function from the misfit between modelled and measuredestimated emergence velocities. ; Step 4: We finally run the model using the corrected bedrock topography and the inverted friction coefficient until the surface topography reaches a steady state. 285 This method allows reaching a thermo-mechanical equilibrium in which surface topography and velocities are in reasonable accordance with the observation which allows a realistic study of the glacier thermal regime.

Transient evolution
Transient simulations awere performed at yearly timestep using the steady state glacier as initial condition. Surface mass 290 balance and enthalpy are updated each year from the surface model described previously and forced by daily temperature reanalysis. We assumed constant basal friction parameter through time.

Thermal regime and ice thickness measured with the GPR
GPR data show that the Rikha Samba Glacier is a polythermal glacier consisting mainly of cold ice but with significant 295 temperate areas (Figs. 2 and 3). The measured maximum thickness was 178 ± 2 m in the middle part of the glacier where the surface slope is relatively gentle. In contrast to the ice of Yala Glacier, another polythermal glacier in the Nepal Himalaya (Sugiyama et al., 2013), temperate ice is also found in the ablation area of Rikha Samba Glacier and only the lowermost part of the glacier where ice thickness is less than 25 m is completely cold in this area. Ice temperature measurements by the thermistor chain support the GPR interpretation of an upper cold ice layer with sub-zero temperatures at the depth of 10 m 300 (annual mean -2 °C) in the ablation area of the glacier (Fig. 1). A notable characteristic of the GPR -based thermal regime is that the occurrence of temperate ice localization seems to be associated with the presence of surface crevasses (Figs. 2 and 3).
The way that scattering is advected by the ice flow shows that the crevassed areas behave as a source term of scattering. It 12 seems unlikely that the observed scattering would be produced by the crevasses field itself since the scatter occurs in the deeper part of the glacier where crevasses are not expected to persist. It would be rather more compatible with the formation of 305 temperate ice which is then advected by ice flow.

Surface mass balance
We run the mass balance model using the 2014 surface topography over the period 1980-2016 using calibrated temperature reanalysis data (ERA-Interim) and (Dee et al., 2011) while assuming a constant precipitation rate. The parameters were constrained by the stake measurements in 2012/2013 (Gurung et al., 2016), meteorological observations (Gurung et al., 2016)), 310 and geodetic mass balance over the periods 1974-1994(Fujita and Nuimura, 2011 (Fig. 4). The parameters are summarized in Table 2.
Balanced conditions for the 2014 geometry are reached for a climate that is 0.7°C colder than the 1980-2016 climate with an Equilibrium Line Altitudeequilibrium line altitude (ELA) of 5770 m a.s.l. (the 1980-2016 ELA is 5880 m a.s.l.; Fujita and Nuimura, 2011). This provides a mean surface mass balance and a melting rate to force the steady state glacier simulation (Fig.  315 4a).
The model provides ais in good agreement with the observations but is not able to reproduce the same mass balance distribution as observed in 2013 (Figs. 4b-d). InterannualProbably, the inter-annual variability of the mass balance produced with our mass balance model is probably not very well represented since we assume a constant precipitation rate. Furthermore, the Rikha Samba Glacier is a summer -accumulation -type of glacier in which precipitation events in summer can significantly affect the 320 mass balance through albedo feedback (Fujita and Ageta, 2000;Fujita, 2008). However, the long-term trend and the mass balance gradient agree with the observations, which is satisfactory for the purposes of this study focusing on the thermal regime.

Enthalpy surface boundary condition
The modeled upper boundary condition field (Fig. 5) revealedreveals a mainly cold surface condition with a temperate band 325 between 5800 and 5900 m a.s.l. where both melting and firn thickness are sufficient to maintain temperate conditions. In the higher part of the accumulation area, the water percolation occurs only in the first two meters due to the limited amount of meltwater resulting in cold temperature at 10 m-depth (see Site 1 in Fig. 5),) whereas, lower down at Site 2, meltwater percolates deep enough to keep temperate conditions all year round at 10 m-depth. In the ablation zone (Site 3), water percolation is limited to the seasonal snow thickness resulting in a cold boundary condition. 330

Modeled steady state glacier
The modeled steady state glacier is in good accordance with measured ice thickness (Fig. 6), measured horizontal velocities (Fig. 7a), and observed crevassed areas (Fig. 7c). The correction made on the bedrock topography following the method described in section 3.4.2 greatly improveds the quality of the modeling in the parts whereareas without radar measurements 13 are inexistent (Fig. 6). A simple interpolation (Fig. 6a) leadswould have led to non-physical ice thickness with unrealistic flux 335 divergence, which areis avoided by our method. The bedrock correction in place where there is radar measurement rarely exceeds 20 m (Fig. 6c) apart from a few exceptions (i.e., 5 points in the 2010 measurements). We preferred to allow the corrected bedrock to differ from radar measurement in order to favour the mass conservation associated with the prescribed surface mass balance. Moreover, radar measurements are uncertain especially in temperate areas where the bedrock interface is weakly determined. The inversion of the basal friction coefficient (Fig. 7b) provides a final steady state where ice flow is in 340 accordance with steady state emergence velocities. The good agreement with horizontal velocities measured at stakes (Fig 7a) shows that our estimated emergence velocities (from surface mass balance) are consistent with the observed ice flow.

Thermal regime: influence of melt water percolation
Modeling results show that water percolation in crevasses strongly affects the steady state thermal structure of the Rikha Samba Glacier leading to large temperate zones even at the glacier bed (Figs. 8b and 9b). It significantly extends the temperate -based 345 parts, which cover almost the entire ablation area. Although we adopted a simple approach for water percolation through crevasses, modeled temperate ice thickness is in fairly good agreement with the GPR data (Fig. 8b). If water percolation through crevasses is neglected, the thermal regime of the glacier forced by mainly cold upper boundary conditions (Fig. 5) would result in a mainly cold based glacier (Fig. 9a). In thissuch case, cold ice advection from the higher part of the glacier is able to compensate for the temperate surface conditions of the lower accumulation zone (Fig. 5), and only two bands of 350 temperate ice are able to reach the bed on both sides of the flow line of the glacier (Fig. 9a). Such a thermal regime is in large disagreement with the observed amount of temperate ice from the GPR data (Fig. 8a). This indicates athe significant role of deep water percolation through cracks in cold ice as suggested by the GPR observations. We show that the use of The observed ( Fig. 8b and 9b) orand modeled ( Fig. 8c and 9c) crevassed area leads to a similar result andwhich validate our approach to model in modelling the location of the localization of crevasses. 355

Transient evolution
Despite the good agreement between the GPR data and the steady state model, a significant difference exists at the highest crevassed field. A temperate area is clearly visible in the GPR data ( Fig. 3a) whereas the steady state thermal regime model predicts cold ice (Fig. 8). Mass balance measurements show that the Rikha Samba Glacier has not been atin an equilibrium state for at least 40 years with an almost constantly negative rate of -0.5 m w.e. a -1 (Fig. 4c4d). This temperate area could be 360 therefore be the signature of a transient response to the climate changeair temperature warming. In order to investigate the potential impact of 40 years of unbalanced state on the glacier thermal regime, we performed a transient simulation starting in 1975 from the modeled steady state (experiment with observed crevasses location imposed from observations) and forced by the ERA-interim reanalysis time-series. This (Fig. 4c). It shows that the upper boundary conditions changed significantly with a cooling of the former accumulation zone in response to firn disappearance and a warming of the highest elevation due to 365 melting increase (Fig. 10b). After a 40-years run, a temperate zone, that diddo not exist at the steady state, had developed in 14 the highest crevassed area  as observed onin the 2015 GPR measurements (Fig. 10d). This result strongly suggests that the presence of temperate ice in this zone is a result of a transient response to the climate changeair temperature warming and increasing surface melting at the higher elevations. ResultsThe results also agree better with the observations, including the thermistor data (Figs. 10a). 370 To assess the sensitivity of the thermal regime to future temperature change, we performed a future simulation of the glacier retreat and thermal change until 2100 with (Fig. 11) and without (Fig. 12) water percolation through crevasses for a linearlinearly increasing temperature trend increase of +1 °C between 2014 and 2100 (+1.7 °C in comparison with the steady state climatic condition). This shows a much faster development of a new temperate area when water percolation in the crevasses is taken into account. In this casesuch a scenario, the glacier becomes almost entirely temperate based by 2050 (Fig.  375 11) whereas it would remainhave remained almost entirely cold if water percolation through the crevasses would behave been neglected (Fig. 12). This highlights the crucial role of deep water percolation through cracks in the thermal regime of the polythermal glaciers. A phenomenon that should be taken into account together with firn heating when modeling the past and future responseschanges of glaciers thermal regimes and retreat of glaciers.

Uncertainty
The modeled thermal regime is sensitive to the basal heat flux and the firn thickness, which are poorly constrained. Sensitivity experiment on those two parameters forThe sensitivity experiment performed for the steady state simulation shows that the amount of modeled basal temperate ice can vary significantly butdepending on the choice of parameters associated with firn thickness and basal heat flux which are poorly constrained. However, the thickness of the modeled temperate ice always 385 remains much less than indicated by the GPR observations if the crevasse influence is not taken into account. This means that uncertainty on those parameters cannotalone do not explain the disagreement between data and model when the role of crevasses is neglected. The mass balance model we used is simplified since seasonality and time variation in precipitation are not taken into account. However, the purpose of the study is not an accurate simulation of Rikha Samba glacier past and future evolution, but a study about its thermal regime. Our study relies on long-term mean value of surface mass balance. This should 390 be adequately captured by our simple model which is calibrated on data.
The mass balance model we used is simplified since we did not take into account the seasonality and time variation in precipitation which can be sources of uncertainties in our study. However, the main focus of this study is the thermal regime of the glacier and, hence, we do not claim that our simulation is accurate regarding the Rikha Samba Glacier past and future evolution. Our study relies on a long-term average of the surface mass balance and our simple calibrated model should be able 395 to adequately captured it.
Since the density of ice is well constrained and there was no snow or firn onover the surveyed part of the glacier at the time of the field measurement in 2015, the main uncertainties of the GPR measurements arise from the GPS positioning of the GPR measurements, the radar wavelength, and scattering of the radar signal. For the point measurements and those parts of the GPR profiles along which the bedrock reflection was clearly identified, the accuracy of the horizontal coordinates is about ±20 m, 400 especially on the steepest surface slopes. In addition, the vertical resolution of the GPR signal is usually considered to be approximately one quarter of the radar signal wavelength, which is about 5.6 m and 33.6 m for the 30 MHz and 5 MHz antennas, respectively. In other words, the vertical resolution of the englacial scattering, interpreted as temperate ice and ice thickness along the continuous GPR profiles, is about 1.4 m, whereas the same for the ice thickness obtained from the point measurements is about 8.4 m. In addition,The manual picking of bed reflectors in the radargrams introduces extra uncertainties 405 up to 20 m, especially in temperate areas where the bed reflectors are weaker (see Fig. 2). Finally, the limited coverage of the radar profiles on the glacier introduces uncertainty in the bedrock topography inferred from the GPR data even after correction using the model. Uncertainties in the surface DEM are expected to be lower than 1 m and do not introduce extra uncertainties in the bedrock topography reconstructed from the measured thicknesses. Our interpretation of the thermal regime based on the englacial radar scattering of the GPR 30 MHz profiles is supported by previously found close agreement between the observed 410 scattering and borehole temperatures without significant difference in observed englacial scattering relative to the expected measurement error at 10, 35 and 50 MHz antenna frequencies (Wilson et al., 2013).
In the modeled bedrock topography, the difficulty arises from the fact that the friction coefficient is unknown and we had to assume a uniform value of basal friction coefficient to correct the bedrock topography from flux divergence anomalies. The friction coefficient we inferred in a second step to force the steady emergence velocity to match the balanced surface mass 415 balance is therefore affected by ice thickness errors. The resulting velocity field is consistent with the surface mass flux conservationbalance but contains uncertainty uncertainties in the respective contribution of the basal sliding and ice deformation. to surface velocities remain. It makes delicate to interpretthe interpretation of the modeled basal friction difficult ( Fig. 7b), which thus has to be seen more as a tuning parameter rather than a parameter revealing physical processes. However, these uncertainties have little influence on the modeled thermal regime since advective processes will be still correctly 420 represented as long as the surface velocities match with the observations (Fig. 7a).

Role of surface runoff through supra-glacial streams
Our study shows that the thermal regime of the Rikha Samba Glacier can be modeled by taking into account melting occurring inthe liquid water derived from local surface melt percolating into crevassed fields and neglecting. Neglect of water inputs coming from the surface runoff through supra-glacial streams seems thus, to be a valid hypothesis. For a polythermal glacier 425 with a cold surface layer, which is a common feature in the ablation areas where firn heating is nonexistent, surface runoff occurs in well-marked and persistent streams at the surface (Ryser et al., 2013). Thoese streams bring water into the crevasse fields through a few localized entry points only. Thus,Hence, as there is only a relatively small contact surface between the streams and cold ice, and only a limited amount of water will refreeze (Fujita et al., 1996).. This is why the influence of surface runoff on the thermal regime of the glacier is likely to remain limited. Similarly, Lüthi et al. (2015) concludehave concluded 430 that moulins have little influence on the thermal regime of the Greenland Ice Sheet since they can be represented as line sources thatwhich provide limited warming ofto the surrounding ice. InOn the contrary, surface melting occurring on the crevasses fieldin crevassed fields is well distributed and can release latent heat on a much larger areas,area, thus having a stronger impact on the thermal regime.

Enhanced influence of climate change on glaciers thermal regime and dynamics 435
The influence of deep latent heat release through melt water percolation in crevasses have been already observed ininferred for Greenland. Lüthi et al. (2015) and Seguinot et al. (2020) observed temperature anomalies in borehole measurements that can only be explained by latent heat released down to a 400 m-depth in the crevassecrevassed fields. Similar conclusions have been made by Hills et al. (2017) although they show that these effects remain localized and may not really influence the thermal regime of the Greenland Ice Sheet aton a large scale. In the case of mountain valley glaciers, crevassed fields can cover a 440 significant fraction of the total glacier area. This, which is combined with a generally much faster ice flow leading to efficient advective processes that transport the heat produced in the crevassed areas. The results of these combined effects are significant and, greatly influenceing the thermal regime at the glacier scale as shown in this study. As already pointed out for ice sheets (Phillips et al. 2010(Phillips et al. , 2013, the timescale of the glacier thermal regime response to climate change is also greatly diminished compared toin comparison with the case where only heat diffusion/advection of surface changes are taken into accountinvolved 445 in the warming process. However, we show here that deep water percolation is likely restricted to the crevassed areas and absent elsewhere in order to reproducemodel the observed thermal structure. This restricts the spatial extent of the process on the glacier as it is dependent on the bedrock topography and related crevasse localization.its related crevasse positions. For instance, the Yala Glacier, another polythermal glacier in the Nepal Himalaya (Sugiyama et al., 2013), is entirely cold in its ablation area likely due to thinner ice and a lack of crevassed areas. Nevertheless, it is likely that the meltwater percolation via 450 crevasses has a significant impact on the thermal regime as highlighted here with the Rikha Samba Glacier, and that it is a common for allof numerous polythermal high altitude glaciers. Future climate changewarming could also lead to a faster thermal regime response than previously thought (Gilbert et al. 2015), especially in the cold accumulation areas where the preindustrial melting rates were not sufficient to form temperate ice.

Conclusions 455
In this study, we use GPR measurements to show that the high elevation Himalayan Rikha Samba Glacier is polythermal. We interpret the field observations of the temperate ice thickness using a 3D thermo-mechanical model constrained by a surface model taking into account water percolation in firn and seasonal snow. We show that the firn/snow heating, heat deformation, and geothermal heat flux alone cannot explain the observed amount of temperate ice. The combininged evidence of the model and observations reveals that valley-type mountain glaciers in cold climates are greatly affected by water percolation into 460 crevassed fields releasingwhich release latent heat into the ice body. It affects the thermal regime at the scale of the whole glacier making temperate ice zones much larger than they would be without this effect. This allowsenables sliding in large areas of the bed and largely affecting the glacier dynamics and ice thickness.
We also show that the thermal regime of the Rikha Samba Glacier is affected by a transient response to the last 40 years climate changeyears' air temperature warming extending the temperate area to the highest part of the glacier. The thermal changes are 465 occurring at a much smaller timescale (~50 years) due to the crevasse effect compared to what it would be by advection/diffusion processes onlyalone (>100 years). Our study reveals the crucial role of deep water percolation through cracks in determining both the steady state and transient thermal structure of the polythermal glacier. We provide a simple approach easily applicable to any glaciers for a more accurate reconstruction of complex thermal structures as the one observed onin the case of the Rikha Samba Glacier. 470 Data availability. As required by ICIMOD publication policy, GPR data used in this study will be made publicly available through ICIMOD's Regional Database System at: http://rds.icimod.org. ICIMOD has an open data policy. The modelling code 475 is based on the open-source code Elmer/Ice available at http://www.elmerfem.org. It includes an example based on the case presented in this paper.