Formation of glacier tables caused by differential ice melting: field observation and modelling

Glacier tables are structures frequently encountered on temperate glaciers. They consist of a rock supported by a narrow ice foot which forms through differential melting of the ice. In this article, we investigate their formation by following their dynamics on the Mer de Glace (the Alps, France). We report field measurements of four specific glacier tables over the course of several days, as well as snapshot measurements of a field of 80 tables performed on a given day. We develop a simple model accounting for the various mechanisms of the heat transfer on the glacier using local meteorological data, which displays a quantitative agreement with the field measurements. We show that the formation of glacier tables is controlled by the global heat flux received by the rocks, which causes the ice underneath to melt at a rate proportional to the one of the surrounding ice. Under large rocks the ice ablation rate is reduced compared to bare ice, leading to the formation of glacier tables. This thermal insulation effect is due to the warmer surface temperature of rocks compared to the ice, which affects the net long-wave and turbulent fluxes. While the short-wave radiation, which is the main source of heat, is slightly more absorbed by the rocks than the ice, it plays an indirect role in the insulation by inducing a thermal gradient across the rocks which warms them. Under a critical size, however, rocks can enhance ice melting and consequently sink into the ice surface. This happens when the insulation effect is too weak to compensate for a geometrical amplification effect: the external heat fluxes are received on a larger surface than the contact area with the ice. We identified the main parameters controlling the ability of a rock to form a glacier table: the rock thickness, its aspect ratio, and the ratio between the averaged turbulent and short-wave heat fluxes.

1 Introduction

A wide variety of spectacular shapes and patterns formed through differential ablation can be found in nature: surface patterns known as rillenkarren result from the dissolution of soluble rocks (Cohen et al., 2016, 2020; Claudin et al., 2017; Guérin et al., 2020), tafoni are cavernous rock domes dug by salt crystallization during wetting–drying cycles (Huinink et al., 2004), scallops and sharp pinnacles are created by convective flows (Huang et al., 2020; Weady et al., 2022), mushroom rocks undergo erosion of their base by strong particle-laden winds (Mashaal et al., 2020), and hoodoos consist of a hard stone protecting a narrow column of sedimentary rock from rain-induced erosion (Young and Young, 1992; Bruthans et al., 2014; Turkington and Paradise, 2005).

On ice and snow surfaces, similar structures can be found: slender snow blades known as penitentes (Mangold, 2011; Bergeron et al., 2006; Claudin et al., 2015), as well as blue ice ripples observed in Antarctica (Bintanja et al., 2001) and on Mars (Bordiec et al., 2020), are caused by sublimation; suncups are bowl-shaped depressions found at the surface of a snow patch (Rhodes et al., 1987; Betterton, 2001; Mitchell and Tiedje, 2010); scallops are regularly spaced patterns of surface indentations at the ice–water or ice–air interface (Bushuk et al., 2019); ice sails (or pyramids) are larger bare-ice structures that form on debris-covered glaciers (Evatt et al., 2017; Fowler and Mayer, 2017); Zen stones found on frozen lakes are pebbles resting on a delicate ice pedestal protected from sun-induced sublimation (Taberlet and Plihon, 2021); and glacier tables (see Fig. 1c) form when a foot of ice resists melting due mostly to thermal insulation provided by a large rock (Agassiz, 1840; Bouillette, 1933, 1934; Smiraglia and Diolaiuti, 2011; Hénot et al., 2021).

On a temperate glacier, the ice ablation rate is influenced by the presence of debris on its surface. Indeed, a dense debris layer covering an ice surface can, when thin enough (less than 0.5 cm), enhance the ice melting compared to a bare-ice surface or, if thick enough, act as an insulation layer and reduce the melting rate (Östrem, 1959). The insulation effect is well captured by complex energy balance models which use meteorological data as input parameters (Reid and Brock, 2010; Collier et al., 2014) and more recently by enhanced temperature index models (Carenzo et al., 2016; Moeller et al., 2016). The melt amplification effect for thin debris layers has been explained by its patchiness (Reid and Brock, 2010) or by its porosity to airflow (Evatt et al., 2015). At a more local scale, patches of dirt or ashes on a glacier are known to lead to the formation of ice structures known as dirt cones, which consist of ice cones covered with a centimetric layer of dirt (Swithinbank, 1950; Krenek, 1958; Drewry, 1972; Betterton, 2001).

Figure 1 (a–c) Formation of a natural granite glacier table (rock 1; see Table 1).

Glacier tables (see Fig. 1c) are structures frequently encountered on the levelled part of temperate glaciers (Agassiz, 1840). They typically consist of a metre-sized rock supported by a column (or foot) of ice. They form over the course of a few days to a few weeks at the end of the spring, from May to June, and may progressively disappear during the summer. Very large tables however, whose size can reach 10 m, can last for several years (Bouillette, 1933, 1934). Glacier tables form because the ablation rate of the ice is lower under the rock than at the air–ice interface. When the ice foot becomes too thin, the table falls, usually on its south side. The ice foot progressively disappears, and the rock can potentially form another table if it is not too late in the season. Similarly to what is known for debris layers, smaller rocks tend not to form tables but can instead increase the melting rate and gradually sink into the ice, creating narrow and deep holes which can reach up to 15 cm in depth (McIntyre, 1984). In a previous work (Hénot et al., 2021), artificial glacier tables were reproduced experimentally in a lab-controlled environment (constant temperature, humidity, absence of wind), at a centimetric scale. The study focused on the initial behaviour of pattern formation using cylindrical “rocks” of various sizes, the aspect ratio and materials, initially resting on a flat ice surface. Although this small-scale study under controlled conditions allowed one to understand the physical mechanisms that could play a role in glacier table formation, it did not encompass the complexity of the energy balance on a natural glacier, in particular the effects of the direct solar irradiation and of the wind.

In this article, we report field observations made on the Mer de Glace (French Alps) of the formation dynamics of glacier tables monitored over the course of a few days, as well as a systematic measurement of already-formed tables on a given day. The Mer de Glace is a temperate glacier, meaning that the ice temperature is always given by the melting temperature of water: T_ice=273 K. We use local meteorological data to fit the ice ablation rate and characterize the heat transfer mechanisms at the surface of the glacier. We then develop a 1D conduction model taking into account the effect of the solar irradiation as well as sensible and long-wave heat fluxes, which is in excellent agreement with the field measurements and which illustrates the synergistic effect between solar irradiation and sensible flux responsible for glacier table formation.

Figure 2 (a) Map of the Mer de Glace in the French Alps (Google Earth image © 2021 Maxar Technologies) locating the measurement site of 3 June 2021 (1), the time-lapse observation sites of June 2021 (2) and June 2019 (3), and the Requin AWS (4). (b) The 3D schematics of a typical glacier table defining the height H(t) of the ice foot and the average thickness h and widths d₁ and d₂ of the rock. (c) Picture of a portion of the glacier taken in June 2021 showing natural glacier tables, from which the data shown in Fig. 3 were obtained.

2 Observation

2.1 Location and definitions

In this article, we report two sets of observations made on the Mer de Glace: a time-resolved camera recording of the formation and evolution of four glacier tables and a field observation of 80 glacier tables.

The lower part of the Mer de Glace (below a 2100 m altitude) is largely covered with granitic debris, with sizes ranging from submillimetric up to several metres. In the following, their dimensions are characterized by their thickness h and their larger and smaller widths d₁ and d₂ (see Fig. 2b). We define an effective width of the rock $d_{\text{eff}}=\mathrm{2}{d}_{\mathrm{1}}{d}_{\mathrm{2}}/(d_{\mathrm{1}}+d_{\mathrm{2}})$, whose expression will be justified in Sect. 3.2, and the aspect ratio $\mathit{\beta }=h/d_{\text{eff}}$. The vertical distance from the bottom of a rock to the surface of the ice far from it is denoted H. This quantity is either positive if the rock forms a table, in which case it corresponds to the height of the ice foot, or negative if the rock sinks into the ice surface, in which case it corresponds to the depth of the hole.

2.2 Time evolution of four large tables

Time-lapse images were obtained using an autonomous solar-powered camera (Enlaps Tikee) placed on three 1.5 m long wood rods set into the ice. Pictures (4608 px×3456 px) were taken every 30 min between 05:00 and 22:00 (LT, UTC+2) over 5 to 7 d until the camera fell on the ice due to the melting around the supporting rods. The movements of the device were corrected by tracking two fixed points on the background of the images. The positions of the top of the rocks were then manually pointed onto each image. The formation of glacier tables created by four granite rocks of various shapes and dimensions was recorded in two time periods: A (7–14 June 2019) and B (4–9 June 2021), situated at locations 3 and 2 respectively on the map of Fig. 2a, as summarized in Table 1. Rocks 1, 2 and 4 were moved on clean, flat ice, in front of the camera, the day before the recording started in order to provide a controlled initial state in which the rocks are lying on a horizontal ice surface (H=0; see Fig. 1a). Rock 3 was already standing on an ice foot, approximately 1 m high, at the beginning of period B (see the Supplement). The evolution of the vertical position of these rocks is plotted in Fig. 4a and b. The position of the ice surface was followed using a scaling rod, embedded into the ice and located approximately 2 m away from the rocks.

The surface temperature of rock 4 was measured every 5 min during period B using thermocouples and a homemade battery-powered device (Arduino MKR ZERO and EVAL-CN0391-ARDZ Shield). The wind speed u_a, air temperature T_a, air specific humidity q_a and solar radiative flux Φ were measured at the Requin automatic weather station (AWS) (see Nadeau et al., 2009, who use identical devices, for a detailed description), at z_m=5 m above ground (located at 45.8846^∘ N, 06.9297^∘ E; see marker 4 in Fig. 2a), 600 m higher than and 3 km away from the measurement site. The time resolution was 1 h in 2019 and 15 min in 2021. During a 2-week time period in July 2021, we also measured the air temperature at location 2, 3 m above the glacier surface (see the Supplement), which was systematically 2.5 ^∘C higher than the AWS data. In the rest of the article, we use the AWS data, to which an offset was added: $T_{\mathrm{a}}=T_{\text{a AWS}}+\mathrm{2.5}$ ^∘C. We do not expect the measured solar flux and the wind speed to be significantly affected by the distance from the measurement site. The incoming long-wave radiation coming from the atmosphere Q_{LW atm ↓} was obtained from the S2M (SAFRAN–SURFEX/ISBA–Crocus–MEPRA) reanalysis, which combines information from numerical weather prediction models and in situ meteorological observation to estimate massif-averaged meteorological data with a 1 h time resolution (Vernay et al., 2022a) (the parameters used to compute the meteorological data in this study are the closest available that are representative of the real field, i.e. massif 3, 2100 m altitude, 20^∘ slope facing north). The meteorological data are displayed in Fig. 4c–l.

2.3 An 800 m² field comprising 80 glacier tables

On 3 June 2021, we systematically measured, in an area of 10 m×80 m located at point 1 in Fig. 2a, the dimensions of granite rocks (thickness and widths) as well as either the height of the ice foot supporting them (H>0) or the depth of penetration in the ice (H<0). The data of H are plotted in Fig. 3 as a function of the thickness h and the widths d₁ and d₂ of the rocks.

Table 1 Characteristics of the rocks studied.

Figure 3 Raw data of the observation made at location 1 (see Fig. 2a) on 3 June 2021: height of the ice foot (H>0) or of the penetration depth of the rock in the ice (H<0) as a function of the rock thickness h (a) and widths d₁ and d₂ (b). On the date of these observations the rocks were still standing on their ice feet and the maximum heights had presumably not been reached yet.

Figure 4 (a, b) Vertical position of the ice surface (blue markers) and of rocks 1 to 4 (coloured markers; see Table 1). The dashed black line corresponds to the model for ice ablation (see Sect. 3.1). The coloured solid lines correspond to the model for rocks (see Sect. 3.2). These models use meteorological data measured at the Requin AWS as input: wind speed u_a (c, d), air temperature T_a (e, f), air specific humidity q_a (g, h) and solar radiation Φ (k, l) as well as the long-wave radiation coming from the atmosphere Q_{LW atm ↓} obtained from the S2M reanalysis (i, j). The dotted lines correspond to the ice surface values. The date format is year-month-day.

3 Model

3.1 Energy balance at the ice surface

In the following, we characterize the heat flux balance on the glacier surface by linking the ablated ice thickness to local meteorological data in the same way as it has previously been carried out in the literature (Hock, 2005; Conway and Cullen, 2013; Fitzpatrick et al., 2017). The net incident heat flux causing the glacier surface to melt with an open environment can be expressed as

$$\begin{array}{}\text{(1)}& Q_{\text{open env.}\to \text{ice}}=Q_{\text{SW}}^{\mathrm{i}}+Q_{\text{LW}}^{\mathrm{i}}+Q_{\mathrm{H}}^{\mathrm{i}}+Q_{\mathrm{E}}^{\mathrm{i}}+Q_{\mathrm{R}}^{\mathrm{i}},\end{array}$$

where $Q_{\text{SW}}^{\mathrm{i}}$ and $Q_{\text{LW}}^{\mathrm{i}}$ are the net short-wave and long-wave radiative fluxes, $Q_{\mathrm{H}}^{\mathrm{i}}$ and $Q_{\mathrm{E}}^{\mathrm{i}}$ are the turbulent heat fluxes corresponding to sensible and latent heat respectively, and $Q_{\mathrm{R}}^{\mathrm{i}}$ is the flux associated with rain. The short-wave radiation (when neglecting the reflected radiation by the surrounding terrain) is given by $Q_{\text{SW}}^{\mathrm{i}}=\mathrm{\Phi }(\mathrm{1}-{\mathit{\alpha }}_{\text{ice}})$, where α_ice is the ice albedo and Φ is the incoming solar radiative flux (direct and diffuse) given by the meteorological S2M reanalysis model. The net long-wave flux comes from a balance between the radiation received by the ice from the atmosphere $Q_{\text{LW atm}\downarrow }^{\mathrm{i}}$ and the radiation emitted by the ice $Q_{\text{LW}\uparrow }^{\mathrm{i}}$ at temperature T_ice. Assuming an emissivity equal to unity 1, $Q_{\text{LW}\uparrow }^{\mathrm{i}}=\mathit{\sigma }{T}_{\text{ice}}^{\mathrm{4}}$, where $\mathit{\sigma }=\mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{4}}$, is the Stefan–Boltzmann constant. This leads to $Q_{\text{LW}}^{\mathrm{i}}=Q_{\text{LW atm}\downarrow }-\mathit{\sigma }{T}_{\text{ice}}^{\mathrm{4}}$. The long wave flux received from the atmosphere $Q_{\text{LW atm}\downarrow }^{\mathrm{i}}$ is given by the meteorological S2M reanalysis model. The turbulent fluxes transmitted to a surface are assumed to be proportional to the wind velocity u_a:

$$\begin{array}{}\text{(2)}& {\displaystyle }{Q}_{\mathrm{H}}^{\mathrm{i}}={\mathit{\rho }}_{\mathrm{a}}{c}_p{C}_{\mathrm{H}}{u}_{\mathrm{a}}(T_{\mathrm{a},z_{\mathrm{m}}}-T_{\text{surface}}),\text{(3)}& {\displaystyle }{Q}_{\mathrm{E}}^{\mathrm{i}}={\mathit{\rho }}_{\mathrm{a}}{\mathcal{L}}_{\mathrm{v}}{C}_{\mathrm{E}}{u}_{\mathrm{a}}(q_{\mathrm{a},z_{\mathrm{m}}}-q_{\text{surface}}),\end{array}$$

where ρ_a=0.98 kg m⁻³ is the air density at atmospheric pressure at a 2100 m altitude p_atm=785 hPa, c_p=1004 $\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{kg}}^{-\mathrm{1}}$ is the specific heat capacity of air, ℒ_v=2.48 MJ kg⁻¹ is the latent heat of evaporation for water, and $T_{\mathrm{a},z_{\mathrm{m}}}$ and $q_{\mathrm{a},z_{\mathrm{m}}}$ are the air temperature and specific humidity respectively at the height z_m above the surface. Assuming saturation near the melting ice surface, the specific humidity is q_surface=q_ice and is given by $\mathrm{0.622}{p}_{\text{vap}}\left(T_{\text{ice}}\right)/p_{\text{atm}}$, where p_vap(T_ice)=611 Pa is the vapour pressure of water at temperature T_ice (here and in the following, p_vap(T) values are obtained following the ITS-90 formula for water; Hardy, 1998). The exchange coefficients C_H and C_E are assumed to be identical and are given by the simplest form found classically in the literature assuming a neutral atmospheric stability (Conway and Cullen, 2013):

$$\begin{array}{}\text{(4)}& C_{\mathrm{H}}=C_{\mathrm{E}}={\displaystyle \frac{k^{\mathrm{2}}}{\left(\ln \right(z_{\mathrm{m}}/z_{\mathrm{0}}){)}^{\mathrm{2}}}},\end{array}$$

where k=0.41 is the von Kármán constant and z₀ is a roughness length characterizing the surface (assumed to be equal for the momentum, temperature and vapour pressure). This length is denoted as z_0 ice for the ice surface. There was no precipitation and thus no associated heat flux on the glacier during the duration of our study ($Q_{\mathrm{R}}^{\mathrm{i}}=\mathrm{0}$). As the incoming heat flux directly causes the ice to melt, the vertical position of the ice surface z_ice(t) and the total ablated thickness of ice since t=0, H_ice(t), are given by

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{z}_{\text{ice}}\left(t\right)-z_{\text{ice}}\left(\mathrm{0}\right)& =-H_{\text{ice}}\left(t\right)=\\ & -{\mathcal{L}}_{\text{fus}}\underset{\mathrm{0}}{\overset{t}{\int }}{Q}_{\text{open env.}\to \text{ice}}\left(t\right)\mathrm{d}t,\end{array}\end{array}$$

where ℒ_fus=303 MJ m⁻³ is the volumetric enthalpy of fusion for ice. This model relies on meteorological data and has two adjustable parameters: the ice surface roughness z_0 ice and its albedo α_ice. These adjustable parameters will be determined from the model fit to the ice ablation rate measured away from the glacier tables in Sect. 4.1.

3.2 Glacier table formation model

In Hénot et al. (2021), an analytical model was developed in order to explain the formation of artificial glacier tables made of cylindrical caps in a controlled environment, in which the heat transfer (natural convection of air and infrared radiation from the enclosure) was modelled through an effective heat exchange coefficient. In the following, we adapt this model to the non-symmetrical geometry of the rocks and we take into account the more complex energy balance of the glacier. Here we only attempt to describe the vertical motion of the rock with respect to the ice surface and we do not consider the lateral melt of the ice foot in this model, but this is briefly discussed in Sect. 5.2.

For the sake of simplicity, the following main assumptions are made in the model: (1) the rocks are considered cuboids (see Fig. 2b). The area A_base of contact with the ice is considered constant during the table formation, although it evidently varies over time. Yet, as discussed in the Supplement, this has little effect on the total heat transfer. The top and bottom areas are thus $A_{\text{top}}=A_{\text{bottom}}=d_{\mathrm{1}}{d}_{\mathrm{2}}$, and the side area is $A_{\text{side}}=\mathrm{2}(d_{\mathrm{1}}+d_{\mathrm{2}})h$. (2) The surface temperature of the rock is assumed to be homogeneous and is denoted as T_rock, allowing the development of a 1D thermal conduction model. (3) The thermal process is assumed to be quasi-static, meaning that all transient thermal effects are neglected. The validity of this assumption is discussed in the Supplement.

The energy balance of a rock, taking into account its 3D structure, is summarized in the schematic of Fig. 5a. The short-wave radiative flux is received by the rock on a surface 〈A_sun〉. This area should depend on the shape and orientation of the rock with respect to the course of the sun, but for the sake of simplicity in the following it is simply assumed to be equal to the rock base surface area 〈A_sun〉=A_base (the validity of this assumption is discussed in the Supplement). The long-wave flux coming from the atmosphere Q_{LW atm ↓} is also assumed to be received only on the top surface of area A_base, but the rock also receives a long-wave flux on its side denoted as Q_{LW env ↓}. Strictly speaking, this quantity depends on the emissivity, surface temperature and view factor of all surfaces (ice, other rocks, terrain) seen from the rock sides (Lienhard, 2019) and thus cannot be simply estimated. In the following it will be kept as an adjustable parameter. The rock also emits a long-wave flux $Q_{\text{LW}\uparrow }^{\mathrm{r}}$ from its external surface (top and side). The sensible and latent fluxes $Q_{\mathrm{H}}^{\mathrm{r}}+Q_{\mathrm{E}}^{\mathrm{r}}$ are received from the air on the external surface of the rock (of area A_base+A_side). Finally a flux Q_rock→ice is transferred to the ice through the contact area A_base.

The quasi-static assumption implies that the flux balance is verified at each time t:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{A}_{\text{base}}{Q}_{\text{rock}\to \text{ice}}=& A_{\text{top}}{Q}_{\text{SW}}^{\mathrm{r}}\\ & +A_{\text{base}}(Q_{\text{LW atm}\downarrow }-Q_{\text{LW}\uparrow }^{\mathrm{r}})\\ & +A_{\text{side}}(Q_{\text{LW env}\downarrow }-Q_{\text{LW}\uparrow }^{\mathrm{r}})\\ & +(A_{\text{base}}+A_{\text{side}})(Q_{\mathrm{H}}+Q_{\mathrm{E}}).\end{array}\end{array}$$

The net short-wave flux is $Q_{\text{SW}}^{\mathrm{r}}=(\mathrm{1}-{\mathit{\alpha }}_{\text{rock}})\mathrm{\Phi }$ where α_rock=0.18 is the albedo of the granite rock (measured by Watson, 1971; see the Supplement). The granite emissivity is taken as equal to 1 (Michalski et al., 2004), leading to $Q_{\text{LW}\uparrow }^{\mathrm{r}}=\mathit{\sigma }{T}_{\text{rock}}^{\mathrm{4}}$. The turbulent fluxes are computed according to an equation similar to Eq. (3), for which the surface terms are those relative to the rock; i.e. T_surface=T_rock, $q_{\text{surface}}=\mathrm{0.622}{p}_{\text{vap}}\left(T_{\text{rock}}\right)/p_{\text{atm}}$ and z₀=z_0 rock. If Q_E>0, water from the air condenses on the rock surface, which justifies the value of q_surface (corresponding to 100 % humidity at the temperature T_rock). However, Q_E<0 would correspond to evaporation and would be limited by the quantity of water present on the rock surface. In the following Q_E is thus forced to 0 if Q_E<0.

The heat flux transmitted from the rock to the ice underneath can be estimated using a 1D conduction model as

$$\begin{array}{}\text{(7)}& Q_{\text{rock}\to \text{ice}}=\mathit{\lambda }{\displaystyle \frac{(T_{\text{rock}}-T_{\text{ice}})}{d_{\text{1D}}}},\end{array}$$

where λ is the thermal conductivity of granite (measured in a previous work – Hénot et al., 2021 – as λ=2.8 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, which is in agreement with the literature; Cho et al., 2009) and d_1D is a length scale for the conduction process in the rock. It was postulated and verified experimentally (see Fig. 5 of Hénot et al., 2021) that this length could be estimated by $d_{\text{1D}}=\mathit{\eta }{V}_{\text{rock}}/A_{\text{ext}}=\mathit{\eta }h/(\mathrm{1}+\mathrm{4}\mathit{\beta })$, where V_rock is the volume of the rock and η=2.5 is a numerical prefactor adjusted from the lab experiments. This expression justifies the definition of d_eff.

From Eqs. (6) and (7) and using the measured meteorological data (Φ(t), Q_{LW atm ↓}(t), u_a(t), T_a(t) and q_a(t)) as input, the rock surface temperature T_rock(t) can be computed at each time t by solving numerically non-linear Eq. (6). The flux Q_rock→ice(t) reaching the ice can then be computed. When positive, this leads to ice melting under the rock, lowering the altitude of the bottom of the rock z_rock:

$$\begin{array}{}\text{(8)}& z_{\text{rock}}\left(t\right)-z_{\text{ice}}\left(\mathrm{0}\right)=-{\mathcal{L}}_{\text{fus}}\underset{\mathrm{0}}{\overset{t}{\int }}{Q}_{\text{rock}\to \text{ice}}\left(t\right)\mathrm{d}t.\end{array}$$

The height of the ice foot supporting the rock is $H\left(t\right)=z_{\text{rock}}\left(t\right)-z_{\text{ice}}\left(t\right)$. The time origin corresponds to the rock lying at the surface of the ice ($H(t=\mathrm{0})=\mathrm{0}$). Note that at this stage, there are two adjustable parameters in the table formation model: Q_{LW env ↓} and z_0 rock (α_rock and η are taken from the literature). However, the roughness size z₀ in the turbulent coefficient is not clearly defined (Hock, 2005), but the resulting model is not very sensitive to its value. Thus in the following, z_0 rock will be taken as equal to that of the ice z_0 ice.

4 Results

4.1 Ice melting

The dashed black lines in Fig. 4a and b correspond to the model described in Sect. 3.1 with parameters α_ice=0.30 and z₀=0.34 mm adjusted on the concatenated data of periods A and B (see the Supplement for a discussion on the sensitivity). The agreement between the model and the field measurements is good: the error on z_ice(t) stays below 3 cm (except at the end of period A when the movement of the camera increases the uncertainty in the measurements). The value of the adjusted ice albedo corresponds to what is reported in the literature for aged surface glacier ice, at the beginning of summer (Brock et al., 2000). The value of the roughness length z₀ also falls in the range previously reported (10⁻⁴–10⁻² m) (Brock et al., 2006; Conway and Cullen, 2013). Figure 5 shows (in white) the time-averaged values of $Q_{\text{SW}}^{\mathrm{i}}$, $Q_{\text{LW}}^{\mathrm{i}}$ and $Q_{\mathrm{H}}^{\mathrm{i}}$, $Q_{\mathrm{E}}^{\mathrm{i}}$ as well as the total flux reaching the ice surface, Q_{open env→ice}, for both time periods. The short-wave flux varies weakly between periods A and B (from 150 to 200 W m⁻²), and the net long-wave flux is comparatively small (−20 to 3 W m⁻²). The total turbulent flux $Q_{\mathrm{H}}^{\mathrm{i}}+Q_{\mathrm{E}}^{\mathrm{i}}$ strongly differs between periods A and B. It decreases from 120 W m⁻² (which represents 50 % of the total flux) in period A to 25 W m⁻² (10 % of the total flux) in period B. This is due to the drastic difference in mean wind speed, from 6.5 m s⁻¹ during period A to 1.0 m s⁻¹ during period B, which explains the difference in the ice ablation rate (of the order of 8 cm d⁻¹ in period A and 5.5 cm d⁻¹ in period B).

Figure 5 (a) Schematics of the heat fluxes considered in the models for the ice surface melting (left) and glacier table formation. The colours denote the surfaces on which the fluxes are received. (b, c) Distribution of the heat fluxes averaged over the duration of observation as predicted by the time-resolved model for the ice surface (white) and for the four rocks studied (coloured).

4.2 Table formation

Rocks 1 and 2 formed glacier tables over the course of approximately a week (see Fig. 4a and b), after which they fell off the ice pedestal. Rock 3 was initially sitting on an ice foot, which grew bigger, while rock 4, although exhibiting dimensions close to those of rock 2 went down almost at the same rate as the bare-ice surface and did not form any significant glacier table (see the Supplement for pictures). The integration of the flux given by the model is shown in solid lines in Fig. 4a and b. The best agreement with the observed data was obtained for a value of the adjustable parameter $Q_{\text{LW env}\downarrow }=\mathrm{340}$ W m⁻². This corresponds to the long-wave flux emitted by a surface with emissivity 1 at a temperature 278.3 K, which seems reasonable as the physical meaning of this quantity is the mean temperature of the environment seen by the sides of the rocks (ice surface, others rocks and terrain surrounding the glacier). The agreement between the measurements and the model prediction is overall good for rocks 1 to 3. The model overestimates slightly the ability to form a table of rock 4, but the error (integrated over 6 d) stays under 3 cm for total ice ablation of 37 cm. Figure 6 shows the measured temperature of the top surface of rock 4 (solid line) as well as the value of T_rock(t) predicted by the model. The two curves qualitatively display the same behaviour. The mean temperature is slightly overestimated by the model on average (by 2 ^∘C), but the typical minimum (0–3 ^∘C) and maximum (22–24 ^∘C) daily temperature remain close.

Figure 6 Temperature T_rock of the top surface of rock 4 measured using a thermocouple (light green) and predicted by the model using the meteorological data of period B (dark line). The date format is year-month-day.

Figure 5 shows for each rock (in colour) the surface- and time-averaged radiative and turbulent flux computed from the model as well as the resulting heat flux transmitted to the ice underneath. The short-wave flux received on the top surface of the rocks is 17 % higher than the one received by the ice due to the difference in albedo. However the net long-wave and turbulent fluxes are significantly reduced due to the fact that on average the rocks' surface temperature 〈T_rock〉 are higher than T_ice. For rocks 1 to 4, the averaged surface temperatures predicted were 8.2, 6.3, 13.0 and 7.2 ^∘C respectively. The turbulent flux reduction is more important during period A than during period B due to the higher mean wind speed. Note that this effect is geometrically amplified (by a factor 1+4β) as these negative fluxes are integrated over the external surface of the rocks, which is larger than the base surface in contact with ice (this is especially true for rock 4, which has the largest aspect ratio). For all the rocks studied, the predicted flux causing ice to melt under the rocks was on average reduced compared to the flux received by the ice surface (from −15 % for rock 4 to −88 % for rock 3), leading to the formation of tables. Note that this amplification effect can also, depending of the shape and size of the rocks, have the opposite consequence of causing a rock to sink into the ice surface, as will be shown in Sect. 5.2.

Figure 7a shows, for each rock, the height of the ice foot, $H\left(t\right)=z_{\text{rock}}\left(t\right)-z_{\text{ice}}\left(t\right)$, as a function of the ablated ice thickness, $H_{\text{ice}}\left(t\right)=z_{\text{ice}}\left(t\right)-z_{\text{ice}}\left(\mathrm{0}\right)$. On average these quantities appear to be proportional. The slope corresponds to a dimensionless growth rate of the glacier table, denoted as $\langle H/H_{\text{ice}}\rangle$: a unity slope would correspond to no melt under the rock (the ice foot height would then correspond exactly to the ablated ice thickness); conversely a zero slope would correspond to a rock descending at the same rate as the ice surface around it and thus never forming a glacier table. The solid lines in Fig. 7a correspond to linear adjustments. The deduced slopes $\langle H/H_{\text{ice}}\rangle$ are plotted for each rock in Fig. 7b as a function of their thickness h (markers). The value predicted by the model is plotted in the same figure for periods A (dashed lines) and B (solid lines) for the values of β corresponding to the four rocks as a function of the thickness of the rock. Moving along a line thus corresponds to a change in scale but not in shape of the rocks. The dimensionless growth rates of the tables are accurately predicted by the model. It is visible that independently of the rock shape and time period, a small enough rock would not form a table and instead sink into the ice surface ($\langle H/H_{\text{ice}}\rangle <\mathrm{0}$).

Figure 7 (a) Evolution of the height H of the ice foot under rocks 1–4 as a function of the ablated ice thickness H_ice. The coloured lines correspond to linear adjustments of the data. The solid black line has a slope of 1. (b) The mean ratio of the height H of the ice foot under a rock and of the total ablated ice thickness H_ice since the beginning of the table formation, as a function of the thickness h of the rock. The markers correspond to the slopes of (a) for each rock. The lines correspond to the application of the model for the meteorological values measured in period A (dashed lines) and B (solid lines) and the aspect ratio β of each rock given in Table 1.

4.3 Analysis of the 80 glacier tables field

In Fig. 3, showing the height of the ice structure formed by 80 rocks 8 d after the melt of the snow layer, is it visible that a critical size exists (4–10 cm in thickness and 9–13 cm in width) above which rocks tend to form glacier tables. The larger the rock, the greater its ability to protect the ice underneath from melting and, assuming all tables started their formation at the same time, the higher the ice foot. Below this critical size, however, the rocks have a tendency to sink into the ice surface. Having successfully modelled the evolution of the four larger glacier tables, the model presented in Sect. 3.2 can be applied to the data obtained on a given day for a field comprising 80 rocks.

The snow layer covering that part of the glacier finished melting on 26 May, leaving the ice exposed. The meteorological data of the time period C (from 26 May until 3 June), alongside the previously adjusted parameters α_ice and z₀, allow us to compute the thickness of ablated ice $H_{\text{ice}}=\mathrm{0.45}\pm \mathrm{0.05}$ m during this period (see Fig. S5 in the Supplement). The observed values of H are plotted in Fig. 8a as a function of the ones predicted by the model for each rock, assuming that the structure formation started on 26 May (H_{26 May}=0). The rocks that have sunk into the ice surface (H<0) are shown in blue. For those, the prediction of the model is, as expected, not accurate. Indeed, a rock that has penetrated into the ice has a greater contact surface area with the ice and a smaller external surface open to the atmosphere than what is taken into account in the model. It is thus cooler, and its (negative) dimensionless growth rate is smaller (in absolute terms) than that given by the model, as observed. The rocks that have formed glacier tables (H>0) are shown in red. The threshold between rocks forming or not forming a table (the point (0,0) in the graph) is well predicted by the model. For the smaller structures (H<0.3 m), the height of the ice foot is also well predicted (with an error smaller than 5 cm). Note that this required no new adjustable parameter and this prediction is the result of a model integrated over 8 d in a time period different from A and B, during which the parameters of the model were adjusted. For ice feet taller than 0.3 m, however, the prediction of the model is not good. For some H>H_ice, and the assumption H_{26 May}=0 is thus obviously not valid. We infer that larger rocks emerged from the snow layer before 26 May and started to form a table earlier.

Figure 8 (a) Observation versus model for the distance H between the bottom of rocks and the ice surface for the data measured on 3 June 2021 (see Fig. 3). The data are represented in red if a table is observed (H>0) and in blue otherwise. The solid black line has a slope of 1, and the dotted lines correspond to the ablated ice thickness H_ice since the melting of the snow layer (see the Supplement). (b) Diagram of the ability of a rock to form a table as a function of its typical width d_eff and the ratio $\langle Q_{\mathrm{H}}^{\mathrm{i}}+Q_{\mathrm{E}}^{\mathrm{i}}\rangle /\langle Q_{\text{SW}}^{\mathrm{i}}\rangle$. The solid lines correspond, for three values of the aspect ratio β, to the delimitation predicted by the model with averaged meteorological data, between a rock forming a table and a rock sinking into the ice surface. The small circles correspond to the data of (a) with the same colour code. The larger circles correspond to rocks 1–4. (c) Distribution of the aspect ratio β of the rocks close to the transition between the two regimes ($\left|H\right|<\mathrm{0.2}$ m). The vertical dashed lines correspond to the limits plotted in (b) with the same colour code.

5 Discussion

5.1 Critical diameter

The critical diameter, corresponding to the transition between the two regimes described in the previous section (the rock sinks into the ice or forms a table), is denoted as d_crit. We chose to consider a critical width instead of a critical thickness as the transition seems sharper when considering this quantity in Fig. 3. This is reminiscent of the results obtained in a controlled environment in the absence of solar radiation, for which the transition occurred for a critical diameter of the cylinders, independently of their aspect ratio (Hénot et al., 2021).

From our model, given a set of meteorological data and an aspect ratio β, a critical width d_crit can be computed. In the following, unlike what has been done until now, typical time-averaged meteorological data (〈Φ〉=210 W m⁻², $\langle Q_{\text{LW atm}\downarrow }\rangle =\mathrm{315}$ W m⁻², 〈T_a〉=7 ^∘C, 〈q_a〉=0.0058 kg kg⁻¹ and 〈u_a〉 ranging from 0 to 20 m s⁻¹) will be used as inputs (we have checked that this pre-averaging of the input data affects the output of the model by less than a few percent). Figure 8b shows a diagram of the ability of a rock to form a table as a function of its effective width d_eff (which in practice can be replaced by the mean width) and the ratio $\langle Q_{\mathrm{H}}^{\mathrm{i}}+Q_{\mathrm{E}}^{\mathrm{i}}\rangle /\langle Q_{\text{SW}}^{\mathrm{i}}\rangle$, which is a quantity proportional to 〈u_a〉 and whose relevance will be discussed in the next paragraph. The four monitored rocks as well as the field data from the 80 rocks found in the 800 m² field are plotted in this diagram. The distribution of aspect ratio β of rocks (46 of them) close to the transition (defined by $\left|H\right|<\mathrm{0.2}$ m) is plotted in Fig. 8c. The bounding values (β=0.16 and β=0.66) as well as the value corresponding to the maximum of the distribution (β=0.27) are used to plot, in Fig. 8b, the values of d_crit predicted by the model. This threshold, which depends slightly on β, shows a good agreement with the transition visible in the field data. This diagram illustrates the sensitivity of the glacier table formation to the meteorological conditions and especially the mean wind speed: d_eff is divided by 2 between time periods A and B.

5.2 Physical discussion

A rock forms a glacier table when the heat flux reaching the ice underneath is reduced compared to that received by the bare-ice surface around it. If, on the contrary, this heat flux is amplified, the rock will sink into the ice surface. Depending on the rock size, both phenomena are observed on temperate glaciers. During summer, the main source of heat flux causing the ice to melt is the short-wave radiation coming from the sun. Due to the lower albedo of the rock, this flux is slightly amplified compared to what is received by a similar area of bare ice and thus does not directly favour the formation of glacier tables. However this induces a heat flux across the thickness of the rock that elevates its surface temperature T_rock (the bottom surface staying at T_ice) and strongly affects the net infrared radiation and turbulent fluxes. As can be seen in Fig. 5, the sign of the sensible flux can even change (if $T_{\text{rock}}>T_{\mathrm{a}}>T_{\text{ice}}$, the wind warms the ice surface but cools the rock) and induces a negative net short-wave flux (the rock flux lost being proportional to $T_{\text{rock}}^{\mathrm{4}}$). At the end, if the rock is large enough (so that its surface temperature can rise enough), all these combined effects lead to the formation of a glacier table. As illustrated by the difference between periods A and B, the wind intensity can vary significantly at the surface of the glacier. The sensible flux being one motor of glacier table formation (with emitted long-wave radiation), this affects the critical size d_crit. The ratio $\langle Q_{\mathrm{H}}^{\mathrm{i}}+Q_{\mathrm{E}}^{\mathrm{i}}\rangle /\langle Q_{\text{SW}}^{\mathrm{i}}\rangle$ used on the abscissa in Fig. 8b demonstrates this effect. When this ratio is small (<0.1), the structure formation is controlled by the long-wave emission of the rock, while a large value (>1) would correspond to a regime dominated by the effect of sensible flux. Note that at this point, a 1D model with d_1D=h and no geometrical amplification would have achieved the same result. However, the 3D structure has its importance: the heat transfer with the environment (short-wave and turbulent) that cools the rock is exchanged on a surface area with is much larger than the base area of a rock. In our model we considered the total external surface of the rock, which is larger by a factor of 1+4β typically of the order of 2–3 than A_base. As detailed in the Supplement, a 1D model would, in order to verify Eq. (6), significantly overestimate T_rock (by 150 ^∘C during sunlight for rock 3, which is not a physically acceptable value). This 3D effect also explains the influence of the rock shape on its ability to form a table.

The formation of glacier tables in the natural environment of a temperate glacier results from a more complex energy balance than the idealized lab-controlled study conducted previously (Hénot et al., 2021) while the main physical ingredient stays the same: the greater temperature of the rock surface, compared to the ice, reduces the incoming heat flux and thus the melt rate under the rock. In the absence of solar radiation, the threshold controlling the ability of a rock to form a table was determined by a trade-off between this effect and a geometrical amplification effect: a rock can act as a fin in amplifying the incoming heat flux, leading to a higher melt rate under the rock than for bare ice. While this effect also plays a role in natural conditions, the threshold is also significantly affected by the direct solar radiation that has two opposite effects: on the one hand, due to lower albedo, a higher flux is received by a rock than a bare-ice surface, but, on the other hand, the temperature gradient induced by the solar flux across a rock elevates its temperature, enhancing the cooling effect due to turbulent and long-wave flux. The first effect dominates for small rocks, while the latter dominates for larger rocks.

While the heat transfer across the rock controls the formation dynamics, the end of life of glacier tables and in particular the maximum height reached by ice feet will likely result from lateral melting of the ice column. This process is expected to be affected by the rock shape and size that would prevent radiative melting due to shading effects, leaving only turbulent flux and long-wave balance. Note also that when the ice pedestal becomes very slender, a heavy rock could also cause the ice column to creep.

6 Conclusions

We studied the formation of four glacier tables over the course of a week, and we measured the characteristics of 80 glacier tables on the Mer de Glace (Alps). We developed a simple model taking into account the infrared and solar radiation and turbulent heat flux received by a rock and by the glacier ice surface using local meteorological data, allowing us to quantitatively reproduce the glacier table formation dynamics and to identify the physical effects at its origin. The table formation is controlled by the ice melting under the rocks, and the ice foot growth rate is proportional to the ice ablation rate at the glacier surface. The size of the rocks is a determinant factor governing table formation: the bigger the rock, the higher and faster the ice foot supporting it will grow. Under a critical size, rocks show the opposite behaviour of sinking into the ice surface. The ability of a rock to form a table is controlled by the balance between two opposing effects: a thermal insulation effect, which depends strongly on the rock size, and a geometrical amplification effect linked to the fact that the external surface on which the rock receives an external heat flux is larger than its contact area with the ice. This second effect becomes dominant only for rocks smaller than a few tens of centimetres. The insulation effect originates from the warmer temperature of the rock surface compared to the ice, which reduces, even sometimes changing the sign of, the net long-wave and turbulent flux, ultimately reducing the heat available for ice melting under the rock. While a rock receives a slightly larger net short-wave flux (the main source of heat at the glacier surface) than ice because of the difference in albedo, this effect is too weak to compensate for the insulation effect for large rocks. The solar radiation, by inducing a strong thermal gradient across the thickness of the rock, raises the rock surface temperature, which also contributes to the insulation effect. To summarize, we identified three main parameters controlling the dimensionless growth rate of the ice foot $\langle H/H_{\text{ice}}\rangle$ that characterize the ability of a rock to form a table: the rock thickness h; its aspect ratio β; and the ratio $\langle Q_{\mathrm{H}}^{\mathrm{i}}+Q_{\mathrm{E}}^{\mathrm{i}}\rangle /\langle Q_{\text{SW}}^{\mathrm{i}}\rangle$, which depends on the meteorological conditions. When this ratio is higher than ≈0.1 (i.e. in windy periods), the turbulent heat exchange lowers the critical diameter d_crit (by up to 50 %), which facilitates glacier table formation. We did not consider in this study the lateral melting of the ice pedestal and the end of life of the structures, but this could be the subject of future work.