Modeling near-surface firn temperature in a cold accumulation zone ( Col du Dôme , French Alps ) : from a physical to a semi-parameterized approach

Analysis of the thermal regime of glaciers is crucial for glacier hazard assessment, especially in the context of a changing climate. In particular, the transient thermal regime of cold accumulation zones needs to be modeled. A modeling approach has therefore been developed to determine this thermal regime using only near-surface boundary conditions coming from meteorological observations. In the first step, a surface energy balance (SEB) model accounting for water percolation and radiation penetration in firn was applied to identify the main processes that control the subsurface temperatures in cold firn. Results agree well with subsurface temperatures measured at Col du Dôme (4250 m above sea level (a.s.l.)), France. In the second step, a simplified model using only daily mean air temperature and potential solar radiation was developed. This model properly simulates the spatial variability of surface melting and subsurface firn temperatures and was used to accurately reconstruct the deep borehole temperature profiles measured at Col du Dôme. Results show that percolation and refreezing are efficient processes for the transfer of energy from the surface to underlying layers. However, they are not responsible for any higher energy uptake at the surface, which is exclusively triggered by increasing energy flux from the atmosphere due to SEB changes when surface temperatures reach 0 C. The resulting enhanced energy uptake makes cold accumulation zones very vulnerable to air temperature rise.


Introduction
The thermal regime of glaciers needs to be modeled to study the impact of climate change on ice flow and intraglacial or subglacial hydrology.Indeed, englacial temperatures control ice viscosity and basal sliding (Paterson, 1994).They also affect the drainage system (Ryser et al., 2013;Gilbert et al., 2012;Skidmore and Sharp, 1999).In addition, most glacier hazard studies require a thorough analysis of the thermal regime of glaciers (Failletaz et al., 2011;Gilbert et al., 2012;Huggel et al., 2004;Haeberli et al., 1989).Although several mathematical models have already been developed to simulate energy transfer within glaciers (Wilson and Flowers, 2013;Aschwanden and Blatter, 2009;Zwinger and Moore, 2009;Lüthi and Funk, 2001;Funk et al., 1994;Hutter, 1982), the boundary conditions of these models are often not physically based and not related to external meteorological data, making future transient simulations impossible.This is due to the fact that applied integration time steps are generally longer than those of surface process timescales such as the percolation/refreezing of surface meltwater.In Wilson and Flowers (2013), volumetric heat flux due to meltwater refreezing is calculated using a degree-day model and released in a given firn layer immediately below the surface.This thickness is however determined arbitrarily and thus depends on the final result.Another approach, used by Zwinger and Moore (2009), is the Wright Pmax model (Wright et al., 2007).This model quantifies meltwater refreezing heat flux on the basis of englacial temperature measurements in the active layer.Although this formulation gives good results and can be used at an annual time step, it requires monthly englacial temperature measurements that are very rare.Even more importantly, it cannot be used for future simulations.Furthermore, investigation of the spatial variability of meltwater refreezing heat flux would require extensive temperature measurements, usually not available in the field.Other studies are limited to steady state simulation and do not focus on the relationship between climate and boundary conditions (Aschwanden and Blatter, 2009;Hutter, 1982) or use only time-dependent surface temperature changes imposed by the Dirichlet condition (Lüthi and Funk, 2001;Funk et al., 1993).
The aim of this study is to develop a model allowing longterm past and future glacier thermal regime simulations by explicitly taking into account the main surface firn processes such as meltwater percolation and refreezing.To achieve this goal, we could have coupled a thermal regime model to a sophisticated physically based snow model (e.g., CRO-CUS, Brun et al., 1989Brun et al., , 1992)).However, this approach requires a considerable amount of meteorological data at a short timescale (hour), which is not appropriate for simulations to be performed over several decades or centuries.This is why we decided to develop a simple surface temperature model suitable for long-term thermal regime simulation based only on daily air temperatures and surface topography parameters.The study site and data are presented in Sect. 2. Numerical models are described in Sect. 3 and results are shown and discussed in Sect. 4. Section 4 also explores the spatial variability of surface and englacial temperatures and applies the model to an example involving deep borehole temperature profiles.The last section presents our conclusions and some future perspectives for this work.
2 Study site and data

Study site
Col du Dôme is located in the Mont Blanc area at an elevation of 4250 m above sea level (a.s.l.).This is a cold accumulation zone (Vincent et al., 2007a) on a saddle with slopes of various aspects (Fig. 1).Snow accumulation rates range from 0.5 m w.e.yr −1 (meters of water equivalent per year) to 3.5 m w.e.yr −1 over short distances (a few hundred meters) and horizontal surface velocities do not exceed 10 m yr −1 (Vincent et al., 2007b).The mass balance seems to have been stable over the last one hundred years (Vincent et al., 2007b).Mean firn temperature is about −10 • C (below the active thermal layer) but has been significantly rising over the last 20 years due to a recent increase in regional air temperature and surface melting (Gilbert and Vincent, 2013).These temperature changes are spatially very variable and dependent on local firn advection velocity, slope, aspect and basal heat flux.

Meteorological data and firn temperatures
An automatic weather station (AWS) located near the center of the saddle (Fig. 1) ran continuously between 3 July and 23 October 2012.The measurements were carried out within the surface boundary layer.Wind speed, air temperature, humidity, incident and reflected short-wave radiation and incoming and outgoing long-wave radiation were recorded as half-hourly means of measurements made every 10 s (see Table 1).Instantaneous values of surface position and wind direction were collected every half hour.The Vaisala hygrothermometer was artificially ventilated in the daytime to prevent measurement errors due to radiation.Five meters from the AWS, 16 thermistors (PT100) were set up in the firn at depths of 0.11 m, 0.24 m, 0.34 m, 0.44 m, 0.65 m, 1 m, 2 m, 3 m, 5 m, 8 m, 12 m and 16 m and at heights of 0.15 m, 0.30 m, 0.45 m and 0.60 m above the surface from 3 July 2012 to 13 June 2013.The purpose of the sensors initially located in the air was to measure the subsurface temperature in the event of snow accumulation.Firn temperature was recorded as half-hourly means of measurements made every minute.Due to surface melting, the first three sensors were found at the same depth at the end of the melting period, making the first 40 cm deep temperature measurements not reliable after September.Air and surface temperatures were also recorded at this location from 3 July 2012 to 13 June 2013.The characteristics and specifications of all sensors are summarized in Table 1.

Deep borehole temperature profiles
Englacial temperature measurements using a thermistor chain (see Table 1) were performed from surface to bedrock in seven boreholes drilled between 1994 and 2011 at three different sites located between 4240 and 4300 m a.s.l.(Fig. 1).Ice thicknesses were 40, 126, and 103 m at sites 1, 2, and 3, respectively (see Gilbert and Vincent (2013)

Near-surface densities
On 27 September 2011, density profiles were measured down to a depth of 4 m at 14 different drilling locations at Col du Dôme using a manual auger device (Fig. 1, blue dots).

Modeling approach
Two distinct one-dimensional and spatially distributed models are used to calculate the near-surface firn temperatures at Col du Dôme.The first is used to calculate the firn surface temperature and melting from a surface energy balance model using the meteorological data recorded by the AWS.
The second is a heat flow model coupled to a water percola-tion model with temperature and melting rate at the surface as input data.This model is used successively with different input data sets obtained first from the energy balance model and then from a parameterized approach.

Model 1: surface energy balance
A surface energy balance (SEB) model coupled to a onedimensional heat flow model was developed in order to determine the energy fluxes within the firn during the measurement period.When heat added by precipitation and penetration of solar radiation are neglected, the SEB equation can be written as (fluxes directed toward the surface are considered positive) (e.g., Oke, 1987) where R is the radiative net balance (W m −2 ), H the turbulent sensible heat flux (W m −2 ), LE the turbulent latent heat flux (W m −2 ), Q the energy flux within the firn (W m −2 ) and Q m the energy flux available for melting.The radiative heat balance is calculated as where S w ↓ is the incident short-wave radiation, S w ↑ the reflected short-wave radiation, L w ↓ the incoming longwave radiation and σ T emitted by the surface (calculated using the modeled surface temperature T s and the black body law with σ = 5.67 × 10 −8 W m −2 K −4 , Essery and Etchevers, 2004).Turbulent fluxes are explicitly calculated using the bulk aerodynamic approach, including stability corrections of the surface boundary layer (Essery and Etchevers, 2004).Scalar and momentum roughness lengths z 0 are assumed to be equal.Since z 0 has been tuned to match the variation of the internal energy variation (see Sect. 4.2), selecting a single value for z 0 or using distinct values for every roughness length (momentum, humidity and temperature) would have changed the values of the roughness parameters, but would not have changed the final results of our simulations.As a consequence, z 0 should be considered more as a tuning parameter than as a true roughness length.Sensible and latent heat flux are calculated from where ρ air and C p air are the density (kg m −3 ) and heat capacity (J K −1 kg −1 ) of air, respectively, Q sat (T s , P s ) is the saturation humidity at snow surface temperature T s and pressure P s , C H is a surface exchange coefficient and L s is the latent heat of sublimation (J kg −1 ).T 1 , Q 1 and U 1 are, respectively, the air temperature (K), the specific humidity (kg kg −1 ) and the wind speed (m s −1 ) at level z 1 (m) above the surface.The exchange coefficient is calculated as a function of the atmospheric stability which is characterized by the bulk Richardson number where ε is the ratio of molecular weights for water and dry air and g the gravitational acceleration (m s −2 ).
Then we have the following (Louis, 1979): where C Hn is the neutral exchange coefficient for roughness lengh z 0 and f h a correction factor for atmospheric stability.k m is the Von Karman constant.
The SEB model is coupled to a subsurface model to determine Q (Eq.1).The one-dimensional heat equation is therefore solved over a 16 m-deep temperature profile.The upper surface boundary condition is determined by the SEB as a flux condition (Neumann) and we assume no heat flux at 16 m depth over the whole simulation period (125 days) (zero-flux boundary condition).This assumption is supported by the fact that over such a short simulation period, basal heat flux has no influence on the modeled temperatures above 10 m depth and consequently it does not have any influence on our results.The time step is set to 5 min (half-hourly meteorological data are linearly interpolated for each time step) and vertical resolution is 4 cm.The heat equation is solved using the Crank-Nicholson scheme.The heat capacity of ice is set to 2050 J K −1 kg −1 (Paterson, 1994).Densities have been measured from 0 to 16 m depth into the borehole drilled to install the thermistors.The heat conductivity is calculated from density data using the relation from Calonne et al. (2011).The initial temperature profile is taken to be similar to the profile measured at the beginning of the simulation (3 July 2012).Solar radiation penetration is accounted for assuming that short-wave radiation exponentially decreases as a function of depth (Colbeck, 1989): where F sw is the radiative flux at depth z and δ (m) is the characteristic length of penetration (m).
During the simulation, at every time step, if the surface temperature T s exceeds 0 • C, T s is systematically reset to 0 • C, and the temperature difference is used to calculate melt energy and the resulting volume of meltwater (Hoffman et al., 2008).This amount of energy is then released into the first underlying cold layer below the surface to simulate water percolation and refreezing.Irreducible water saturation (see Sect. 3.2) is taken into account in each model layer.Indeed, when the water content within a model layer exceeds the irreducible water saturation, the remaining water is allowed to reach the neighboring underlying layer.All parameters and variables of model 1 are summarized in Table 2.

Model 2: coupled water percolation and heat transfer model
In this second model, Dirichlet conditions are assumed at the surface for temperature and the surface boundary condition for percolation is treated as a water flux.Because the vertical resolution of the second model is coarser than that of the previous one (10 cm instead of 4 cm), we assume that the first surface layer is thick enough to absorb all the solar radiation.This assumption allows us to take into account the solar radiation penetration, which is here converted into a temperature rise of the surface layer only.In this way, we add a temperature T rad (K) to T s .Water is assumed to percolate through homogeneous snow.In this way we solve for water saturation in snow using gravity flow theory (Colbeck and Davidson, 1973) adapted for dry and cold snow.All symbols used in the following equations (Eqs.12-20) are explained in Table 3 together with their respective units and values when Table 2. Parameters and variables used in the SEB (model 1), with their respective values when available.

Symbol Values and units
Ice heat capacity c p 2050 J kg −1 K −1 Air heat capacity c pair 1.005 × 10 3 J kg −1 K −1 Solar radiation characteristic length of penetration δ 2.5 × 10 −2 m Ratio of molecular weights for water and dry air available.We define the effective water saturation S * : where S is the water saturation in the snow, and S r is the irreducible water saturation that is permanently retained by capillary forces.If S < S r there is no water flow and if S > S r , S * gravitationally advects and we have where is the porosity, n a constant set to 3.3 (Colbeck and Davidson, 1973), ρ w the water density (kg m −3 ), g acceleration due to gravity (m s −2 ), µ the viscosity of water (kg m −1 s −1 ), K the intrinsic snow permeability (m 2 ) and R a negative source term coming from liquid water refreezing (s −1 ).Snow permeability is calculated as a function of firn density ρ f and mean grain size d (m) using the relationship from Shimizu (1970): The heat advection/diffusion equation is solved and coupled to the water saturation at each time step: where z is the depth (m), t the time (s), T the firn temperature (K), c p the heat capacity of ice (J kg −1 K −1 ), ρ f the firn density (kg m −3 ), v z the vertical advection velocity (m s −1 ), k the thermal conductivity of firn (W m −1 K −1 ), Q lat the latent heat released by refreezing meltwater (W m −3 ) per unit of time and volume, m the mass of water that freezes per unit volume (kg m −3 ) during the time increment t (s) and L the latent heat of fusion (J kg −1 ).Density variations due to meltwater refreezing are neglected.Indeed, melting events at such elevation are very rare (not exceeding 10 days per year) and are always punctuated by snowfall events.These lead to a very small density increase.As proposed by Illangasekare et al. (1990), we impose that m is only a fraction of the maximum water freezing (m max ) allowed by the conservation of heat because in general the time step used in modeling will be small compared to the velocity of freezing processes, so we define where w is the snow water content (kg m −3 ), T 0 the ice melting point (273.15K), dt the time step (s) and τ f the freezing calibration constant (s −1 , τ f < 1 / dt).
The model is run in 30 min time steps.We assume a zero flux boundary condition at 16 m depth over the entire simulation period (125 days).Actually, this flux at 16 m is not insignificant (approximately of 4.0 × 10 −1 W m −2 ), but sensitivity tests performed with various constant fluxes at 16 m have shown that assuming a zero flux at 16 m does not change the modeled temperature above 10 m depth.The initial temperature profile is taken to be similar to the profile measured at the beginning of the measurement period (3 July 2012) and the density profile is set to the one measured on 3 July 2012.This density profile is assumed to be constant over the whole simulation period.The problem is solved numerically using the Elmer/ice finite element model (Gagliardini et al., 2013) based on the Elmer open-source multi-physics package (see http://www.csc.fi/elmer for details) that will make it possible to work in three dimensions and perform thermo-mechanical coupling easily in future studies.Vertical resolution varies from 10 cm near the surface to 40 cm at 16 m depth.

Meteorological conditions in summer 2012
AWS half-hourly measurements from 3 July 2012 to 23 October 2012 are reported in Fig. 2 for wind speed, incoming short-and long-wave radiation, relative humidity, air temperature and firn surface temperature.These measurements show 13 days of positive air temperature events, with the warmest event from 15 to 20 August.There were 42 cloudy days and 69 clear sky days.Based on surface elevation measurements, we estimate that three snow falls occurred on 29 August (16 cm), 11 September (30 cm) and 27-28 September (15 cm) (Fig. 3).Mean wind speed was 5.4 m s −1 , mainly from the southwest and west (> 50 %) with some strong wind events from the northeast.Wind speed and wind direction were unavailable during two periods: 12 July to 6 August and 7-12 October, representing 27 % of the whole measuring period (Fig. 2).Wind speed was therefore set to its mean value during these periods for energy balance modeling.This approximation is supported by the fact that our simulated surface temperatures during this period remain in good agreement with measurements.

Surface energy balance (model 1)
The calculated temperature and surface melt obtained from the SEB model are reported in Fig. 3. Given that all meltwater refreezes within the cold firn pack, energy is conserved and the modeled energy input should match the energy content variation of the firn pack.For this purpose, the roughness length for momentum z 0 is tuned so that the measured firn internal energy matches the modeled energy input (Fig. 4) for every simulation.comes from turbulent fluxes.The values of radiation penetration characteristic length δ and water residual saturation S r are chosen to minimize RMSE and bias between modeled and measured surface and subsurface temperatures.We find δ = 0.025 m that corresponds to the value reported by Brandt and Warren (1993) and S r = 0.005, in good agreement with the model 2 study, where water percolation is investigated in more detail (see Sect. 4.3).The corresponding value of z 0 that best matches the measured firn internal energy is 4.0 mm.This value is close to the expected value found in the literature at high-altitude cold sites like in Brock et al. (2006) or Wagnon et al. (2003), where roughness lengths have all been chosen equal and also used as tuning parameters.
The calculated hourly surface melt and corresponding cumulative melt are displayed in Fig. 3d together with the measured distance between the surface and the SR50 ultrasound sensor expressed as a surface height change in m w.e. using a constant surface density (set to 380 kg m −3 ) and corrected according to local air temperature.Two main melting events (25 July to 4 August, and 16 to 26 August) can be identi-fied and agree fairly well with the surface height variation measurements.The mismatch observed between cumulative surface melt and surface height variations after 1 September may be attributed to wind erosion and fresh snow settlement.Modeled surface temperatures also agree well (r 2 = 0.87, RMSE = 1.6 K, bias = −0.02K, n = 5326 half hours) with measured surface temperatures inferred from long-wave radiation emissions (Fig. 3a), further supporting the ability of the model to simulate surface melting efficiently.Fig. 3a, b and c compare measured (red line) and modeled surface temperatures at 0, 11 and 24 cm depth, respectively, accounting for solar radiation penetration (black line) or not (blue line).When no melting occurs, the cold bias between measurement and model is efficiently attenuated by taking into account solar radiation penetration, showing that this effect cannot be neglected at this high elevation cold site, like in cold firn in Greenland (Kuipers Munneke et al., 2009).
Figure 4a compares the integrated thermal energy of the firn pack (red line) obtained from temperature measurements with the modeled energy input (E b , blue line) obtained from the model and the modeled firn thermal energy.Firn thermal energy E firn is calculated from where D(m) is the temperature measurement depth, which is constant over the simulation period, i.e., 16 m at Col du Dôme. Figure 4a shows that energy is conserved except during the second melting event.In this case, the energy input calculated from the SEB exceeds the firn pack thermal energy.Indeed, part of the energy transferred to the firn pack from the surface is stored as latent heat because some firn layers have reached the melting point (Fig. 3e).Consequently, liquid water is stored in the firn (Fig. 4b, blue line).This energy is released when the water refreezes from 17 August to 6 September.This explains why the modeled energy input is once again equal to the firn pack thermal energy several days after the melting event (Fig. 4a).
As illustrated in Fig. 4a, each melting event results in a large energy increase in firn.This means that the surface energy balance is modified during these events.Figure 5 focuses on two consecutive days of melting (18 and 19 August).During these two days, we compare the energy flux of two cases: (i) a virtual case for which the surface temperature can exceed 0 • C and no melting can occur (dashed lines) and (ii) the true case for which the surface temperature cannot exceed 0 • C and melting occurs.Note that case (i) is a virtual case where meteorological records have been considered unchanged even though surface temperature is allowed to exceed 0 • C. In reality, meteorological variables such as air temperature, air vapor pressure, etc. would be modified along with the turbulent fluxes if surface temperature rises above 0 • C. Nevertheless, we believe that this simple comparative  approach is useful for qualitatively understanding the impact of melting on the SEB over cold surfaces.
In the daytime, for case (i), the short-wave radiative balance is efficiently compensated by the other fluxes due to increasing surface temperature, implying on one hand an enhanced energy loss due to increased outgoing long-wave radiation and on the other hand unstable conditions in the surface boundary layer, leading to negative values for sensible heat flux and a strong negative latent heat flux.For case (ii), net short-wave radiation largely dominates the other fluxes because surface temperature cannot exceed 0 • C, thereby limiting heat loss through long-wave radiation and maintaining stable conditions inside the surface boundary layer that reduces turbulent fluxes.In this case, a large amount of energy is transferred to the firn.At night, energy fluxes remain unchanged for both cases.Consequently, the energy uptake during melting events is due to the fact that surface temperature is maintained at 0 • C by thermodynamic equilibrium between the liquid and solid phases.
We conclude that each melting event is associated with a significant energy transfer to the firn pack and the duration of the event therefore has a very strong impact on the total energy balance of the firn pack during summer.With expected higher air temperatures, melting events will become more frequent and the energy will be transferred more efficiently to the firn in the accumulation zones of glaciers.In our case, the strong energy uptake during melting events is less due to melt water percolation and refreezing than to an energy gain in the SEB due to peculiar conditions of the lower atmosphere-firn surface continuum triggered by a 0 • C firn surface.In other words, for similar atmospheric conditions, firn surfaces are able to absorb more energy when the surface temperature is at 0 • C than when it is negative, explaining why warming of cold accumulation zones of glaciers is more efficient when melting conditions are encountered than when they are not.

Subsurface temperature and water content
(model 2)

Comparison with data
Measured subsurface temperatures during summer 2012 are shown in Fig. 3e.The influence of surface melting events is clearly visible in the firn temperature data (Fig. 3e).Indeed, we observe a step change in the time evolution of firn temperature between 2 and 5 m deep on 25-27 July and 16-18 August.These temperature increases at these dates are too sharp and rapid to come from diffusive processes and are likely due to additional energy supplied by refreezing meltwater.These observations reveal two striking features: (i) water percolates into cold firn until 4 to 5 m deep and (ii) liquid water crosses the cold firn without totally refreezing.Model 2 uses the surface temperature and surface melting flux calculated by the SEB model as input data.The only free parameters of the model are the percolation parameters that are not constrained (S r and τ f ) and the constant temperature correction T rad of surface temperature accounting for heating due to solar radiation penetration.These parameters were adjusted to τ f = 2.0 × 10 −5 s −1 , S r = 0.005 and T rad = 2.0 K, respectively, to match the measured temperatures at all depths (Fig. 6a-e).Figure 4b shows that the modeled water content using model 2 agrees well with the expected value based on the energy balance (see Sect. 4.1).However, the value of S r is one order of magnitude lower than past published values (S r = 0.03 to 0.07) (Illangasekare et al., 1990).We think that water does not percolate within the cold firn in a uniform manner but more locally using small preferential pathways (Harrigton et al., 1996), explaining why less water is retained by capillarity in this case.Higher values of S r lead to more liquid water being stored within the first meters of firn below the surface, which prevents water from percolating deep enough to explain the observed firn temperatures.The freezing calibration constant τ f makes it possible to simulate water percolation in cold firn where only part of the liquid water refreezes between two time steps (30 min here).This parameter is well constrained by our temperature measurements.Indeed, if the value of τ f is too high, the water will not manage to percolate through cold firn and will not influence the temperature field deep enough compared to measurements.Conversely, if the value of τ f is too low, the firn temperature never reaches 0 • C and the energy released by meltwater refreezing is distributed over too large a thickness.The value of T rad shows that it is not possible to simulate subsurface temperature in snow from surface temperature measurements only, because solar radiation partly penetrates into the firnpack.Using a constant correction on T s seems to give satisfactory results.
From these results, we can conclude that our subsurface firn temperature model is able to reproduce the water content and the subsurface temperature field accurately (Figs.4b  and 7).However, this model cannot be applied for simulations over several decades or centuries because the halfhourly data it requires are not available over such long periods.Consequently, a simplified approach has been developed with boundary conditions parameterized using daily air temperature data.

Simplified approach
For numerical reasons, in order to use a daily time step, the vertical resolution has been reduced from a few centimeters to about 50 cm.In addition, the previous percolation approach is not meaningful at this space and timescale.For this reason, we now use a constant percolation velocity.The value of 1.0 × 10 −6 m s −1 provides the best agreement between the results obtained with this simplified approach and those with the previous model using the Colbeck and Davidson formulation (Fig. 6).All other parameters remained unchanged.
In this way, the boundary conditions are now parameterized using the daily mean and maximum air temperatures.Daily surface temperature is assumed to vary in the same way as daily mean air temperature.This is confirmed by one year of hourly simultaneous measurements of surface (infra-red camera) and air temperature at this site (Fig. 7).From these measurements, the mean difference between daily mean air and surface temperature is estimated at 3.4 K (Fig. 7).In order to take into account the enhanced energy uptake during melting events, melting is assumed to occur when daily maximum air temperature reaches 0 • C. The daily surface melt flux is calculated using the following degree-day model (Hock, 1999): where M is the daily amount of surface melt (m w.e.d −1 ), a the melt factor (m w.e.d −1 K −1 ), T 0 the melting point (273.15K) and T max the daily maximum air temperature.
Air temperature data from Lyon-Bron meteorological station located ∼ 200 km west of the studied site were selected given that it is one of the longest meteorological series in this region (Gilbert et al., 2012).Comparison between Lyon-Bron air temperature and on-site surface temperature at daily time steps between 3 July 2012 and 13 June 2013 leads to an altitudinal gradient between Lyon air temperature and Col du Dôme surface firn surface temperature of −5.94 × 10 −3 K m −1 .The altitudinal gradient varies significantly between summer and winter (−5.6 × 10 −3 K m −1 in winter and −6.5 × 10 −3 K m −1 in summer).As shown previously, solar radiation penetration cannot be neglected to model firn temperature.As a consequence, the altitudinal gradient selected (or tuned) to reconstruct surface temperature explicitly includes this effect.The constant mean value of −5.9 × 10 −3 K m −1 allows us to match modeled and measured subsurface temperatures.The melt factor is adjusted to simulate the total melt calculated from the SEB during summer 2012 (3 July 2012-23 October 2012) and set to 3.3 × 10 −4 m w.e.d −1 K −1 .Calculated firn temperatures using this simplified approach agree well with in situ measurements (Fig. 6c, d, f) and properly account for the step change observed during melting events.This reveals that the simple model using daily temperature data from a remote station provides satisfactory results and can be used to simulate long-term firn subsurface temperature variations in a cold accumulation zone.

Melting spatial variability
In order to quantify melting spatial variability, 14 4 m-deep density profiles were measured in the field on 27 September 2011 (Figs. 8 and 9).The amount of meltwater was quantified using the density anomaly related to meltwater refreez- ing in the firn.This anomaly was quantified in every section of the firn core by comparing the measured density with a computed density obtained from an empirical firn densification model (Herron and Langway, 1980) (Fig. 9, blue line).Gilbert et al. (2010) successfully applied this method to a high-altitude site in the Andes to quantify local melting.Uncertainty in melt quantification (error bars in Fig. 11) was computed considering uncertainty from mass, length and circumference measurements of each core section.This Col du Dôme site has been monitored for snow accumulation since 1994 using a well-distributed stake network (Vincent et al., 2007b) that provides a good assessment of the spatial variability of snow accumulation at Col du Dôme.In this way, we were able to extrapolate spatially the snow accumulation monitored since 2010 in the vicinity of the measured density profiles (Fig. 8, red stars).Knowing the accumulation rate, we were then able to date every density profile accurately (Fig. 9, dashed lines).Only the density anomaly above the horizon identified on 28 October 2010 has been considered here to make sure that we only take into account the refreezing meltwater of summer 2011.The quantified melt at each site is plotted in Fig. 10 as a function of the mean annual incoming potential solar radiation (PSR) of the corresponding site, taking into account the shading effect of the surrounding relief and atmospheric transmissivity (Fig. 8).Atmospheric transmissivity was determined by comparison of measured short-wave radiation during clear sky conditions with calculated theoretical potential solar radiation and set to 0.78.Although uncertainties related to these measurements are high, we observe a good correlation between PSR and surface melt (Fig. 10 Herron and Langway (1980).
melt values obtained from the density anomaly method at site 8, 5 m from the AWS.In order to model surface melt over the whole domain using the SEB model, the measured incoming short-wave radiation at AWS is varied artificially according to PSR before re-running the SEB model.In this way we obtain a relationship between surface melt and PSR for summer 2011 from the SEB model.A quadratic function provides a reasonable fit for the melt-PSR relationship (Fig. 10).The non-linearity can be explained by the fact that increasing short-wave radiation enhances the frequency and the duration of melting events, which in turn shifts the energy balance towards positive values (see Sect. 4.1).
In order to take into account the effect of PSR on the melting intensity in the degree-day formulation, we revise Eq. ( 22).PSR is now taken into account in Eq. ( 23) proposed by Hock (1999) and melt is calculated by where a PSR is the melt factor as a function of PSR.
Using daily maximum air temperature inferred from Lyon-Bron daily maximum air temperature (using the same gradient as the one obtained between Lyon-Bron air temperature and Col du Dôme surface temperatures in 2012-13) and melt calculated using SEB in summer 2011, we calculate a PSR for different PSR values.We found a quadratic relation-ship (Fig. 10): In this way, this relationship can be applied to calculate the surface melt and the firn temperature can be calculated using the above-described simplified model for the whole area.However, this relationship is not likely to be transferable to another site.Indeed, melt factors depend on site characteristics that influence the surface energy balance such as mean albedo, wind speed, humidity and surface roughness.Thus, the relationship between melt factor and PSR needs to be recalibrated for every studied site.

Application of the simplified model at multi-decennial scale to reconstruct deep borehole temperature profiles
The locations of borehole sites 1, 2 and 3 are indicated on the map in Fig. 1 and their respective annual PSR values are assessed at 220, 220 and 192 W m −2 .Our simplified 1-D model (model 2 using daily temperatures from Lyon-Bron station and topographic parameters) was applied to the three sites and run over the period 1907-2012.The vertical advection profile is calculated from the measured surface advection velocity and assumed to vary linearly with depth (Gilbert and Vincent, 2013).Basal heat flux is set to 15 × 10 −3 W m −2 for sites 1 and 2 and set to 30 × 10 −3 W m −2 for site 3 according to measured basal temperature gradients.The basal heat flux is specified at 150 m depth in the bedrock because it can be considered to be approximately constant at this depth for the timescale of the simulation.The comparison with simulations using a basal heat flux specified at 800 m depth shows that this assumption only influences modeled glacier temperatures below 100 m depth and with a temperature difference that does not exceed 0.2 K.The thermal properties of the rock (gneiss and granite) are taken from Lüthi and Funk (2001) with a thermal conductivity of 3.2 W m K −1 m −1 , a heat capacity of 7.5 × 10 2 J kg −1 K −1 and a density of 2.8 × 10 3 kg m −3 .The initial temperature profile in 1907 is calculated as a steady state temperature profile.Constant surface steady state temperature for each site is given by Gilbert and Vincent (2013).
For each site, the altitudinal temperature gradient is adjusted to match observed borehole temperatures.We found values of 5.95 × 10 −3 , 5.70 × 10 −3 and 5.83 × 10 −3 K m −1 for sites 1, 2 and 3 respectively.The respective melt factors for sites 1, 2 and 3 are derived from PSR values (220, 220 and 192 W m −2 ) and equal to 3.5 × 10 −4 , 3.5 × 10 −4 , and 2.0 × 10 −4 m w.e.d −1 K −1 .The results plotted in Fig. 11 show that temperature differences between the three boreholes can be explained by differences in both surface melt (in agreement with PSR differences) and vertical advection velocities.As already seen, PSR has a strong influence on firn temperature of a cold glacier mainly through surface melt and must be accounted for to simulate the thermal regime of cold glaciers.Indeed, site 2 experiences higher surface melting rates than site 3, leading to a stronger firn temperature rise over the first decade of this century.Differences between the two sites are amplified by stronger vertical advection velocities at site 2. The use of the melt factor specified on the basis of the PSR relationship leads to modeled temperature in good agreement with measurements for sites 1 and 2 but also to an englacial warming overestimation of about 0.5 • C for site 3, which is acceptable for thermal regime simulation at the glacier scale.This overestimation at site 3 could be corrected by the use of melting factor values of 1.0 × 10 −4 mw.e. d −1 K −1 instead of 2.0 × 10 −4 mw.e. d −1 K −1 (Fig. 11, blue dashed line).

Conclusions
To simulate the transient thermal regime of a cold glacier, we developed a model based on a surface energy balance and water percolation parameterization.The results agree well with in situ firn temperature measurements for summer 2012.
Our energy balance study highlights the impact of the duration of melting events on cold firn temperature.Once the surface temperature reaches 0 accumulation zones become extremely sensitive to climate change.Indeed, when surface temperature reaches 0 • C, the sum of the radiative and turbulent fluxes shifts towards more positive values, which means that more energy is transferred to underlying layers through water percolating and then refreezing within cold layers.Consequently, air temperature rise has a very strong impact on temperature profiles in cold glaciers by increasing the frequency and duration of melting events at high elevations.Note that water percolation and refreezing are very efficient processes for energy transfer from the surface to subsurface firn layers; however, they are not responsible for a higher energy uptake at the surface, which is mainly due to the fact that surface temperature is limited to 0 • C. Our study also highlights the influence of solar radiation penetration on firn temperature.Indeed, using measured surface temperature as a Dirichlet condition amounts to neglecting radiation penetration, which leads to significant underestimation of sub-surface temperature.
To perform numerical simulations over several centuries, instead of applying a sophisticated surface energy balance requiring an excessive amount of data, we also developed a simplified approach based only on daily air temperature and topography.Results from this simplified approach agree well with in situ firn temperature measurements and with results coming from SEB modeling.In addition, if the melt factor is known, this simplified model allows us to reconstruct deep borehole temperature profiles that agree well with measurements.
Our measurements show that the spatial variability of melting is highly dependent on the potential solar radiation.This confirms the results obtained by Suter (2002) at these high elevations and we propose a relationship between the melt factor and potential solar radiation at Col du Dôme.
The study of water percolation into the firn shows that gravity flow theory (Colbeck and Davidson, 1973) is sufficient to reproduce the observed subsurface temperature by taking into account water flow and refreezing.However, residual saturation in firn must be set to a low value given that water does not percolate uniformly into the cold firn.This may be due to the formation of impermeable ice layers capable of driving water through preferential pathways, leaving some parts of the firn pack dry.This results in an apparent residual saturation much lower than expected.However, the gravity flow approach is meaningless at a daily timescale and the use of a constant velocity flow is recommended instead.Although this water percolation scheme is far from reality, it provides good results for subsurface temperature modeling.The same model could be applied in temperate accumulation zones; however, it would be necessary to determine how the water drains at the firn-ice transition.
As the climate is expected to change in the future (IPCC, 2007), cold glacier temperatures will be modified.The response of subsurface firn temperature to air temperature rise will be largely amplified by an increase in the duration and frequency of melt events.This will lead to strong changes in ice temperature fields.Coupling our surface model with a thermo-mechanical model will make it possible to study the transient response of cold glaciers to climate change and investigate glacier hazards related to thermal regime changes, such as cold hanging glacier stability.

Fig. 1 .
Fig. 1.Map of Col du Dôme (Mt Blanc range, France) showing the locations of the automatic weather station (AWS; blue square), density measurement sites (blue dots) and borehole drilling sites (stars).

Fig. 2 .
Fig. 2. Meteorological data recorded by the AWS in summer 2012 and surface temperature calculated using measured infra-red surface emissions.

Fig. 5 .
Fig. 5. Energy fluxes during two consecutive days of melting (18 and 19 August).Comparison between two cases: (i) no phase change is taken into account and firn temperature can artificially exceed 0 • C (dashed line); (ii) melting is taken into account (lines and filled blue area).

Fig. 6 .
Fig. 6.Measured (red line) and modeled (blue and black lines) firn temperature at different depths (a, b, c, d).Modeled firn subsurface temperature during summer 2012 using the energy balance model at a 30 min time step (e, blue line in a, b, c, d) compared to modeled temperature using the Lyon-Bron temperature at a daily timescale (f, black line in c, d).

Fig. 7 .
Fig. 7. Comparison between daily measured air (T air ) and surface (T s ) temperature between 3 July 2012 and 13 June 2013 at Col du Dôme.

Fig. 8 .
Fig. 8. Incoming potential solar radiation at Col du Dôme (color scale).Numbers and blue dots indicate the locations of sites where density measurements were performed.Red stars indicate the locations of sites where snow accumulation has been monitored since 2010.Black square is the location of the AWS.
). Between 27 June and 26 October 2011, the AWS ran continuously, allowing us to run the SEB model over the whole summer of 2011 and in turn quantify the corresponding on-site melt values.Results agree well with the A. Gilbert et al.: Modeling near-surface firn temperature in a cold accumulation zone

Fig. 10 .
Fig. 10.Melt obtained from density anomalies (black stars) as a function of potential solar radiation, compared to melt obtained using the surface energy balance model (red dashed line).Corresponding melt factor values in the degree-day model (Eq.24) are plotted in blue.

Fig. 11 .
Fig. 11.Measured (dashed lines and dots) and modeled (solid lines) temperature profiles at the three drilling sites (site 1 = red, site 2 = black, site 3 = green) for different dates.The blue dashed line is the temperature profile modeled at site 3 in 1999 and 2012 for a melt factor set to 1.0 × 10 −4 m w.e.d −1 K −1 .

Table 1 .
List of different sensors of AWS with their specifications.
1 Quantities are recorded as half-hourly means of measurements made every 10 s except for wind direction and accumulation/ablation, which are instantaneous values every 30 min. 2 Data measurements were interrupted from 12 July 2012 to 6 August 2012 and from 7 October 2012 to 12 October 2012 3 Measurement period was 3 July 2012 to 13 June 2013 4 Sensor depths are described in Sect.2.2.1.

Table 3 .
Parameters and variables used in the percolation and heat transport model (model 2), with their respective values when available.
Density profiles (red line) measured on 27 September 2011.Dashed lines correspond to different surface horizons identified at different dates.Blue lines are reference density profiles determined using the empirical model proposed by