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

A. Gilbert, C. Vincent, D. Six, Patrick Wagnon, L. Piard, Patrick Ginot

► To cite this version:

A. Gilbert, C. Vincent, D. Six, Patrick Wagnon, L. Piard, et al. Modeling near-surface firn temperature in a cold accumulation zone (Col du Dôme, French Alps): from a physical to a semi-parameterized approach. The Cryosphere, 2014, 8 (2), pp.689-703. 10.5194/tc-8-689-2014. hal-04704190

HAL Id: hal-04704190

https://hal.science/hal-04704190v1

Submitted on 23 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

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

A. Gilbert<fnref n="1" />,<fnref n="2" />, C. Vincent<fnref n="1" />,<fnref n="2" />, D. Six<fnref n="1" />,<fnref n="2" />, P. Wagnon<fnref n="2" />,<fnref n="3" />,<fnref n="4" />,<fnref n="5" />, L. Piard<fnref n="1" />,<fnref n="2" />, and P. Ginot<fnref n="2" />,<fnref n="3" />,<fnref n="4" />

*Correspondence to:* A. Gilbert (gilbert@lggeobs.ujf-grenoble.fr)

Received: 11 October 2013 – Published in The Cryosphere Discuss.: 19 November 2013

Revised: 11 February 2014 – Accepted: 4 March 2014 – Published: 17 April 2014

**Abstract.** 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.

## 1 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

<footnote n="1">¹CNRS, LGGE (UMR5183), 38041 Grenoble, France</footnote>

<footnote n="2">²Univ. Grenoble Alpes, LGGE (UMR5183), 38041 Grenoble, France</footnote>

<footnote n="3">³IRD, LGGE (UMR5183), 38041 Grenoble, France</footnote>

<footnote n="4">⁴IRD, LTHE (UMR5564), 38041 Grenoble, France</footnote>

<footnote n="5">⁵ICIMOD, GPO Box 3226, Kathmandu, Nepal</footnote>
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 melt-water 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 long-term 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., CROCUS, Brun et al., 1989, 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

## 2.1 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⁻¹ (meters of water equivalent per year) to 3.5 m w.e. yr⁻¹ over short distances (a few hundred meters) and horizontal surface velocities do not exceed 10 m yr⁻¹ (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.

<figure><img src="image_4.png" /><figcaption>**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).</figcaption></figure>

## 2.2 Field measurements

### 2.2.1 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.
<table><caption>**Table 1.** List of different sensors of AWS with their specifications.</caption><thead><tr><th>Quantity<sup>1</sup></th><th>Sensor type</th><th>Sensor height (m)</th><th>Accuracy according to manufacturer</th></tr></thead><tbody><tr><td>Air temperature, °C</td><td>Vaisala HMP155 – aspirated</td><td>1.5</td><td>±0.12 °C</td></tr><tr><td>Relative humidity, %</td><td>Vaisala HMP155 – aspirated</td><td>1.5</td><td>±1 % in [0–90 %]<br/>±1.7 % in [90–100 %]</td></tr><tr><td>Wind speed<sup>2</sup>, m s<sup>-1</sup></td><td>Young 05103</td><td>2.5</td><td>±0.3 m s<sup>-1</sup></td></tr><tr><td>Wind direction<sup>2</sup>, deg</td><td>Young 05103</td><td>2.5</td><td>±3 deg</td></tr><tr><td>Incident short-wave radiation, W m<sup>-2</sup></td><td>Kipp & Zonen CNR4<br/>0.3 &lt; 1 &lt; 2.8 µm</td><td>1.0</td><td>&lt; 5 % for daily sums</td></tr><tr><td>Reflected short-wave radiation, W m<sup>-2</sup></td><td>Kipp & Zonen CNR4<br/>0.3 &lt; 1 &lt; 2.8 µm</td><td>1.0</td><td>&lt; 5 % for daily sums</td></tr><tr><td>Incoming long-wave radiation, W m<sup>-2</sup></td><td>Kipp & Zonen CNR4<br/>4.5 &lt; 1 &lt; 42 µm</td><td>1.0</td><td>&lt; 10 % for daily sums</td></tr><tr><td>Outgoing long-wave radiation, W m<sup>-2</sup></td><td>Kipp & Zonen CNR4<br/>4.5 &lt; 1 &lt; 42 µm</td><td>1.0</td><td>&lt; 10 % for daily sums</td></tr><tr><td>Snow surface level, cm</td><td>Campbell Sci SR50</td><td></td><td>±1 cm</td></tr><tr><td>Firn temperature<sup>3</sup>, °C</td><td>PT100 (16 sensors<sup>4</sup>)</td><td></td><td>±0.05 °C</td></tr><tr><td>Surface temperature<sup>3</sup>, °C</td><td>Campbell Sci IR 120</td><td></td><td>±0.2 °C</td></tr><tr><td>Air temperature<sup>3</sup>, °C</td><td>PT100 – not aspirated</td><td>2.0</td><td>±0.05 °C</td></tr><tr><td>Deep borehole temperature, °C</td><td>YSI 44031 and PT100 (site 3, 1999)</td><td></td><td>±0.1 °C</td></tr></tbody></table>

### 2.2.2 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) for more details). Site 1 was measured in January 1999 and March 2012; site 2 in June 1994, April 2005 and March 2010; site 3 in January 1999 and March 2012.

### 2.2.3 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).

## 3 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 percolation

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.

### 3.1 Model 1: surface energy balance

A surface energy balance (SEB) model coupled to a one-dimensional 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)

$$R + H + LE + Q = Q_m \quad (\text{in W m}^{-2}), \qquad (1)$$

where R is the radiative net balance (W m⁻<fnref n="2" />), H the turbulent sensible heat flux (W m⁻<fnref n="2" />), LE the turbulent latent heat flux (W m⁻<fnref n="2" />), Q the energy flux within the firn (W m⁻<fnref n="2" />) and Qm the energy flux available for melting. The radiative heat balance is calculated as

$$R = S_w \downarrow + S_w \uparrow + L_w \downarrow - \sigma T_s^4 \quad (\text{in W m}^{-2}), \qquad (2)$$

where $S_w \downarrow$ is the incident short-wave radiation, $S_w \uparrow$ the reflected short-wave radiation, $L_w \downarrow$ the incoming long-wave radiation and $\sigma T_s^4$ the outgoing long-wave radiation

<footnote n="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.</footnote>

<footnote n="2">² Data measurements were interrupted from 12 July 2012 to 6 August 2012 and from 7 October 2012 to 12 October 2012</footnote>

<footnote n="3">³ Measurement period was 3 July 2012 to 13 June 2013</footnote>

<footnote n="4">⁴ Sensor depths are described in Sect. 2.2.1.</footnote>
emitted by the surface (calculated using the modeled surface temperature $T_s$ and the black body law with $\sigma = 5.67 \times 10^{-8} \text{ W m}^{-2} \text{ 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

$$H = \rho_{\text{air}} C_{\text{air}} C_{\text{H}} U_1 (T_1 - T_s) \quad (3)$$

$$LE = L_s \rho_{\text{air}} C_{\text{H}} U_1 (Q_1 - Q_{\text{sat}}(T_s, P_s)), \quad (4)$$

where $\rho_{\text{air}}$ and $C_{\text{air}}$ are the density (kg m$^{-3}$) and heat capacity (J K$^{-1}$ kg$^{-1}$) of air, respectively, $Q_{\text{sat}}(T_s, P_s)$ is the saturation humidity at snow surface temperature $T_s$ and pressure $P_s$, $C_{\text{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

$$R_{\text{ib}} = \frac{gz_1}{U_1^2} \times \left( \frac{T_1 - T_s}{T_1} + \frac{Q_1 - Q_{\text{sat}}(T_s, P_s)}{Q_1 + \varepsilon/(1 - \varepsilon)} \right), \quad (5)$$

where $\varepsilon$ 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):

$$C_{\text{H}} = f_{\text{h}} C_{\text{Hn}} \quad (6)$$

with

$$C_{\text{Hn}} = k_m \left( \ln \left( \frac{z_1}{z_0} \right) \right)^{-2} \quad (7)$$

$f_{\text{h}} =$

$$\left\{ \begin{array}{ll} (1 + 10 R_{\text{ib}})^{-1} & R_{\text{ib}} \ge 0(\text{stable}) \\ (1 - 10 R_{\text{ib}} \left( \frac{1 + 10 C_{\text{Hm}}(-R_{\text{ib}})^{1/2}}{f_{\text{z}}} \right))^{-1} & R_{\text{ib}} < 0(\text{unstable}) \end{array} \right. \quad (8)$$

$$f_z = 1/4 \left( \frac{z_0}{z_1} \right)^{1/2}, \quad (10)$$

where $C_{\text{Hn}}$ is the neutral exchange coefficient for roughness length $z_0$ and $f_{\text{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):

$$F_{\text{Sw}}(z) = (S_{\text{w}} \downarrow - S_{\text{w}} \uparrow) * \exp(-z/\delta), \quad (11)$$

where $F_{\text{sw}}$ is the radiative flux at depth $z$ and $\delta$ (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.

## 3.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_{\text{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><caption>**Table 2.** Parameters and variables used in the SEB (model 1), with their respective values when available.</caption><thead><tr><th></th><th>Symbol</th><th>Values and units</th></tr></thead><tbody><tr><td>Ice heat capacity</td><td>c<sub>p</sub></td><td>2050 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Air heat capacity</td><td>c<sub>pair</sub></td><td>1.005 &times; 10<sup>3</sup> J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Solar radiation characteristic length of penetration</td><td>&delta;</td><td>2.5 &times; 10<sup>-2</sup> m</td></tr><tr><td>Ratio of molecular weights for water and dry air</td><td>&epsilon;</td><td>0.622</td></tr><tr><td>Snow conductivity</td><td>k</td><td>W K<sup>-1</sup> m<sup>-1</sup></td></tr><tr><td>Von Karman constant</td><td>k<sub>m</sub></td><td>0.4</td></tr><tr><td>Latent heat of fusion</td><td>L</td><td>3.34 &times; 10<sup>5</sup> J kg<sup>-1</sup></td></tr><tr><td>Latent heat of sublimation</td><td>L<sub>s</sub></td><td>2.834 &times; 10<sup>6</sup> J kg<sup>-1</sup></td></tr><tr><td>Specific humidity at level z<sub>1</sub></td><td>Q<sub>1</sub></td><td>-</td></tr><tr><td>Saturation humidity</td><td>Q<sub>sat</sub></td><td>-</td></tr><tr><td>Air density</td><td>&rho;<sub>air</sub></td><td>7.9 &times; 10<sup>-1</sup> kg m<sup>-3</sup></td></tr><tr><td>Firn density</td><td>&rho;<sub>f</sub></td><td>kg m<sup>-3</sup></td></tr><tr><td>Residual water saturation</td><td>S<sub>r</sub></td><td>0.005</td></tr><tr><td>Firn temperature</td><td>T</td><td>K</td></tr><tr><td>Ice melting point</td><td>T<sub>0</sub></td><td>273.15 K</td></tr><tr><td>Air temperature at level z<sub>1</sub></td><td>T<sub>1</sub></td><td>K</td></tr><tr><td>Snow surface temperature</td><td>T<sub>s</sub></td><td>K</td></tr><tr><td>Wind speed at level z<sub>1</sub></td><td>U<sub>1</sub></td><td>m s<sup>-1</sup></td></tr><tr><td>Roughness length</td><td>z<sub>0</sub></td><td>4.0 &times; 10<sup>-3</sup> m</td></tr><tr><td>Level of measurements above the surface</td><td>z<sub>1</sub></td><td>1.5 m</td></tr><tr><td>Liquid water content</td><td>w</td><td>kg m<sup>-3</sup></td></tr></tbody></table>

available. We define the effective water saturation $S^*$:

$$S^* = \frac{S - S_r}{1 - S_r}, \qquad (12)$$

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

$$\phi(1 - S_r) \frac{\partial S^*}{\partial t} + n \rho_w g \mu^{-1} K S^{*(n-1)} \frac{\partial S^*}{\partial z} = R(\text{in s}^{-1}), \quad (13)$$

where $\Phi$ is the porosity, $n$ a constant set to 3.3 (Colbeck and Davidson, 1973), $\rho_w$ the water density (kg m$^{-3}$), $g$ acceleration due to gravity (m s$^{-2}$), $\mu$ 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 $\rho_f$ and mean grain size $d$ (m) using the relationship from Shimizu (1970):

$$K = 0.077 d^2 \exp(-7.8 \rho_f). \qquad (14)$$

The heat advection/diffusion equation is solved and coupled to the water saturation at each time step:

$$\rho_f c_p \left( \frac{\partial T}{\partial t} + v_z \frac{\partial T}{\partial z} \right) = \frac{\partial}{\partial z} \left( k \frac{\partial T}{\partial z} \right) + Q_{\text{lat}} \left( \text{in W m}^{-3} \right) \quad (15)$$

$$Q_{\text{lat}} = \frac{mL}{\Delta t} = R (L \rho_w \Phi (1 - S_r)) \left( \text{in W m}^{-3} \right), \quad (16)$$

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}$), $\rho_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_{\text{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 $\Delta 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_{\text{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

$$w = S \Phi \rho_w \quad (\text{in kg m}^{-3}) \qquad (17)$$

$$m_{\text{max}} = \rho_f c_p (T_0 - T)/L \quad (\text{in kg m}^{-3}) \qquad (18)$$

$$m = \begin{cases} w, & w < dt \tau_f m_{\text{max}} \\ dt \tau_f m_{\text{max}}, & w > dt \tau_f m_{\text{max}} \end{cases} \qquad (19)$$

where $w$ is the snow water content (kg m$^{-3}$), $T_0$ the ice melting point (273.15 K), $dt$ the time step (s) and $\tau_f$ the freezing calibration constant (s$^{-1}$, $\tau_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
<table><caption>**Table 3.** Parameters and variables used in the percolation and heat transport model (model 2), with their respective values when available.</caption><thead><tr><th></th><th>Symbol</th><th>Values and units</th></tr></thead><tbody><tr><td>Ice heat capacity</td><td>c<sub>p</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Snow mean grain size</td><td>d</td><td>1.0 × 10<sup>-3</sup> m</td></tr><tr><td>Acceleration due to gravity</td><td>g</td><td>9.81 m s<sup>-2</sup></td></tr><tr><td>Snow intrinsic permeability</td><td>K</td><td>m<sup>2</sup></td></tr><tr><td>Snow conductivity</td><td>k</td><td>W K<sup>-1</sup> m<sup>-1</sup></td></tr><tr><td>Latent heat of fusion</td><td>L</td><td>3.34 × 10<sup>5</sup> J kg<sup>-1</sup></td></tr><tr><td>Dynamic viscosity of water at 0 °C</td><td>μ</td><td>1.79 × 10<sup>-3</sup> kg m<sup>-1</sup> s<sup>-1</sup></td></tr><tr><td>Mass of refrozen water by unit of volume</td><td>m</td><td>kg m<sup>-3</sup></td></tr><tr><td>Maximum mass of water that can refreeze by unit of volume</td><td>m<sub>max</sub></td><td>kg m<sup>-3</sup></td></tr><tr><td>Water flow constant</td><td>n</td><td>3.3</td></tr><tr><td>Porosity</td><td>Φ</td><td></td></tr><tr><td>Latent heat source</td><td>Q<sub>lat</sub></td><td>W m<sup>-3</sup></td></tr><tr><td>Firn density</td><td>ρ<sub>f</sub></td><td>kg m<sup>-3</sup></td></tr><tr><td>Water density</td><td>ρ<sub>w</sub></td><td>1000 kg m<sup>-3</sup></td></tr><tr><td>Water source</td><td>R</td><td>s<sup>-1</sup></td></tr><tr><td>Water saturation</td><td>S</td><td></td></tr><tr><td>Effective water saturation</td><td>S*</td><td></td></tr><tr><td>Residual saturation</td><td>S<sub>r</sub></td><td>0.005</td></tr><tr><td>Freezing calibration constant</td><td>τ<sub>f</sub></td><td>2.0 × 10<sup>-5</sup> s<sup>-1</sup></td></tr><tr><td>Firn temperature</td><td>T</td><td>K</td></tr><tr><td>Ice melting point</td><td>T<sub>0</sub></td><td>273.15 K</td></tr><tr><td>Surface temperature correction</td><td>T<sub>rad</sub></td><td>2.0 K</td></tr><tr><td>Firn advection velocity</td><td>v<sub>z</sub></td><td>m s<sup>-1</sup></td></tr><tr><td>Water content</td><td>w</td><td>kg m<sup>-3</sup></td></tr></tbody></table>

insignificant (approximately of $4.0 \times 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.

# 4 Results and discussion

## 4.1 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 \text{ 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.

## 4.2 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. The energy balance is only adjusted for turbulent fluxes (via $z_0$) given that the radiative fluxes are directly measured in the field and the main uncertainty
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Meteorological data recorded by the AWS in summer 2012 and surface temperature calculated using measured infra-red surface emissions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** (a) Measured (calculated from long-wave radiation emissions, red line) and modeled surface temperature during summer 2012 at the AWS with $\delta = 0.0$ m and $S_r = 0.0$ (blue line) and $\delta = 0.025$ m and $S_r = 0.005$ (black line). (b, c) Measured and modeled temperature at 11 and 24 cm depth using SEB model (d) Modeled cumulated (black line) and hourly (blue line) surface melting from energy balance compared to measured snow height (red line). (e) Firn temperature measurements as a function of depth and time during summer 2012.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 4. (a) Modeled cumulative energy input in firn from surface energy balance ($E_b$, blue line), measured and modeled firn thermal energy $E_{fim}$ (red and black lines, respectively). (b) Energy difference between the energy input calculated from the SEB ($E_b$) and firn thermal energy ($E_{fim}$) compared to modeled water content in firn using model 2 (integrated on the vertical axis).</figcaption></figure>

comes from turbulent fluxes. The values of radiation penetration characteristic length $\delta$ and water residual saturation $S_f$ are chosen to minimize RMSE and bias between modeled and measured surface and subsurface temperatures. We find $\delta = 0.025$ m that corresponds to the value reported by Brandt and Warren (1993) and $S_f = 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⁻³) and corrected according to local air temperature. Two main melting events (25 July to 4 August, and 16 to 26 August) can be identified

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.02 K, $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_{fim}$ is calculated from

$$E_{fim} = \int_{0}^{D} c_{p}(z) \rho_{f}(z) T(z) dz \quad (\text{in J m}^{-2}), \quad (21)$$

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
<figure><img src="image_8.png" /><figcaption>**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).</figcaption></figure>

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.

## 4.3 Subsurface temperature and water content (model 2)

### 4.3.1 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 $\tau_f$) and the constant temperature correction $T_{rad}$ of surface temperature accounting for heating due to solar radiation penetration. These parameters were adjusted to $\tau_f = 2.0 \times 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
<figure><img src="image_9.png" /><figcaption>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).</figcaption></figure>

prevents water from percolating deep enough to explain the observed firn temperatures. The freezing calibration constant $\tau_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 $\tau_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 $\tau_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 half-hourly 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.

### 4.3.2 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 \times 10^{-6} \text{ 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):

$$M = \begin{cases} (T_{\text{max}} - T_0)a, & T_{\text{max}} > T_0 \\ 0, & T_{\text{max}} \le T_0 \end{cases} \quad (\text{in m w.e. d}^{-1}), \quad (22)$$
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Comparison between daily measured air ($T_{\text{air}}$) and surface ($T_{\text{s}}$) temperature between 3 July 2012 and 13 June 2013 at Col du Dôme.</figcaption></figure>

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.15 K) and $T_{\text{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 \times 10^{-3}$ K m$^{-1}$. The altitudinal gradient varies significantly between summer and winter ($-5.6 \times 10^{-3}$ K m$^{-1}$ in winter and $-6.5 \times 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 \times 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 \times 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.

## 4.4 Spatial variability of melting and subsurface temperature

### 4.4.1 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 refreezing in the firn.

<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

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). 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
<figure><img src="image_12.png" /><figcaption>Fig. 9. 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 Herron and Langway (1980).</figcaption></figure>

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

$$M = \tag{23}$$

$$ \begin{cases} (T_{\max} - T_0)a_{\text{PSR}}(x, y), & T_{\max} > T_0 \\ 0, & T_{\max} \le T_0 \end{cases} \quad \text{(in m w.e. d}^{-1}) $$

where $a_{\text{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_{\text{PSR}}$ for different PSR values. We found a quadratic relation-

ship (Fig. 10):

$$ \begin{aligned} a_{\text{PSR}} = & 3.3 \times 10^{-8} \times \text{PSR}^2 - 8.23 \times 10^{-6} \\ & \times \text{PSR} + 5.62 \times 10^{-4}. \end{aligned} \tag{24} $$

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.

## 4.5 Application of the simplified model at multi-decadennial 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⁻². 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 \times 10^{-3} \text{ W m}^{-2}$ for sites 1 and 2 and set to $30 \times 10^{-3} \text{ W m}^{-2}$ for site 3 according to measured basal temperature gradients. The basal
<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

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⁻¹ m⁻¹, a heat capacity of 7.5 × 10² J kg⁻¹ K⁻¹ and a density of 2.8 × 10³ kg m⁻³. 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⁻³, 5.70 × 10⁻³ and 5.83 × 10⁻³ K m⁻¹ 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⁻²) and equal to 3.5 × 10⁻⁴, 3.5 × 10⁻⁴, and 2.0 × 10⁻⁴ m w.e. d⁻¹ K⁻¹. 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

<figure><img src="image_14.png" /><figcaption>**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⁻⁴ m w.e. d⁻¹ K⁻¹.</figcaption></figure>

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⁻⁴ mw. e. d⁻¹ K⁻¹ instead of 2.0 × 10⁻⁴ mw. e. d⁻¹ K⁻¹ (Fig. 11, blue dashed line).

## 5 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 °C in summer, cold
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 re-
freezing 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 re-
freezing 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 radia-
tion penetration on firn temperature. Indeed, using measured
surface temperature as a Dirichlet condition amounts to ne-
glecting radiation penetration, which leads to significant un-
derestimation 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.

*Acknowledgements.* This study was funded by the AQWA European program (212250). We are grateful to G. Picard for providing the SEB code and we thank Météo-France for providing Lyon-Bron air temperature data. We thank H. Blatter and R. H. Giesen for helpful comments and suggestions, and we are grateful to H. Harder, who revised the English of this manuscript.

Edited by: M. van den Broeke

<figure><img src="image_15.png" /></figure>

The publication of this article is
financed by CNRS-INSU.

References

Aschwanden, A. and Blatter, H.: Mathematical modeling and nu-
merical simulation of polythermal glaciers, J. Geophys. Res.,
114, F01027, doi:10.1029/2008JF001028, 2009.

Brandt, R. E. and Warren, S. G.: Solar heating rates and temperature profiles in Antarctic snow and ice, J. Glaciol., 39, 99–110, 1993.

Brock, B. W., Willis, I. C., and Sharp, M. J., Measurement and parameterization of aerodynamic roughness length variations at Haut Glacier d'Arolla, Switzerland, J. Glaciol., 52, 281-297, 2006.

Brun, E., Martin, E., Simon, V., Gendre, C., and Coléou, C.: An energy and mass model of snow cover suitable for operational avalanche forecasting, J. Glaciol., 35, 333-342, 1989.

Brun, E., David, P., Sudul, M., and Brunot, G.: A numerical model to simulate snowcover stratigraphy for operational avalanche forecasting, J. Glaciol., 38, 13–22, 1992.

Calonne, N., Flin, F., Morin, S., Lesaffre, B., Rolland du Roscoat, S., and Geindreau, C.: Numerical and experimental investigations of the effective thermal conductivity of snow, Geophys. Res. Lett., 38 L23501, doi:10.1029/2011GL049234, 2011.

Colbeck, S. C.: Snow-crystal growth with varying surface tempera-
tures and radiation penetration, J. Glaciology, 35, 23–29, 1989.

Colbeck, S. C. and Davidson, G.: Water percolation through ho-
mogenous snow, The Role of Snow and Ice in Hydrology, Proc.
Banff. Symp., 1, 242-257, 1973.

Essery, R. L. H. and Etchevers, P.: Parameter sensitivity in simulations of snowmelt, J. Geophys. Res., 109, D20111, doi:10.1029/2004JD005036, 2004.
Faillettaz, J., Sornette, D., and Funk, M.: Numerical modeling of a gravity-driven instability of a cold hanging glacier: reanalysis of the 1895 break-off of Altelsgletscher, Switzerland, J. Glaciol., 57, 817–831, 2011.

Funk, M., Echelmeyer, K., and Iken, A.: Mechanisms of fast flow in Jakobshavns Isbrae, Greenland; Part II: Modeling of englacial temperatures, J. Glaciol., 40, 569–585, 1994.

Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a new-generation ice sheet model, Geosci. Model Dev., 6, 1299–1318, doi:10.5194/gmd-6-1299-2013, 2013.

Gilbert, A. and Vincent, C.: Atmospheric temperature changes over the 20th century at very high elevations in the European Alps from englacial temperatures, Geophys. Res. Lett., 40, 2102–2108, doi:10.1002/grl.50401, 2013.

Gilbert, A., Wagnon, P., Vincent, C., Ginot, P., and Funk, M.: Atmospheric warming at a high elevation tropical site revealed by englacial temperatures at Illimani, Bolivia (6340 m a.s.l., 16° S, 67° W), J. Geophys. Res., 115, D10109, doi:10.1029/2009JD012961, 2010.

Gilbert, A., Vincent, C., Wagnon, P., Thibert, E., and Rabatel, A.: The influence of snowcover thickness on the thermal regime of Tête Rousse Glacier (Mont Blanc range, 3200 m a.s.l.): consequences for water storage, outburst flood hazards and glacier response to global warming, J. Geophys. Res., 117, F04018, doi:10.1029/2011JF002258, 2012.

Haeberli, W., Alean, J.-C., Muller, P., and Funk, M.: Assessing risks from glacier hazards in high mountain regions: some experiences in the Swiss Alps, Ann. Glaciol., 13, 96–102, 1989.

Harrington, R. F., Bales, R. C., and Wagnon, P.: Variability of melt-water and solute fluxes from homogeneous melting snow at the laboratory scale, Hydrol. Process., 10, 945–953, 1996.

Herron, M. and Langway Jr., C.: Firn densification: An empirical model, J. Glaciol., 25, 373–385, 1980.

Hock, R.: A distributed temperature-index ice-and snowmelt model including potential direct solar radiation, J. Glaciol., 45, 101–111, 1999.

Hoffman, M. J., Fountain, A. G., and Liston, G. E.: Surface energy balance and melt thresholds over 11 years at Taylor Glacier, Antarctica, J. Geophys. Res., 113, F04014, doi:10.1029/2008JF001029, 2008.

Huggel, C., Haeberli, W., Kääb, A., Bieri, D., and Richardson, S.: An assessment procedure for glacial hazards in the Swiss Alps, Can. Geotech. J., 41, 1068–1083, 2004.

Hutter, K.: A mathematical model of polythermal glaciers and ice sheets, Geophys. Astro. Fluid, 21, 201–224, 1982.

Illangasekare, T. H., Walter, R. J., Meier, M. F., and Pfeffer, W. T.: Modeling of meltwater infiltration in subfreezing snow, Water Resour. Res., 26, 1001–1012, 1990.

Intergovernmental Panel on Climate Change (IPCC), Climate Change 2007: The Physical Science Basis, Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., Averyt, K. B., Tignor, M., and Miller, H. L., Cambridge Univ. Press, Cambridge, UK, 2007.

Louis, J.-F.: A parametric model of vertical eddy fluxes in the atmosphere, Bound. Lay. Meteorol., 17, 187–202, 1979.

Lüthi, M. and Funk, M.: Modeling heat flow in a cold, high-altitude glacier: Interpretation of measurements from Colle Gnifetti, Swiss Alps, J. Glaciol., 47, 314–323, doi:10.3189/172756501781832223, 2001.

Kuipers Munneke, P., van den Broeke, M. R., Reijmer, C. H., Helsen, M. M., Boot, W., Schneebeli, M., and Steffen, K.: The role of radiation penetration in the energy budget of the snowpack at Summit, Greenland, The Cryosphere, 3, 155–165, doi:10.5194/tc-3-155-2009, 2009.

Oke, T. R.: Boundary Layer Climates, 2nd Edn., 435 pp., Routledge, New York, 1987.

Paterson, W. S. B.: The Physics of Glaciers, 3rd Edn., 480 pp., Pergamon, New York, 1994.

Ryser, C., Lüthi, M., Bladow, N., Suckro, S., Funk, M., and Bauder, A.: Cold ice in the ablation zone: Its relation to glacier hydrology and ice water content, J. Geophys. Res.-Earth, 118, 693–705, doi:10.1029/2012JF002526, 2013.

Shimizu, H.: Air permeability of deposited snow, Contributions from the Institute of Low Temperature Science, A22, 1–32, 1970.

Skidmore, M. L. and Sharp, M. J.: Drainage system behaviour of a High-Arctic polythermal glacier, Ann. Glaciol., 28, 209–215, 1999.

Suter, S.: Cold firn and ice in the Monte Rosa and Mont Blanc areas: spatial occurrence, surface energy balance and climatic evidence, PhD thesis, ETH Zürich, 2002.

Vincent, C., Le Meur, E., Six, D., Possenti, P., Lefebvre, E., and Funk, M.: Climate warming revealed by englacial temperatures at Col du Dôme (4250 m, Mont Blanc area), Geophys. Res. Lett., 34, L16502, doi:10.1029/2007GL029933, 2007a.

Vincent, C., Le Meur, E., Six, D., Funk, M., Hoelzle, M., and Preunkert, S.: Very high-elevation Mont Blanc glaciated areas not affected by the 20th century climate change, J. Geophys. Res., 112, D09120, doi:10.1029/2006JD007407, 2007b.

Wagnon, P., Sicart, J.-E., Berthier, E., and Chazarin, J.-P.: Winter-time high-altitude surface energy balance of a Bolivian glacier, Illimani, 6340 m above sea level, J. Geophys. Res., 108, 4177, doi:10.1029/2002JD002088, 2003.

Wilson, N. J. and Flowers, G. E.: Environmental controls on the thermal structure of alpine glaciers, The Cryosphere, 7, 167–182, doi:10.5194/tc-7-167-2013, 2013.

Wright, A. P., Wadham, J. L., Siegert, M. J., Luckman, A., Kohler, J., and Nuttall, A. M.: Modeling the refreezing of meltwater as superimposed ice on a high Arctic glacier: A comparison of approaches, J. Geophys. Res., 112, F04016, doi:10.1029/2007JF000818, 2007.

Zwinger, T. and Moore, J. C.: Diagnostic and prognostic simulations with a full Stokes model accounting for superimposed ice of Midtre Lovénbreen, Svalbard, The Cryosphere, 3, 217–229, doi:10.5194/tc-3-217-2009, 2009.