Surface heat fluxes at coarse blocky Murtèl rock glacier (Engadine, eastern Swiss Alps)

Abstract

We estimate the surface energy balance (SEB) of the Murtèl rock glacier, a seasonally snow-covered permafrost landform with a ventilated coarse blocky active layer (AL) located in the eastern Swiss Alps. We focus on the parameterisation of the turbulent heat fluxes. Seasonally contrasting atmospheric conditions occur in the Murtèl cirque, with downslope katabatic jets in winter and a strongly unstable atmosphere over the heated blocky surface in summer. We use a novel comprehensive sensor array both above the ground surface and in the coarse blocky AL to track the rapid coupling by convective heat and moisture fluxes between the atmosphere, the snow cover, and the AL for the time period September 2020–September 2022. The in situ sensor array includes a sonic anemometer for eddy-covariance flux above-ground and sub-surface long-wave radiation measurements in a natural cavity between the AL blocks. During the thaw seasons, the measurements suggest an efficient ($\sim$ 90 %) export of the available net radiation by sensible and latent turbulent fluxes, thereby strongly limiting the heat available for melting ground ice. Turbulent export of heat and moisture drawn from the permeable AL contributes to the well-known insulating effect of the coarse blocky AL and partly explains the climate resiliency of rock glaciers. This self-cooling capacity is counteracted by an early snow melt-out date, exposing the low-albedo blocky surface to the intense June–July insolation and causing reduced evaporative cooling due to exacerbated moisture scarcity in the near-surface AL during dry spells. With climate change, earlier snowmelt and increased frequency, duration, and intensity of heat waves and droughts are projected. Regarding the parameterisation of the turbulent fluxes, we estimated the year-round turbulent fluxes using a modified scheme. The monthly SEB is closed within 20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ except during the snowmelt months and under katabatic drainage winds in winter. Detected sensible turbulent fluxes from nocturnal ventilation processes, although a potentially important ground cooling mechanism, are within our 20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ uncertainty because nighttime wind speeds are low. Wintertime katabatic wind speeds needed to be scaled to close the SEB, which hints at the limits of parameterisations based on the Monin–Obukhov similarity theory in complex mountain terrain and katabatic drainage winds. The present work contributes to the process understanding of the SEB and climate sensitivity of coarse blocky landforms.

1 Introduction

Coarse blocky landforms such as rock glaciers, block fields, and talus slopes are covered by a thick clast-supported debris mantle typically $\sim$ 1–5 m thick. These “cold rocky landforms” exhibit lower ground temperatures compared with adjacent fine-grained or bedrock areas, a phenomenon referred to as “undercooling” or “algific” . This special ground thermal regime is owed to the interactive effects of energy exchange processes between the atmosphere and the ground arising from the debris mantle properties . Large blocks and typically sparse fine materials near the surface create a vast, highly connected and thus permeable pore space. Depending on the stability of the air column in the debris mantle controlled by temperature gradients, a variable part of the near-sub-surface debris mantle is ventilated and effectively participates in the turbulent heat exchange with the atmosphere. Heat and moisture stored in the ventilated near-surface sub-layer of the debris mantle are rapidly mobilised and contribute to turbulent fluxes at short ($\le$ hourly) timescales. These non-conductive sub-surface heat transfer processes related to storage and phase changes in heat, vapour, water, and ice produce the peculiar micro-climate observed in debris mantles. As a consequence, to calculate the surface energy balance, no unambiguous, clearly defined and fixed interface of energy conversion (“active surface”, ) separating radiative–convective processes in the atmosphere from (dominantly) diffusive processes in the sub-surface is available . These active surfaces do not necessarily coincide for different processes (radiation conversion, wind-forced and buoyancy-driven air convection, interception of precipitation, water flow) or quantities (solar radiation, momentum, water), and therefore, they vary in time and are not necessarily at ground surface. On seasonally snow-covered sites, snow controls ground–atmosphere heat exchange processes at the surface by its effects on albedo and thermal insulation .

These complex interactions between debris mantle, snow cover, and atmosphere demand continuous high-resolution sub-surface measurements in addition to the “classical” above-surface weather station measurements and ground temperature. However, gathering the necessary measurements of sub-surface parameters is challenging in remote mountain terrain, and only a few comprehensive data sets beyond ground temperatures exist in mountain permafrost, exceptions being and . They deployed a heat flux plate, ultrasound probes, conductometer, vapour traps, and reflectometer probes to characterise the ground hydro-thermal regime of a steep, permafrost-underlain scree slope. Instead, non-conductive heat transfer processes in block fields, talus slopes, and debris cover of glaciers have been inferred from their effect on the ground thermal regime by means of (high-resolution) temperature measurements , in cases aided by gas/smoke tracer experiments or thermal infrared imaging .

These challenges have been described in energy balance studies conducted on the Murtèl rock glacier, situated in a cirque in the Upper Engadine (eastern Swiss Alps). In this permafrost landform, the debris mantle overlies a perennially frozen rock glacier core, is seasonally frozen, and roughly coincides with the thermally defined active layer (AL). We will use the term “coarse blocky AL” throughout the text to point to both material property and thermal state. In a micro-climatological study by , significant deviations from a closed surface energy balance (SEB) of up to 78 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in winter and $-\mathrm{130}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in summer were found. These deviations exceed methodological uncertainties. The authors suggested that the apparent heat sink in summer and heat source in winter arose from processes unaccounted for in their calculations due to the lack of measurements, namely advective and convective heat transport in the coarse blocky AL. addressed the seasonal SEB imbalance by adopting a volumetric energy balance approach consisting of adding a porous interfacial buffer layer able to store and release heat and including radiative and sensible turbulent heat transfer in the coarse blocky AL. The integration of additional sub-surface heat transfer mechanisms and storage components reduced the deviations to 26 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in summer and $-\mathrm{29}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in winter. The remaining seasonal deviations were interpreted as latent heat effects of freezing and thawing of ice in the AL and at the permafrost table. This agreed well with long-term melt rates derived from photogrammetric and three-dimensional borehole deformation leading to subsidence estimates ($\sim$ 5 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$) .

In this work, we estimate the SEB of the Murtèl rock glacier using in situ measurements on and within its ventilated coarse blocky AL. We revisit the micro-climatological studies by and adding new measurements from a novel sensor array. Since the deviations present in these works were largely caused by uncertainties in the turbulent fluxes, here we focus on the parameter and parameterisations required for their calculation. We address two questions: a general one and a technical one. The general question is how large the individual surface heat fluxes on the Murtèl rock glacier are and how these are seasonally distributed. From a quantitative process understanding, we gain insight into how the insulating effect of the coarse blocky AL works and how the rock glacier, a thermally conditioned permafrost landform, responds to climate change. The technical question is which turbulent flux parameterisations and input parameters are appropriate for the study of a seasonally snow-covered ventilated landform in complex mountain terrain. Difficulties in calculating the turbulent fluxes on Murtèl arise at two spatial scales: the meso-scale relief (${\mathrm{10}}^{\mathrm{2}}$–${\mathrm{10}}^{\mathrm{3}}$ m) and sub-landform scale (${\mathrm{10}}^{-\mathrm{1}}$–${\mathrm{10}}^{\mathrm{1}}$ m). The meso-scale relief (${\mathrm{10}}^{\mathrm{2}}$–${\mathrm{10}}^{\mathrm{3}}$ m) modifies wind speeds and atmospheric stability in the Murtèl cirque, which must be accurately reflected by the flux parameterisations. In winter, persistent katabatic winds develop on the snow-covered rock faces and slopes and converge in the cirque. These katabatic jets determine the near-surface wind velocities and the vertical turbulent exchange. In summer, surrounding rock spurs weaken the regional valley wind in the sheltered Murtèl cirque. Thermals rising from the strongly heated debris surface create an unstable atmosphere with comparatively low horizontal wind speeds. We test different stability corrections commonly used for debris-covered glaciers. At sub-landform scale (${\mathrm{10}}^{-\mathrm{1}}$–${\mathrm{10}}^{\mathrm{1}}$ m), the coarse blocky AL is ventilated and participates in the convective exchange of momentum, heat, and moisture with the atmosphere, unless the ground is covered by a sufficiently thick snow cover. Heat and moisture can be drawn from an interfacial buffer layer coupled with the atmosphere. The active surface does not necessarily coincide with the ground or snow surface. Where should we appropriately measure the meteorological variables, namely surface temperature and surface humidity, needed to estimate the turbulent fluxes – on the ground/snow surface or at some depth? We describe these turbulent processes at the interface between atmosphere and uppermost coarse blocky AL ($\sim$ 1.5 m depth) using in situ wind speed measurements and link these to measurements of eddy-covariance sensible fluxes. We then use the gained process understanding to define the appropriate “surface” temperature and humidity for the Bowen and bulk aerodynamic approaches to parameterise the turbulent fluxes.

Our work contributes to the quantitative process understanding of surface energy fluxes on a ventilated coarse blocky landform situated in a complex mountain terrain. The quantification of individual surface heat fluxes and near-surface storage terms will benefit the modelling of past and present mountain permafrost distribution and will help to anticipate the response of coarse blocky permafrost landforms to climate change.

2 Study site and past research

2.1

Murtèl rock glacier

The studied Murtèl rock glacier (WGS 84: 46°25${}^{\prime }$47${}^{\prime \prime }$ N, 9°49${}^{\prime }$15${}^{\prime \prime }$ E; CH1903+/LV95: 2'783'080, 1'144'820; 2620–2700 m a.s.l.; Fig. ) is located in a north-facing periglacial area of Piz Corvatsch in the Upper Engadine (eastern Swiss Alps), a slightly continental rain-shadowed high valley (Fig. ). Mean annual air temperature (MAAT) is $-\mathrm{1.7}$ °C; mean annual precipitation is $\sim$ 900 mm . This tongue-shaped, single-unit (monomorphic, ) active rock glacier is $\sim$ 250 m long and $\sim$ 150 m wide, surrounded by steep rock faces and directly connected to a talus slope (2700–2850 m a.s.l.). Crescent-shaped furrows ($\sim$ 3–5 m deep) and ridges with steep and, in some places, near-vertical slopes dissect the slightly north-northwestward-dipping surface ($\sim$ 10–12°) and create a pronounced furrow-and-ridge micro-topography in the lowermost part of the rock glacier. The snow cover is thicker and lasts longer in furrows than on ridges, influencing the ground thermal regime at small scale . The coarse-grained and clast-supported debris mantle is only 1–2 m thick in the colder furrows, while it is 3–5 m thick on the rest of the rock glacier. The ground ice table is accessible in a few places. Characteristic clast size ranges from 0.1 to 2 m edge length, with a few rockfall-deposited boulders of $\sim$ 3–5 m. Fine material ($\le$ sand) is virtually absent from most of the surface; its volume fraction increases with depth (inverse grading; ). Rain and percolating meltwater quickly disappear, and the surface appears dry. Beneath the coarse blocky debris mantle, roughly coinciding with the thermally defined AL, lies the perennially frozen ice-supersaturated rock glacier core. Drill cores have revealed sand- and silt-bearing massive ice (3–28 m depth, ice content over 90 % by volume), although boreholes drilled within $\sim$ 30 m distance suggest some lateral small-scale heterogeneity . Surface creep rates are $\sim$ 10 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$, show a coherent creep pattern, and have been slightly accelerating in the past decade . Water emerges seasonally from several springs or seeps at the foot of the rock glacier front and flows to the rock glacier forefield of till-veneered bedrock not underlain by permafrost (lower boundary of discontinuous permafrost) .

The site is snow covered for 7–9 months annually with an up to 1–2 m thick snowpack. Persistent katabatic winds develop in the topographically shaded cirque (no direct insolation in November–February) and redistribute snow from the windswept ridges into the furrows, eroding the snow around large blocks . These are preferential spots for snow funnels that form after the first snowfall in early winter. Oscillating airflow, resembling breathing, has been observed in these openings through the snow cover. This airflow allows for some vertical heat exchange to occur between the coarse blocky AL and the atmosphere, even during the winter months. As a result, the insulating effect of the snow cover is reduced to some degree, although the extent of this reduction has not been quantified .

Figure 1

Location of Murtèl rock glacier in the Upper Engadine, a high valley in the eastern Swiss Alps. Inset map: location and extent (black rectangle) of regional map within Switzerland (source: Swiss Federal Office of Topography swisstopo).

[Figure omitted. See PDF]

2.2 Past micro-climatological research

The Murtèl rock glacier and the Murtèl–Chastelets periglacial debris slope have been intensely investigated since ca. 1970. One of the longest continuous mountain permafrost temperature time series worldwide (since 1987) and atmospheric measurements from an automatic weather station (AWS; since 1997) run by the Swiss Permafrost Monitoring Network (PERMOS) have turned this site into a “natural laboratory” for mountain permafrost research (summarised by ). Relevant statistical and process-oriented energy balance studies are the ones by , , , , , , , , , , , and , as well as at least 10 unpublished master theses.

Apart from the two studies on Murtèl by and mentioned above, perhaps the most similar investigation to the present work in terms of ground properties is , performed on a ventilated coarse blocky permafrost site in southern Norway (Juvvasshøe); its meso-scale landscape, however, is a windswept flat mountain top dissimilar to the Murtèl cirque.

3 Measurements and data processing

3.1 Sensor placement

A total of 50 m away from the existing Murtèl PERMOS cluster (; “AWS” in Fig. c), additional sensors were installed above the ground surface and within natural cavities of the porous coarse blocky AL in August 2020. This PERMA-XT sensor cluster comprised snow and atmospheric sensors above the ground surface, AL sensors distributed in natural cavities between the blocks (Table ), and two automatic time-lapse cameras in the RGB colour and thermal infrared spectral range. We also used a four-component radiation sensor from the PERMOS cluster for our analysis; its specifications are reported in .

The above-ground sensors were located on a rock glacier ridge and included air temperature and humidity sensors, a barometer, a sonic ranger for snow height, and a sonic anemometer (CSAT; 3.87 m a.g.l.) for eddy-covariance measurements mounted on a custom-made sensor pylon. On ground level, an unheated tipping-bucket rain gauge measured liquid precipitation. Four unshielded snow thermistors at 0 (ground level), 25, 50, and 100 cm a.g.l. measured the vertical snow temperature profile. Sensor specifications are presented in Table .

The below-ground sensors were distributed in natural cavities at different micro-topographical positions (ridges, slopes, and furrows) around the meteo pylon within a 30 m distance and at different depths beneath ground surface. Five horizontally lying and vertically hanging thermistor strings and five thermo-anemometers (TP01/WS01) measured temperature and an airflow speed proxy at 5 and 30 min intervals, respectively. Most sub-surface sensors were concentrated in a 3 m deep and 0.5–1.5 m wide instrumented cavity (Fig. e). A thermistor string measured the vertical temperature profile of the cavity air (TK1/1–5), complemented by two hygrometers near the surface and in the mid-cavity (HV5; $-\mathrm{0.7}$ and $-\mathrm{2}$ m). Five thermistors (TK6/1–5) were drilled 5 cm into the blocks at depths corresponding to the TK1 thermistors. Three thermo-anemometers recorded wind speed at three levels: close to the surface ($-\mathrm{0.35}$ m), mid-cavity ($-\mathrm{1.5}$ m), and in a narrow extension at $-\mathrm{2.1}$ m (not used in this study). Finally, a back-to-back pair of pyrgeometers mounted at mid-cavity level (CGR3; $-\mathrm{1.55}$ m) measured the upward and downward long-wave radiation in the cavity. Detailed sensor specifications are presented in Table . Since accurate distances were required for the calculation of vertical gradients and fluxes, we triangulated the relative height of the sensors in the instrumented cavity with a laser distance meter and goniometer (Leica DISTO X310).

All sensors were solar-powered and wired to data loggers connected to the Internet via a mobile network. Hourly data transmission, when allowed by battery voltage, enabled timely detection of technical failures and intervention. To prevent uncontrolled power shortages during the no-insolation winter period, a protocol progressively sent power-demanding sensors into power-saving mode (50 % duty cycle or deactivated completely) as battery voltage decreased. This affected the most power-demanding sonic anemometer and, more rarely, the thermo-anemometer. Power-saving mode was active preferentially during nighttime. Such incomplete data sets are biased to daytime measurements, and thus daily average might not be representative (Sect. ). For all other sensors, power from diffuse and snow-reflected radiation was sufficient for continuous operation.

Figure 2

Sketch map of Murtèl–Corvatsch cirque (a) with two active rock glaciers, Murtèl and Marmugnun. (b) Photo of the above-surface PERMA-XT installations. (c) Wind patterns are seasonally varying in the sheltered cirque. Panels (d) and (e) show locations of sensors on Murtèl rock glacier and in the instrumented main cavity, respectively (Table ). The term “cavity roof” used throughout the text refers to the uppermost $\sim$ 1 m of the instrumented cavity.

[Figure omitted. See PDF]

Table 1

PERMA-XT sensor specifications.

Quantity [unit]	Manufacturer	Sensor type	Accuracy
Sensors above ground (atmospheric and snow sensors in Fig. d)
Air temperature $T_{\mathrm{a}}$ [°C]	CSI${}^{\mathrm{a}}$	107 temperature probe${}^{\mathrm{b}}$	$\pm \mathrm{0.01}$ °C
Relative humidity (rH for $q_{\mathrm{a}}$) [%]	CSI	HygroVUE10 hygrometer${}^{\mathrm{b}}$	$\pm \mathrm{3}$%; $\pm \mathrm{0.1}$ °C
Barometric pressure $P$ [Pa]	CSI/SETRA	CS100 barometer	$\pm \mathrm{1.5}$ hPa
Eddy-covariance flux${}^{\mathrm{c}}$	CSI	CSAT3B three-dimensional sonic anemometer
Liquid precipitation [$\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}$]	CSI	SBS500 tipping-bucket rain gauge (unheated)	$\pm \mathrm{30}$ % (undercatch)
Snow temperature [°C]	TE Connectivity${}^{\mathrm{d}}$	44031RC NTC thermistors	$\pm \mathrm{0.1}$ °C
		(0, 25, 50, 100 cm a.g.l., unshielded)
Snow height $h_{\mathrm{S}}$ [cm]	CSI	SR50A sonic ranging sensor	$max\mathit{\{}\pm \mathrm{1}$ cm, $\pm \mathrm{0.4}\mathit{\%}\mathit{\}}$
Automatic camera	MOBOTIX	M16B IP camera (RGB)
Sensors below ground (active-layer sensors in Fig. d, e)
Air temperature $T_{\mathrm{a}}\left(z\right)$ [°C]	TE Connectivity${}^{\mathrm{d}}$	44031RC NTC thermistor chain TK1/1–5	$\pm \mathrm{0.1}$ °C
Relative humidity(rH for $q_{\mathrm{sa}}$) [%]	CSI	HygroVUE5 hygrometer	$\pm \mathrm{0.3}$ °C; $\le \pm \mathrm{4}$ %
Rock temperature $T_{\mathrm{r}}\left(z\right)$ [°C]	TE Connectivity${}^{\mathrm{d}}$	44031RC NTC thermistor chain TK6/1–5	$\pm \mathrm{0.1}$ °C
		(drilled 5 cm into the blocks)
Long-wave radiation $L_{\mathrm{al}}$[$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]	Kipp & Zonen	CGR3 pyrgeometer (4.5–42 $\mathrm{\mu }$m, FoV 150°)	$<\mathrm{4}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
Airflow speed proxy [–]	Hukseflux	TP01 thermal properties sensor (formerly WS01)	$\le \pm \mathrm{25}$ % or
		(artificial leaf used as hot-film anemometer${}^{\mathrm{e}}$)	$\le \pm \mathrm{0.2}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$

Measurement range and accuracy by manufacturer/vendor. Specifications of PERMOS sensor available in and . ${}^{\mathrm{a}}$ CSI: Campbell Scientific, Inc. ${}^{\mathrm{b}}$ Sensors in radiation shield RAD10E. ${}^{\mathrm{c}}$ CSAT measures three-dimensional high-frequency sonic wind speed and derives sonic buoyancy flux. ${}^{\mathrm{d}}$ Snow temperature setup and thermistor strings manufactured by Waljag GmbH. ${}^{\mathrm{e}}$ Semi-quantitative airflow speed derived from measured heat flux.

3.2 Data processing

In the present work, we analysed data of 2 years from 1 September 2020 to 30 September 2022 except for the eddy-covariance data. The sonic anemometer (CSAT) was operational from 5 November 2020 onwards.

3.2.1 Surface radiation

All four components of the radiative heat fluxes at the surface – i.e. incoming and outgoing short-wave ($S$) and long-wave ($L$) radiation – were measured directly by the PERMOS micro-meteorological station with two back-to-back pairs of Kipp & Zonen CM3 pyranometers (0.3–3 $\mathrm{\mu }$m) and CG3 pyrgeometers (5–50 $\mathrm{\mu }$m) . We applied the same radiation corrections as in : (1) instrumental corrections of both pyranometers with nighttime short-wave radiation measurements (pyranometer offset correction), (2) correction for snow cover of the upward-looking pyranometer ($S^{\downarrow }=S^{\uparrow }/\mathrm{0.86}$ when $S^{\uparrow }>S^{\downarrow }$), (3)  correction of the incoming long-wave radiation for interference with solar radiation, and (4) rejection of incoming short-wave radiation at low sun angles (if albedo exceeds unity).

3.2.2 Eddy-covariance data

We estimated the sensible turbulent flux $Q_{\mathrm{H}}$ using the eddy-covariance method, based on the covariance of vertical wind speed fluctuations $w^{\prime }$ and temperature fluctuations $T^{\prime }$, measured by a CSI CSAT and averaged over 30 min (; see Sect. ):

1$$Q_{\mathrm{H}}^{\text{eddy}}={\mathit{\rho }}_{\mathrm{a}}{C}_p\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }},$$ where $\left({\mathit{\rho }}_{\mathrm{a}}{C}_p\right)$ refers to the volumetric heat capacity. We calculated buoyancy fluxes from the eddy raw data (10 Hz sampling rate) using Campbell Scientific's software EasyFlux™ and following the conventional pre-processing chain . Processing steps were trend and outlier removal, despiking , coordinate rotation using the double rotation method , and spectral corrections .

Due to the lack of fast-response moisture measurements (no high-frequency gas analyser), the eddy-covariance method yielded the sonic buoyancy flux $\stackrel{\mathrm{‾}}{w^{\prime }{T}_{\mathrm{snc}}^{\prime }}$, which, for $\left|Q_{\mathrm{H}}\right|>\left|Q_{\mathrm{LE}}\right|$, is close to but not equal to the sensible heat flux $\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}$ . We assess the discrepancy using the Bowen ratio method modified by based on (“SND correction”): 2$$\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}=\stackrel{\mathrm{‾}}{w^{\prime }{T}_{\mathrm{snc}}^{\prime }}{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{0.51}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}{C}_p}{L_{\mathrm{v}}\mathit{Bo}}}\right)}^{-\mathrm{1}},$$ where $T_{\mathrm{snc}}$ is the sonic temperature, ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}$ is the air temperature (averaged over 30 min), and $\mathit{Bo}$ refers to the Bowen ratio (Sect. ). $C_p$ and $L_{\mathrm{v}}$ are constants defined in Sect. .

3.2.3 Snow height and snow water equivalent

We measured snow height with a sonic distance sensor (Table ) on a rock glacier ridge in addition to the nearby PERMOS snow height measurements. Raw measurements were compensated for variations in the speed of sound with air temperature, following the manufacturer's guidelines .

We used the semi-empirical $\mathrm{\Delta }$SNOW model to convert the measured snow height to snow water equivalent (SWE). This parsimonious model requires snow height and its temporal changes as the sole input. Calibration data used to develop their model were gathered in the Swiss and Austrian Alps in climatic regions similar to the Engadine.

3.2.4 Precipitation

We measured liquid precipitation with an unheated tipping-bucket rain gauge (Table ) mounted on a rock glacier ridge.

Snow or sleet (a mix of snow and rain) can fall in any month of the year. The distinction from rain is important because of the latent heat of melting $L_{\mathrm{m}}$ that ice carries and by far exceeds the sensible heat of rainwater ($c_{\mathrm{w}}\mathrm{\Delta }T$). In the comparatively dry climate of the Engadine, wet-bulb temperature can be a better discriminator than dry-bulb temperature . Wet-bulb temperature $T_{\mathrm{wb}}$ [K] is defined in :

3$$e-e^{\ast }\left(T_{\mathrm{wb}}\right)=-\mathit{\gamma }(T-T_{\mathrm{wb}})\left[\mathrm{Pa}\right],$$ where $\mathit{\gamma }=\left(c_{\mathrm{p}}P\right)/\left(\mathrm{0.622}{L}_{\mathrm{v}}\right)$ [$\mathrm{Pa}{\mathrm{K}}^{-\mathrm{1}}$] is the psychrometric “constant”. Equation () is solved iteratively. Based on hourly images from the automatic time-lapse camera, we defined a conservative wet-bulb temperature threshold of 2 °C that separates snow or sleet ($T_{\mathrm{wb}}<\mathrm{2}\mathrm{°}\mathrm{C}$) from rain.

3.2.5 Sub-surface long-wave radiation

The in-cavity net long-wave radiation $Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}$ was calculated from the upwards $L_{\mathrm{al}}^{\uparrow }$ and downwards $L_{\mathrm{al}}^{\downarrow }$ long-wave radiation components measured with a back-to-back pyrgeometer pair installed in the instrumented cavity at a depth of 1.55 m beneath ground level (Table ):

4$$Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}=(L_{\mathrm{al}}^{\uparrow }-L_{\mathrm{al}}^{\downarrow }),$$ which is not to be confused with the long-wave radiation $L$ measured at 2 m a.g.l. (Sect. ).

We corrected raw outputs $L_{\mathrm{raw}}^{\uparrow /\downarrow }$ from the two pyrgeometers in the instrumented cavity by accounting for the long-wave radiation emitted by the instruments themselves : 5$$L_{\mathrm{al}}^{\uparrow /\downarrow }=L_{\mathrm{raw}}^{\uparrow /\downarrow }+\mathit{\sigma }{T}_{\text{CGR3}}^{\mathrm{4}},$$ where $T_{\text{CGR3}}$ refers to the pyrgeometer housing temperature. Large ($>\mathrm{0.5}$ °C) or rapid changes in housing temperature differences between the two back-to-back mounted pyrgeometers indicated dust or water deposition on the upward-facing pyrgeometer window. Such disturbed measurements appeared in the high-resolution (10 min) data but did not significantly affect the daily net long-wave radiation balance in the sheltered cavity.

3.2.6 Sub-surface airflow speed

We refer to the sub-surface “wind” in the cavity as “airflow” to differentiate it from the atmospheric wind. We used the Hukseflux WS01/TP01 sensor to perform airflow speed measurements in the AL cavities. This sensor included a heated foil that measures a cooling rate expressed as a convective heat transfer coefficient $h_{\text{WS}}$ [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$] related to airflow speed . We did not convert $h_{\text{WS}}$ to airflow speed but process this variable as a semi-quantitative indicator assuming $u\propto h_{\text{WS}}$, sufficient for our purpose. WS01/TP01 do not resolve the direction of the airflow (hence the term “speed” instead of “velocity”). Deposition and evaporation of liquid water can disturb measurements ($h_{\text{WS}}$ increase), as revealed by wrapping the heated foils in moist tissues. WS01 measurements during precipitation events were filtered out. Repeated zero-point checks were performed throughout the snow-free season by enclosing the heated foil in small, dry plastic bags for a few hours, ensuring stagnant conditions with zero airflow speed. Neither drift nor temperature dependency beyond measurement uncertainty was detected.

4 Surface energy balance calculation

4.1

Energy balance of the Murtèl near-surface AL

The point-scale surface energy balance at seasonally snow-covered sites accounts for net radiation $Q^{\ast }$ [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$], composed of short- and long-wave radiation components; turbulent fluxes, composed of sensible heat $Q_{\mathrm{H}}$ [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$] and latent heat $Q_{\mathrm{LE}}$ [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]; melt energy of snow $Q_{\mathrm{M}}$ at the surface [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]; energy from precipitation $Q_{\mathrm{P}}$ [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]; and heat flux $Q_{\mathrm{G}}$ [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$] into and from the ground when snow-free, replaced by conductive heat flux across the snow cover $Q_S$ [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$] when snow covered: 6$$\begin{array}{rl}& \underset{\text{net radiation}{Q}^{\ast }}{\underbrace{S^{\downarrow }+S^{\uparrow }+L^{\downarrow }+L^{\uparrow }}}+\underset{\text{turbulent fluxes}}{\underbrace{Q_{\mathrm{H}}+Q_{\mathrm{LE}}}}+Q_{\mathrm{M}}+Q_{\mathrm{P}}+Q_{\mathrm{G}/\mathrm{S}}\\ & =\mathrm{0}\left[\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\right].\end{array}$$ Fluxes are counted as positive if these provide energy to the reference surface, i.e. the terrain or snow surface. Unlike , we consider the reference surface as an infinitely thin skin layer without storage. Fluxes must be balanced at all times (Eq. ). Ground $Q_{\mathrm{G}}$ and snow $Q_{\mathrm{S}}$ heat fluxes measured beneath the surface are extrapolated to the reference surface by means of the calorimetric correction.

4.2 Flux parameterisations

We estimated the fluxes (terms in Eq. ) as follows: all radiative fluxes were derived from on-site measurements, i.e. net radiation $Q^{\ast }$ both above the surface (PERMOS data) and within the AL (PERMA-XT data); turbulent fluxes $Q_{\mathrm{H}}$ and $Q_{\mathrm{LE}}$ were estimated using the Bowen energy balance and the bulk aerodynamic methods and directly from eddy-covariance measurements. Snowmelt $Q_{\mathrm{M}}$ and precipitation $Q_{\mathrm{P}}$ heat fluxes were estimated using the calorimetric method from SWE estimates and on-site measured rainfall rates, respectively. The heat flux in the snowpack $Q_{\mathrm{S}}$ is the calorimetrically corrected conductive heat flux at the base of the snowpack. The ground heat flux $Q_{\mathrm{G}}$ (flux in the near-surface AL) was estimated analogously to the calorimetrically corrected net long-wave radiation measured in situ in the AL at 1.5 m depth.

4.2.1 Surface radiative fluxes

The net radiation $Q^{\ast }$ is the sum of the measured and corrected radiation components (Sect. ):

7$$Q^{\ast }=\underset{S^{\ast }}{\underbrace{S^{\downarrow }+S^{\uparrow }}}+\underset{L^{\ast }}{\underbrace{L^{\downarrow }+L^{\uparrow }}},$$ where $S^{\uparrow }$ is the outgoing short-wave radiation, $S^{\downarrow }$ is the incoming short-wave radiation, and net short-wave radiation $S^{\ast }$ is the sum of both; correspondingly, $L^{\uparrow }$ represents the outgoing long-wave radiation, $L^{\downarrow }$ is the incoming long-wave radiation, and net long-wave radiation $L^{\ast }$ is the sum of both.

4.2.2

Surface turbulent heat fluxes $Q_{\mathrm{H}}$, $Q_{\mathrm{LE}}$

We estimated the turbulent flux using three different methods: (1) the Bowen energy balance method , (2) the bulk aerodynamic method , and (3) directly with the eddy-covariance method from CSAT measurements (in the lack of a fast-response vapour analyser, only the sensible turbulent flux $Q_{\mathrm{H}}$ is estimated; Sect. ).

Bowen energy balance method

The Bowen ratio $\mathit{Bo}$ is defined as the ratio of sensible to latent heat flux and reflects the partitioning of the turbulent fluxes into sensible and latent components: 8$$\mathit{Bo}:={\displaystyle \frac{Q_{\mathrm{H}}}{Q_{\mathrm{LE}}}}={\displaystyle \frac{K_{\mathrm{h}}}{K_{\mathrm{w}}}}{\displaystyle \frac{C_p\mathrm{\Delta }T}{L_{\mathrm{v}}\mathrm{\Delta }q}}\approx {\displaystyle \frac{C_p(T_{\mathrm{a}}-T_{\mathrm{s}})}{L_{\mathrm{v}}(q_{\mathrm{a}}-q_{\mathrm{s}})}},$$ where $K_{\mathrm{h}}$ represents the eddy diffusivities for sensible heat and $K_{\mathrm{w}}$ is water vapour. Invoking the similarity principle and assuming $K_{\mathrm{h}}\approx K_{\mathrm{w}}$ , the Bowen ratio can be calculated from the isobaric specific heat capacity of air ($C_p=C_{\mathrm{pd}}(\mathrm{1}+\mathrm{0.84}q)$ with $C_{\mathrm{pd}}=\mathrm{1005}$ $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) , the latent heat of vaporisation $L_{\mathrm{v}}$ (if $T_{\mathrm{s}}\ge \mathrm{0}$ °C) or sublimation $L_{\mathrm{s}}$ (if $T_{\mathrm{s}}<\mathrm{0}$ °C) ($\mathrm{2.48}\times {\mathrm{10}}^{\mathrm{6}}$ $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$, $\mathrm{2.83}\times {\mathrm{10}}^{\mathrm{6}}$ $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$), and the gradients of temperature $\mathrm{\Delta }T$ and specific humidity $\mathrm{\Delta }q$ (specified below). The sensible and latent turbulent fluxes are then expressed as a function of the available energy $Q^{\ast }+Q_{\mathrm{M}}+Q_{\mathrm{G}/\mathrm{S}}$ ($Q_{\mathrm{P}}$ considered negligible; Sect. ): 9$$\begin{array}{rl}& Q_{\mathrm{H}}^{\mathit{Bo}}={\displaystyle \frac{Q^{\ast }+Q_{\mathrm{M}}+Q_{\mathrm{G}/\mathrm{S}}}{\mathrm{1}+\mathrm{1}/\mathit{Bo}}},\\ & Q_{\mathrm{LE}}^{\mathit{Bo}}={\displaystyle \frac{Q^{\ast }+Q_{\mathrm{M}}+Q_{\mathrm{G}/\mathrm{S}}}{\mathrm{1}+\mathit{Bo}}}.\end{array}$$

The ground and snow surface temperature $T_{\mathrm{s}}$ [°C] was calculated from the measured long-wave radiation components $L$ and the Stefan–Boltzmann law via 10$$\mathit{\epsilon }\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}=L^{\uparrow }-(\mathrm{1}-\mathit{\epsilon })L^{\downarrow },$$ where $\mathit{\epsilon }$ is the surface emissivity; $\mathit{\sigma }$ represents the Stefan–Boltzmann constant ($\mathit{\sigma }=\mathrm{5.670}\times {\mathrm{10}}^{-\mathrm{8}}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{4}}$); and $L^{\uparrow }$ and $L^{\downarrow }$ represent the measured outgoing and incoming long-wave radiation, respectively. We took the emissivity value of $\mathit{\epsilon }=\mathrm{0.96}$ from .

The surface specific humidity $q_{\mathrm{s}}$ [$\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}$] is either the specific humidity of the air in the near-surface coarse blocky AL $q_{\mathrm{sa}}$ (measured) or the saturated snow surface $q_{\mathrm{ss}}^{\ast }(T_{\mathrm{s}},P_{\mathrm{a}})$, depending on whether or not the snow cover is thick enough to suppress convective air exchange between the AL and the atmosphere (Sect. ). 11$$q_{\mathrm{s}}=\left\{\begin{array}{ll}{q}_{\mathrm{sa}},& \text{for}{h}_{\mathrm{S}}{\le }^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}\text{(snow-free or open snow cover)}\\ q_{\mathrm{ss}}^{\ast },& \text{for}{h}_{\mathrm{S}}{>}^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}\text{(decoupling snow cover)}\end{array}\right.$$ The critical snow height ${}^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}$ is reached when the AL is so strongly decoupled from the atmosphere that the convective fluxes across the snow cover are no longer detectable by our SEB estimations. We determine the snow height ${}^{\text{SEB}}{\widehat{h}}_S^{\text{crit}}$ using the AL air temperature, specific humidity, and airflow speeds (Sect. ). The specific humidity of the air $q_{\mathrm{a}}$ was calculated from the measured air temperature $T_{\mathrm{a}}$ and relative humidity at the measurement level $z_{\mathrm{m}}$ of 2 m a.g.l.

Bulk aerodynamic method

In the bulk parameterisation, sensible $Q_{\mathrm{H}}^{\text{bulk}}$ and latent $Q_{\mathrm{LE}}^{\text{bulk}}$ turbulent fluxes are driven by the gradients of temperature $\mathrm{\Delta }T=T_{\mathrm{a}}-T_{\mathrm{s}}$ [K, °C] and specific humidity $\mathrm{\Delta }q=q_{\mathrm{a}}-q_{\mathrm{s}}$ [$\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}$], respectively, and horizontal wind speed $u$ [$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$]: $$\begin{array}{}\text{12}& {\displaystyle }{Q}_{\mathrm{H}}^{\text{bulk}}& {\displaystyle }={\mathit{\rho }}_{\mathrm{a}}{C}_{p}u(T_{\mathrm{a}}-T_{\mathrm{s}})C_{\mathrm{bt}},\text{13}& {\displaystyle }{Q}_{\mathrm{LE}}^{\text{bulk}}& {\displaystyle }={\mathit{\rho }}_{\mathrm{a}}{L}_{\{\mathrm{v},\mathrm{s}\}}u(q_{\mathrm{a}}-q_{\mathrm{s}})C_{\mathrm{bt}},\end{array}$$ where ${\mathit{\rho }}_{\mathrm{a}}$ [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$] represents air density. The flux–gradient relationship is modified by atmospheric stability, accounting for enhanced turbulent fluxes in an unstable atmosphere and suppressed turbulent fluxes in a stable atmosphere, by means of the bulk exchange factors $C_{\mathrm{bt}}$ (bulk turbulent heat and vapour transfer coefficients) [unitless, –]. Different formulations of the stability functions exist. Here, we use the modified scheme , motivated by its use in the GEOtop model, a distributed hydrological model designed for complex terrain , and other studies (e.g. ). We compare the modified scheme with the iterative Monin–Obukhov scheme and the widely used Businger–Dyer scheme . The latter has been used in the previous SEB studies on Murtèl by and and in many SEB estimates on debris-covered glaciers , despite discrepancies (overestimated fluxes) at strongly unstable atmosphere noted by . The issue with the Businger–Dyer parameterisation is that near-surface wind speeds on topographically sheltered rough terrain tend to be lower for a given atmospheric stability compared with the conditions under which the empirical Businger–Dyer parameterisation was originally developed (flatland in Kansas, USA; ; ). This approach is therefore problematic in complex terrain, including our study site. Finally, we test the parameterisation without stability correction (“bulk c0”), which corresponds to the special case of a neutral atmosphere. Detailed explanations of the parameterisation schemes of the turbulent heat transfers are described in the appendix (Appendix ).

Accurate estimates of turbulent fluxes rely on representative values for the roughness lengths for momentum $z_{\mathrm{0}\mathrm{m}}$, heat $z_{\mathrm{0}\mathrm{h}}$, and moisture $z_{\mathrm{0}q}$ . We calculated the roughness length for momentum $z_{\mathrm{0}\mathrm{m}}$ from the eddy-covariance data following and : 14$$z_{\mathrm{0}\mathrm{m}}=(z_{\mathrm{u}}-d)\exp \left\{-k_{\ast }{\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{h}}}{u_{\ast }}}-{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{u}}-d}{L}}\right)\right\},$$ where ${\stackrel{\mathrm{‾}}{u}}_{\mathrm{h}}$ [$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$] represents the mean horizontal wind speed, $u_{\ast }$ [$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$] is the friction velocity, $k_{\ast }$ [0.40] is the von Kármán constant, $z_{\mathrm{u}}$ [m] is the measurement height of wind speed, $L$ [m] is the Obukhov length, $\mathit{\psi }$ is the integrated stability function (Appendix ), and $d$ is the zero-plane displacement height. We set the (unknown) zero-plane displacement height $d$ to zero . The sensor height $z_{\mathrm{u}}\left(t\right)$ is the variable distance above the ground or snow surface: $z_{\mathrm{u}}\left(t\right):=z_{\mathrm{CSAT}}-h_{\mathrm{S}}\left(t\right)$. The recommended filtering is to use only ${\stackrel{\mathrm{‾}}{u}}_{\mathrm{h}}>\mathrm{3}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ under near-neutral conditions $|z_{\mathrm{u}}/L|<\mathrm{0.05}$ . We compared our eddy-covariance-derived roughness length for momentum $z_{\mathrm{0}\mathrm{m}}$ with the aerodynamically derived values reported in . For simplicity, the scalar lengths for heat $z_{\mathrm{0}\mathrm{h}}$ and humidity $z_{\mathrm{0}q}$ are often considered equal to the momentum roughness length $z_{\mathrm{0}\mathrm{m}}$ or 1–3 orders of magnitudes smaller . Here, we could not independently estimate the scalar roughness lengths $z_{\mathrm{0}\mathrm{h}}$ and $z_{\mathrm{0}q}$ due to the lack of humidity-corrected sonic temperature or high-frequency humidity measurements. We assumed an equal roughness length for heat and moisture, $z_{\mathrm{0}\mathrm{h}}=z_{\mathrm{0}q}$, and used the unknown ratio ${\widehat{R}}_{z\mathrm{0}}:=z_{\mathrm{0}\mathrm{h}}/z_{\mathrm{0}\mathrm{m}}=z_{\mathrm{0}q}/z_{\mathrm{0}\mathrm{m}}$ as a calibration parameter.

Finally, we compared the sensible turbulent flux with the katabatic model of specifically developed for conditions of katabatic or nocturnal drainage winds. This is different from all the above-mentioned parameterisations as this predicts a quadratic rather than near-linear relation between $Q_{\mathrm{H}}$ and the driving temperature difference $\mathrm{\Delta }T$. Since we lacked vertical profile observations of potential temperature required for this parameterisation, we could not test this bulk method. Instead, we checked the validity of the quadratic $Q_{\mathrm{H}}$–$\mathrm{\Delta }T$ relation on Murtèl (cf. ) and show wind profile measurements collected by in the years 1997–2000 (Appendix ).

4.2.3

Snowmelt heat flux $Q_{\mathrm{M}}$, snow heat flux $Q_{\mathrm{S}}$, and snowpack sensible heat storage $H_{\mathrm{S}}$

The snowmelt heat flux $Q_{\mathrm{M}}$ is the heat flux consumed by the melting snowpack. We estimated $Q_{\mathrm{M}}$ from the daily change in SWE [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$] derived from the snow height data $h_{\mathrm{S}}$ (sonic ranger data) with the semi-empirical parsimonious $\mathrm{\Delta }$SNOW model : 15$$Q_{\mathrm{M}}\approx L_{\mathrm{m}}{\displaystyle \frac{\mathrm{\Delta }\left(\text{SWE}\right)}{\mathrm{\Delta }t}},\text{ if }{T}_{\mathrm{b}}=\mathrm{0}\mathrm{°}\mathrm{C},$$ where $L_{\mathrm{m}}$ [336 $\mathrm{kJ}{\mathrm{kg}}^{-\mathrm{1}}$] represents the latent heat of fusion for ice. We considered runoff-generating snowmelt to occur in spring when temperature at the base of the snowpack $T_{\mathrm{b}}$ reaches the melting point and ignored melting–refreezing events within the snowpack.

The snow heat flux $Q_{\mathrm{S}}$ is the conductive heat flux across the single-layer snowpack. We calculated ${\stackrel{\mathrm{̃}}{Q}}_{\mathrm{S}}$ in the lowermost 25 cm of the snowpack, ignoring transient effects : 16$$\begin{array}{rl}{\stackrel{\mathrm{̃}}{Q}}_{\mathrm{S}}& =-k_{\mathrm{S}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\approx -k_{\mathrm{S}}{\displaystyle \frac{T_{\mathrm{s}}-T_{\mathrm{b}}}{h_{\mathrm{S}}}}\approx -k_{\mathrm{S}}{\displaystyle \frac{T_{\mathrm{25}}-T_{\mathrm{0}}}{\mathrm{\Delta }z}},\\ & \text{ if }{h}_{\mathrm{S}}>\mathrm{30}\mathrm{cm}\text{and}\mathrm{0}\text{otherwise},\end{array}$$ where $\mathrm{d}T/\mathrm{d}z$ is a linearised temperature gradient in the lowermost 25 cm of the snowpack. To ensure that both thermistors are snow covered, the minimum snow height $h_{\mathrm{S}}$ required is 30 cm. We related the snow thermal conductivity $k_{\mathrm{S}}$ to the snow density via the empirical equation $k_{\mathrm{S}}:=k_{\mathrm{S}}\left({\mathit{\rho }}_{\mathrm{S}}\right)=\mathrm{2.93}({\mathrm{10}}^{-\mathrm{6}}{\mathit{\rho }}_{\mathrm{S}}^{\mathrm{2}}+{\mathrm{10}}^{-\mathrm{2}})$ developed for Murtèl , where we estimated the bulk density of the snowpack with the $\mathrm{\Delta }$SNOW model via ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{S}}=\text{SWE}/h_{\mathrm{S}}$ .

${\stackrel{\mathrm{̃}}{Q}}_{\mathrm{S}}$ was calculated near the base of the snowpack instead of at the snow reference surface to which Eq. () refers. With this approach, we used the more stable snow density and thermal conductivity near the snowpack base, which is less affected by compaction than the near-surface layers that receive fresh snow. The heat flux ${\stackrel{\mathrm{̃}}{Q}}_{\mathrm{S}}$ was extrapolated to the snow surface by adding the sensible heat storage changes in the snowpack $\mathrm{\Delta }{H}_{\mathrm{S}}$ above $z_{\mathrm{S}}=\mathrm{25}$ cm (“calorimetric correction”) . We estimated the changes in cold content $\mathrm{\Delta }{H}_{\mathrm{S}}$ as 17$$\mathrm{\Delta }{H}_{\mathrm{S}}\approx c_{\mathrm{S}}\left(m_{\mathrm{S}/\mathrm{A}}\right)\mathrm{\Delta }{T}_{\mathrm{S}}/\mathrm{\Delta }t\approx c_{\mathrm{S}}\text{SWE}⟨\mathrm{\Delta }{T}_{\mathrm{S}}⟩/\mathrm{\Delta }t,$$ where $\left(m_{\mathrm{S}/\mathrm{A}}\right)$ represents the mass of the snowpack per area, i.e. the SWE [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$]; $c_{\mathrm{S}}$ [$\mathrm{2.050}\mathrm{kJ}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] is the specific heat capacity of snow or ice at 0°C; and $⟨\mathrm{\Delta }{T}_{\mathrm{S}}⟩$ [°C] represents the layer-averaged snow temperature changes. The calorimetrically corrected snow heat flux $Q_{\mathrm{S}}$ was then calculated as 18$$Q_{\mathrm{S}}={\stackrel{\mathrm{̃}}{Q}}_{\mathrm{S}}+\mathrm{\Delta }{H}_{\mathrm{S}}.$$

4.2.4

Precipitation heat flux $Q_{\mathrm{P}}$

The rainfall heat flux $Q_{\mathrm{Pr}}$ was estimated via 19$$Q_{\mathrm{Pr}}=C_{\mathrm{w}}r(T_{\mathrm{P}}-T_{\mathrm{s}}),$$ where $C_{\mathrm{w}}={\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{w}}$ (4.18 $\mathrm{MJ}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{K}}^{-\mathrm{1}}$) is the water volumetric heat capacity and $r$ [${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$] is the rainfall rate intercepted at the surface. Precipitation temperature $T_{\mathrm{P}}$ was approximated using the wet-bulb temperature $T_{\mathrm{wb}}$, calculated from air temperature and relative humidity (Eq. ). We assumed precipitation in the form of rain if $T_{\mathrm{wb}}\ge \mathrm{2}$ °C. Water contributions from upslope flowing onto the rock glacier and liquid precipitation falling into the snowpack were not accounted for.

The heat flux $Q_{\mathrm{Ps}}$ due to rapid melting of sleet or shallow summertime snow ($T_{\mathrm{wb}}<\mathrm{2}$ °C, $h_{\mathrm{S}}<\mathrm{10}$ cm) was roughly estimated as 20$$Q_{\mathrm{Ps}}\sim L_{\mathrm{m}}{\mathit{\rho }}_{\mathrm{w}}{r}_{\mathrm{s}},$$ where $r_{\mathrm{s}}$ [${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$] represents the amount of solid precipitation that melts in the pluviometer per time period $\mathrm{\Delta }t$ [s]. Our measurement setup was not designed to accurately record the precipitation rate during mixed rain- and snowfall or the fraction of ice crystals and liquid water in the total precipitation. The heat flux $Q_{\mathrm{Ps}}$ is intended as an order-of-magnitude estimate.

4.2.5

Near-surface ground heat flux $Q_{\mathrm{G}}$ and sensible heat storage $H_{\mathrm{al}}^{\mathit{\theta }}$

Analogously to the calorimetrically corrected snow heat flux $Q_{\mathrm{S}}$, the ground heat flux $Q_{\mathrm{G}}$ is the sum of the measured net long-wave radiation in the instrumented cavity $Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}$ (Eq. ) and the sensible heat storage changes $\mathrm{\Delta }{H}_{\mathrm{al}}^{\mathit{\theta }}$ of the virtual layer between the ground surface and the depth of the pyrgeometer: 21$$Q_{\mathrm{G}}=Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}+\mathrm{\Delta }{H}_{\mathrm{al}}^{\mathit{\theta }}.$$

Rock mass subject to changing temperatures constitutes a heat source or sink . The sensible heat storage change in the (dry) blocks $\mathrm{\Delta }{H}_{\mathrm{al}}^{\mathit{\theta }}$ is proportional to the rate of change in rock temperature, assuming a constant volumetric heat capacity $\partial \mathit{\left\{}\right(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{al}}\left){\mathit{\rho }}_{\mathrm{r}}{c}_{\mathrm{r}}\mathit{\right\}}/\partial t=\mathrm{0}$. In the lack of rock temperatures over the 2-year time period analysed in this work, we use in-cavity air temperatures ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}$ averaged over a thermal adjustment timescale $\mathrm{\Delta }{t}_{\mathrm{r}}$ (1 d) as a surrogate (Sect. ; Appendix ). The rate of sensible heat storage and release of the blocks was then approximately calculated as 22$$\begin{array}{rl}\mathrm{\Delta }{H}_{\mathrm{al}}^{\mathit{\theta }}& =\underset{-z_{\mathrm{s}}}{\overset{\mathrm{0}}{\int }}{\displaystyle \frac{\partial }{\partial t}}\mathit{\left\{}\right(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{al}}\left){\mathit{\rho }}_{\mathrm{r}}{c}_{\mathrm{r}}{T}_{\mathrm{r}}\right(z\left)\mathit{\right\}}\mathrm{d}z\\ & \approx (\mathrm{1}-{\mathit{\varphi }}_{\mathrm{al}}){\displaystyle \frac{⟨{\mathit{\rho }}_{\mathrm{r}}{c}_{\mathrm{r}}⟩}{\mathrm{\Delta }{t}_{\mathrm{r}}}}\sum _i\mathit{\{}⟨{\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}(z_i,t+\mathrm{\Delta }{t}_{\mathrm{r}})⟩\\ & -⟨{\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}(z_i,t)⟩\mathit{\}}\mathrm{\Delta }{z}_i,\end{array}$$ where $z_{\mathrm{s}}=\mathrm{155}$ cm is the distance from the pyrgeometer pair to the ground surface (reference level); ${\mathit{\rho }}_{\mathrm{r}}$ is the rock density (2690 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) (Corvatsch granodiorite; ); $c_{\mathrm{r}}$ is the specific heat capacity (790 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$); ${\mathit{\varphi }}_{\mathrm{al}}=\mathrm{0.4}$ is the AL porosity; and $T_{\mathrm{r}}\left(z\right)$ and $T_{\mathrm{a}}\left(z\right)$ are the vertical rock and in-cavity air temperature profile [°C], respectively. In the discretised formulation, temperatures $⟨\stackrel{\mathrm{‾}}{T}\left(z_i\right)⟩$ are layer-wise averages in the $i$th layer with thickness $\mathrm{\Delta }{z}_i$ (denoted by “$⟨\cdot ⟩$”) derived from the thermistor string TK1/1 and the radiometric surface temperature $T_{\mathrm{s}}$.

5 Results

5.1 Meteorological conditions

The weather in each season differed markedly between the 2 years analysed in the present work (2020–2022; Fig. ). The winter 2020–2021 was colder than the 2021–2022 one (November–April: average temperature: $-\mathrm{6.2}$ °C vs. $-\mathrm{5.3}$ °C; minimum daily average temperature: $-\mathrm{16.5}$ °C vs. $-\mathrm{15.1}$ °C) and richer in terms of snow amount (November–April: average snow height measured on a wind-swept ridge: 76 cm vs. 54 cm) and duration (early onset of snow cover: 5 October vs. 3 November; later melt-out: mid-June vs. mid-May). Summer 2021 was cool and wet compared with the hot and dry summer 2022; temperatures were lower (July–August: average: 6.9 °C vs. 9.3 °C) with frequent passage of synoptic fronts, often bringing cold air ($\le \mathrm{3}$ °C; minimum daily average temperature: 0.7 °C vs. 5.6 °C) and mixed precipitation (sleet). Snowfall occurred in a few days throughout the summer and melted within hours. A few snow patches survived over the summer after melt-out of the winter snowpack in mid-June. In contrast, the hot and dry summer 2022 was marked by three heat waves (in June, July, and August) and daily minimum temperatures not below 5 °C. Several dry spells occurred during this season; the longest one was an 11 d long dry spell within the 5–19 July heat wave. Almost no precipitation was recorded between 20 June and 1 August despite some convective precipitation events recorded on the nearby Piz Corvatsch cable car station (MétéoSuisse station at 3294 m a.s.l.). Discharge data of the rock glacier outflow (own measurements, not shown), camera images, and field observations (fresh debris flow deposits, flooding of furrows) revealed rainwater funnelled onto the rock glacier. Data gaps for this period were filled with MétéoSuisse precipitation data from the nearby station “Piz Corvatsch”.

Wind speeds (measured by PERMOS) in the sheltered Murtèl cirque were generally low (hourly means: 1–3 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$; peak wind speed $\sim \mathrm{5}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$; Fig. c). The wind pattern was often marked by a strong diurnal cycle. Peak wind speed was reached during the night in winter (strong katabatic wind blowing downslope from southeast) and in the afternoon in summer (regional valley wind known as the “Maloja wind” blowing from west-northwest–west-southwest overruled a local anabatic wind). Summer nights were calm or with weak katabatic winds (wind speed: $\le \mathrm{2}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ from west–southeast).

Figure 3

Meteorological conditions. (a) Air temperature (daily mean) and precipitation (daily sum). (b) Snow height and SWE. Rain and sleet (mixed precipitation) separated based on a wet-bulb temperature threshold of 2 °C. Precipitation data at Piz Corvatsch from MétéoSuisse.

[Figure omitted. See PDF]

5.2 Snow and ground thermo-hygric conditions

During the snow-rich winter 2020–2021, a thick and insulating snow cover (exceeding 70 cm) sealed the cavity from the atmosphere. The cavity air was kept isothermal (within $\pm \mathrm{0.2}$ °C; Fig. ; Appendix ) and isohume at saturation at a much higher moisture level than that in the (colder) atmosphere (Fig. ; Appendix ). Virtually no dry and cold air from the atmosphere was mixed into the closed cavity system. A stable winter equilibrium temperature (WEqT; bottom temperature of snow cover (BTS, ) was reached in March 2021 ($-\mathrm{4.59}\pm \mathrm{0.05}$ °C). In contrast, during the snow-poor winter 2021–2022, the rapidly fluctuating temperature and humidity (relative and specific) indicated a connection across the thin snow cover and an exchange with the dry, cold air from the atmosphere. Following snow melt-out (end of zero curtain) and re-connection with the atmosphere, the cavity air began to warm and “desaturate” from the surface downwards (July 2021; June 2022).

In summer, the near-surface cavity roof was generally warmer and more humid compared with both the atmosphere (despite slightly lower relative humidity) and the deeper cavity (specific humidity profiles are shown in Appendix Fig. ). Daily average specific humidity in the cavity roof and that in the atmosphere were strongly correlated ($R^{\mathrm{2}}=\mathrm{0.868}$). The near-surface cavity consistently presented a moisture surplus in relation to that of the atmosphere throughout the year. This surplus persisted even during dry spells. Summertime in-cavity temperature gradients were more stable the higher the surface temperature was. In summer 2021, frequent passages of cold fronts with rapid atmospheric cooling destabilised the in-cavity air column by reducing the temperature gradients. In summer 2022, dry spells impacted the sub-surface moisture conditions down to 2 m within 3–6 d after the last precipitation event: water infiltrated within minutes, near-surface relative humidity started to decrease within hours, and the cavity air drying front (evaporation front at the isohume of rH $=$ 100 %) receded to greater depths within days. At a depth of 2 m, saturation was lost $\sim$ 5 d after the last precipitation event. Consequently, the in-cavity humidity gradients reversed, indicating a switch from downwards to upwards humidity transport. This was most pronounced in July 2022 during the intense mid-July dry spell accompanying a heat wave (Appendix Figs. , ). In contrast, near the surface, the measured gradient in specific humidity between the near-surface AL and the atmosphere and, therefore, $Q_{\mathrm{LE}}$ remained largely constant regardless of the different weather conditions and the humidity gradient within the AL (Sect. ).

Figure 4

Indicators of AL–atmosphere coupling. (a) Relationship between the daily temperature amplitude in the cavity roof (normalised by the 2 m temperature amplitude) and the corresponding snow depth $h_{\mathrm{S}}$ measured on a wind-swept rock glacier ridge (PERMA-XT station) and a broad flat area (PERMOS station). (b) Normalised airflow speed ($u/u_{\max }$) at different locations as a proxy for sub-surface ventilation and convective coupling. Measurements below the level of detection (LoD) were considered zero. Airflow speed decreases with increasing snow height relative to the maximum speeds under snow-free conditions. Onset of decoupling varies with sensor location. WS/5 is beneath a large wind-exposed block on a ridge where snow funnels remain open longer than those on flat terrain (WS/3, WS/4).

[Figure omitted. See PDF]

We determine the snow height necessary to close the snow cover and to decouple the AL from the atmosphere with the sub-surface AL air temperature, specific humidity (not shown, correlated with air temperature), and airflow speeds (Fig. ). With increasing snow depth, the differences in air temperature and specific humidity between the cavity and the atmosphere increased (Appendix Figs. , ), the correlation between in-cavity and atmospheric signal was lost, and rapid (hourly–daily) fluctuations in the in-cavity temperature and specific humidity weakened (Fig. a). Normalised near-surface AL temperature amplitudes indicate the degree of convective coupling: amplitudes similar to that in the atmosphere mean coupled and strongly attenuated amplitudes mean decoupled (Fig. a; additional winter 2022–2023 data shown). The decoupling proceeded gradually with snow depth (sketched by the schematic envelope), and the thresholds are only approximate values. Also, the sub-surface airflow speed was attenuated gradually according to the micro-topographic setting that controls the local snow depth (Fig. b): on gently sloping and less rough terrain, the vertical connection between the coarse blocky AL and the atmosphere was reduced at snow heights of 5–20 cm and lost at $\sim \mathrm{50}$–60 cm of snow (WS/3 and WS/4); on wind-swept ridges with wind erosion, a much thicker snow cover was necessary to shut down the last snow funnel ($\sim \mathrm{80}$ cm at WS/5 at the upwind side of a big block on a ridge).

5.3 Surface energy balance

Monthly SEB (Eq. ) is dominated by the short-wave $S^{\ast }$ and long-wave $L^{\ast }$ radiation components, followed by the sensible $Q_{\mathrm{H}}$ and latent $Q_{\mathrm{LE}}$ turbulent heat fluxes and the ground heat flux $Q_{\mathrm{G}}$ (Fig. ; Table ). In winter, turbulent fluxes compensate for the energy lost by net radiation. In summer, turbulent fluxes export roughly 90 % of the available net radiation. Snowmelt $Q_{\mathrm{M}}$ absorbs practically the entire net radiation ($Q_{\mathrm{M}}\approx -Q^{\ast }$) in the respective melt-out month (June 2021; May 2022), and the sensible and latent turbulent fluxes either are then small or roughly cancel each other (Fig. ). The ground heat flux $Q_{\mathrm{G}}$ is downwards (negative) during snowmelt (infiltration of meltwater into the frozen coarse blocky AL releases latent heat) and net negative during the thaw season. The sensible rain heat flux $Q_{\mathrm{P}}$ is a negligible SEB component. In short, $S^{\ast }\gtrsim L^{\ast }\gg \left|Q_{\mathrm{H}}\right|>\left|Q_{\mathrm{LE}}\right|>\left|Q_{\mathrm{G}}\right|⋙Q_{\mathrm{P}}$. The parameters and constants used in the calculations are tabulated in Appendix  (Table ).

Figure 5

(a) Monthly energy balance components. Turbulent fluxes calculated using the modified parameterisation (cL${}^{\prime }$). (b) Daily and monthly surface albedo $A$ as a snow-cover indicator.

[Figure omitted. See PDF]

Table 2

Season-averaged heat fluxes and heat-flux ratios.

Heat flux	Hydrological year 2020–2021			Hydrological year 2021–2022
[$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]	Oct–Feb	Mar–May	Jun–Sep	Oct–Feb	Mar–May	Jun–Sep
$S^{\ast }$	$\mathrm{10.9}$	$\mathrm{56.3}$	$\mathrm{155.3}$	$\mathrm{17.0}$	$\mathrm{76.2}$	$\mathrm{193.4}$
$L^{\ast }$	$-\mathrm{34.3}$	$-\mathrm{47.8}$	$-\mathrm{63.2}$	$-\mathrm{44.9}$	$-\mathrm{50.4}$	$-\mathrm{81.5}$
$Q^{\ast }$	$-\mathrm{23.4}$	$\mathrm{8.5}$	$\mathrm{91.7}$	$-\mathrm{27.9}$	$\mathrm{25.9}$	$\mathrm{111.9}$
$Q_{\mathrm{H}}$	$\mathrm{30.0}$	$\mathrm{13.4}$	$-\mathrm{31.3}$	$\mathrm{37.4}$	$\mathrm{18.9}$	$-\mathrm{71.5}$
$Q_{\mathrm{LE}}$	$-\mathrm{2.9}$	$-\mathrm{17.7}$	$-\mathrm{17.5}$	$-\mathrm{11.7}$	$-\mathrm{11.3}$	$-\mathrm{27.2}$
$Q_{\mathrm{S}}$	$-\mathrm{0.4}$	$-\mathrm{0.7}$	$\mathrm{0.0}$	$\mathrm{0.2}$	$-\mathrm{0.4}$	$\mathrm{0.0}$
$Q_{\mathrm{G}}$ ($Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}$)	$-\mathrm{0.1}$	$-\mathrm{1.5}$	$-\mathrm{5.5}$ ($-\mathrm{5.2}$)	$\mathrm{0.5}$	$-\mathrm{1.2}$	$-\mathrm{8.7}$ ($-\mathrm{9.0}$)
$Q_{\mathrm{M}}$	$\mathrm{0.0}$	$-\mathrm{4.3}$	$-\mathrm{22.6}$	$\mathrm{0.0}$	$-\mathrm{20.4}$	$-\mathrm{0.1}$
SEB imbalance	$\mathrm{3.3}$	$-\mathrm{2.3}$	$\mathrm{14.8}$	$-\mathrm{1.5}$	$\mathrm{11.5}$	$\mathrm{4.4}$
Ratios
$\mathit{Bo}:=Q_{\mathrm{H}}/Q_{\mathrm{LE}}$	$-\mathrm{10.31}$	$-\mathrm{0.76}$	$\mathrm{1.79}$	$-\mathrm{3.21}$	$-\mathrm{1.67}$	$\mathrm{2.63}$
$Q_{\mathrm{G}}/Q^{\ast }$	$\mathrm{0.004}$	$-\mathrm{0.08}$	$-\mathrm{0.06}$	$-\mathrm{0.02}$	$-\mathrm{0.05}$	$-\mathrm{0.08}$

Turbulent fluxes calculated using the modified (cL${}^{\prime }$) bulk parameterisation.

5.3.1 Surface radiation

Short- and long-wave radiation are by far the largest SEB components (Fig. a). Consequently, the well-known effect of the snow cover on albedo is the single largest control of net radiation and hence of the entire SEB (Fig. b).

In mid-winter, the meso-scale relief controls the solar radiation budget. The shaded north-facing Murtèl cirque receives no direct insolation from November to February. Depending on cloud cover and amount of incoming diffuse or terrain-reflected short-wave radiation, the net radiation is negative and dominated by the long-wave radiation budget ($L^{\ast }\ll S^{\ast }<\mathrm{0}$). Net radiation $Q^{\ast }$ was less negative in the cloudy and precipitation-rich winter 2020–2021 (December–February; average: $-\mathrm{19}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) than in the sunny and dry winter 2021–2022 ($-\mathrm{40}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) due to greater incoming long-wave radiation $L^{\downarrow }$ emitted by the clouds (233 vs. 215 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$). A case in point is the exceptionally warm and sunny November 2020, which received in total 33 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ less short- and long-wave radiation than November 2021. The monthly mean radiation deficits are up to $-\mathrm{53}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, which, together with the steep snow-covered slopes, is a favourable setting for strong katabatic winds.

5.3.2 Turbulent heat fluxes

The Bowen ratio partitions the available energy from net radiation, the ground heat flux, and the snowmelt heat flux into sensible and latent turbulent fluxes. Throughout the seasons of both years, the heat fluxes were as follows (Fig. ): during wintertime, from November to March, without direct insolation in the shaded cirque, the energy loss of 20–50 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ due to long-wave emission and terrain shading was largely compensated by the sensible heat flux $Q_{\mathrm{H}}$. Latent heat flux $Q_{\mathrm{LE}}$ at the cold snow surface was small, with mostly resublimation (moist winter 2020–2021) and sublimation fluxes (dry winter 2021–2022) within $-\mathrm{6}\pm \mathrm{10}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (Fig. a). With the latent heat flux being small and in a variable direction, the computed Bowen ratio showed a large scatter ($\left|\mathit{Bo}\right|>\mathrm{1}$; Fig. b). In the drier winter 2021–2022, $Q_{\mathrm{H}}$ compensated for the heat lost by sublimation in addition to the radiative heat loss. During spring, from April to May, before the snowmelt period, the magnitude of the fluxes remained similar but the direction reversed. In summer, after complete snow melt-out, from June/July to October, the turbulent fluxes were larger to compensate for the large radiation surplus. In the more rainy summer 2021, sensible and latent heat exports were similarly large ($\mathit{Bo}=\mathrm{1.25}\pm \mathrm{2.71}$), whereas in the drier summer 2022, heat export was more dominated by sensible flux ($\mathit{Bo}=\mathrm{1.48}\pm \mathrm{2.47}$). Towards early autumn, with decreasing temperatures, the sensible flux lost importance relative to the latent flux ($\mathit{Bo}$ decreased and approached zero). In September/October, with mixed precipitation (sleet) and first snow falls at still warm conditions near 0 °C, the latent heat export by sublimation from the snow cover intermittently dominated over the sensible heat uptake ($\mathit{Bo}$ between $\mathrm{0}$ and $-\mathrm{1}$). With further cooling, the moisture supply became limited again, and sensible heat uptake took over to offset the increasing radiation deficit ($\mathit{Bo}<-\mathrm{1}$).

We compare the estimates of turbulent fluxes (monthly averages) with the net radiation (Fig. ) and among each other (Fig. ). In the snow-free summer months, the bulk fluxes tend to slightly underestimate the fluxes, while the Businger–Dyer fluxes (cB&D) clearly overestimate them. The measured eddy-covariance flux (buoyancy flux) falls short of closing the SEB, roughly by 50 %–75 % depending on how large the (unmeasured) eddy latent flux would be. The imbalance (“missing” flux to close the SEB; ) is, in any case, substantial.

Figure 6

(a) Sensible and latent (dotted boxes) turbulent surface fluxes. Available energy from net radiation $Q^{\ast }$, snowmelt $Q_{\mathrm{M}}$, and ground heat flux $Q_{\mathrm{G}}$ (hatched boxes) is partitioned into sensible $Q_{\mathrm{H}}^{\mathit{Bo}}$ and latent $Q_{\mathrm{LE}}^{\mathit{Bo}}$ turbulent fluxes according to the Bowen ratio. The Bowen SEB is closed by design. Bulk fluxes differ according to the stability function and do not necessarily close the SEB. (b) Bowen ratio (Eq. ). November–March: large scatter in winter because $Q_{\mathrm{LE}}^{\mathit{Bo}}$ is small and changes direction (unstable calculation). June/July–October: $\mathit{Bo}\sim$ 0.5–2.5, reflecting the wet and cool summer 2021 and the dry and hot summer 2022 (Fig. ). Monthly average based on 16–31 valid values per month.

[Figure omitted. See PDF]

Figure 7

Comparison of turbulent flux estimates from different parameterisations (daily and monthly averages; see Fig.  for time series). Comparison metrics: root mean square error (RMSE [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]), mean bias error (MBE [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]), and Pearson correlation coefficient ($r$ [–]) on daily estimates.

[Figure omitted. See PDF]

5.3.3 Snowmelt heat flux

The latent heat of the melting snowpack $Q_{\mathrm{M}}$ largely consumes the available net radiation $Q^{\ast }$ in the respective snowmelt months; the snowmelt heat flux is roughly as large as typical summer sensible turbulent heat fluxes $Q_{\mathrm{H}}$ (Fig. ).

5.3.4 Precipitation heat flux

The sensible rain heat flux $Q_{\mathrm{Pr}}$ exerts a negligible cooling effect of $\le -\mathrm{2}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ at the daily timescale (not shown). Short but intense rainfall (thunderstorms) with fluxes up to $-\mathrm{100}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in 10 min is averaged out because such high-precipitation events are short. The heat fluxes $Q_{\mathrm{Ps}}$ arising from mixed precipitation or a shallow snow cover that melts within hours (“summer snow”; 10–30 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ daily average) are similar to $Q_{\mathrm{G}}$ and not negligible. Such events typically occur in early autumn (September–October) but also occurred throughout the wet and cool summer 2021 (June–October). Their timing often coincides with episodes of rapid ground cooling (Fig. ).

5.3.5 Snow and ground heat fluxes

Snow heat flux $Q_{\mathrm{S}}$ is in the range of 2 to $-\mathrm{4}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ and ground heat flux $Q_{\mathrm{G}}$ in the range of 30 to $-\mathrm{20}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ at the daily timescale (Fig. ). The snowpack is a thermal insulator. Within-snowpack conductive fluxes are 2–3 orders of magnitude smaller than typical SEB components. The conductive snow heat flux is negligible compared with the overall SEB and across-snowpack convective fluxes (snow funnels; Figs. , ). The simplifications from the single-layer snowpack that ignore vertically varying snow density do not detract from this finding.

Heat storage changes in the snowpack (Fig. a) and uppermost coarse blocky AL (Fig. b) are asymmetric in opposite directions. In the snowpack, downwards heat transfer or warming (storage gain) is more intense (larger maxima) but less frequent, likely due to non-conductive fluxes from refreezing meltwater (synchronised with warm spells $T_{\mathrm{a}}>\mathrm{0}$ °C or rapid warming to the zero curtain). In the near-surface AL, upwards heat transfer or cooling (storage loss) is more intense but less frequent. Different weather conditions in the two summers analysed are reflected by $Q_{\mathrm{G}}$: the passage of cold fronts in summer 2021 led to more frequent convective cooling, often accompanied and enhanced by sleet; conversely, ground warming is pronounced during dry spells in summer 2022. Similarly for the two winters, virtually no storage changes ($Q_{\mathrm{G}}\approx \mathrm{0}$) occurred in the snow-rich winter 2020–2021, whereas temperature and sensible heat storage fluctuations continued in the snow-poor winter 2021–2022, albeit strongly attenuated compared with those in the snow-free season (Fig. ; Sect. ). The non-zero $Q_{\mathrm{G}}$ in winter 2021–2022 beneath the thin snow cover reflects convective heat exchange across the snow funnels since the conductive snow heat flux $Q_{\mathrm{S}}$ is too small to account for $Q_{\mathrm{G}}$.

Figure 8

(a) Snow heat flux $Q_{\mathrm{S}}$ and changes in cold content of the snowpack $\mathrm{\Delta }{H}_{\mathrm{S}}$. (b) Ground heat flux $Q_{\mathrm{G}}$, net long-wave radiation in the instrumented cavity $Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}$, and sensible heat storage changes $\mathrm{\Delta }{H}_{\mathrm{al}}^{\mathit{\theta }}$ (uppermost 1.5 m of the coarse blocky AL). Beneath the insulating snow cover in winter 2020–2021, the ground thermal regime is stable, and $Q_{\mathrm{G}}\approx \mathrm{0}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. In contrast, the rapid fluctuations in $Q_{\mathrm{G}}$ during the snow-poor winter 2021–2022 are faster than those during the conduction time, indicating convective processes as a driver. Heat input was larger in summer 2022 than in summer 2021.

[Figure omitted. See PDF]

5.4 Aerodynamic roughness length

Valid values of roughness length for momentum $z_{\mathrm{0}\mathrm{m}}$ (Eq. ) scatter over 2 orders of magnitude in the range of ${\mathrm{10}}^{-\mathrm{2}}$ to ${\mathrm{10}}^{\mathrm{0}}$ m. Few (2.1 %) data points met the quality criteria $u_{\ast }>\mathrm{0.1}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ under near-neutral conditions $\left|\mathit{\zeta }\right|<\mathrm{0.05}$ or $\stackrel{\mathrm{‾}}{u}>\mathrm{3.0}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ (Fig. a). The bin-wise average (median) is 0.19 m (0.23 m) for snow-free conditions ($h_{\mathrm{S}}<\mathrm{10}$ cm) and slightly decreases with increasing snow height up to 0.07 m (0.08 m) at 90–100 cm of snow (Fig. b). Averages were calculated from the average of the logarithmised values . Snow heights between 10 and 50 cm or exceeding 100 cm are rarely observed in the 2-year data set, hence the observation gap. $z_{\mathrm{0}\mathrm{m}}$ is itself a function of sensor distance and hence of snow height $h_{\mathrm{S}}$ (Eq. ). To control for possible spurious correlation in the $z_{\mathrm{0}\mathrm{m}}$–$h_{\mathrm{S}}$ relation, we eliminated the confounding variable $h_{\mathrm{S}}$ and tested with the constant measurement height $z_{\mathrm{CSAT}}=\mathrm{3.78}$ m. This did not significantly affect the $z_{\mathrm{0}\mathrm{m}}$–$h_{\mathrm{S}}$ relation.

Figure 9

Calculated momentum roughness length $z_{\mathrm{0}\mathrm{m}}$ vs. (a) horizontal wind speed ${\stackrel{\mathrm{‾}}{u}}_{\mathrm{h}}$ and (b) snow height $h_{\mathrm{S}}$. Of all values (grey dots), 2.1 % met the quality criteria (blue dots). Momentum roughness length decreases from 0.19 to 0.07 m ($\times {\mathrm{10}}^{-\mathrm{0.43}}$) with snow height increasing from 0 to 100 cm.

[Figure omitted. See PDF]

6 Discussion

6.1 Surface fluxes and uncertainties

6.1.1 Fluxes

Monthly SEBs are closed within the calculation uncertainties of $\le \mathrm{20}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (Fig. , blue bars; Table ). This represents a substantial improvement compared with and and is due to the novel sensor array used in the present work. The largest SEB imbalances (deviation from closure) occur during the transition between seasons, when ground thermal conditions linger around 0 °C, during extreme meteorological conditions, and in mid-winter (December–March). These are the snowmelt months (June 2021, May 2022), early autumn (September), and the July 2022 heat wave and accompanying dry spell that strongly impacted the ground thermal and moisture regime (Figs. –). Larger deviations that occurred in the snow-rich mid-winter 2020–2021 are reduced by a “katabatic wind correction” for the variable anemometer height above the snow surface (Sect. ).

Figure 10

SEB imbalance (monthly averages). Turbulent fluxes calculated using modified (cL${}^{\prime }$) bulk parameterisation. Correcting excessively high wind speed measured in the low-level katabatic jet at thick snow cover improves the SEB in the snow-rich winter 2020–2021 (katabatic wind correction; Sect. , Eq. ).

[Figure omitted. See PDF]

6.1.2 Parameter sensitivity

Uncertainties arise from terrain or snow cover variability across the rock glacier, spatially variable properties of the coarse blocky AL (e.g. emissivity, porosity, intrinsic permeability), and instrumental measurement errors, among other factors. We assessed the impact of the largest sources of parameter uncertainty of the heat fluxes (Table ; sensor accuracy from ). We estimate our overall accuracy as $\sim \mathrm{20}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ at the daily timescale due to uncertainties in $Q_{\mathrm{G}}$ (uncertainties in sensible heat storage changes $\mathrm{\Delta }{H}_{\mathrm{al}}^{\mathit{\theta }}$) and $Q_{\mathrm{H}}$ (uncertainty in surface temperature $T_{\mathrm{s}}$ propagated from the emissivity $\mathit{\epsilon }$; Eq. ). These uncertainties play an important role in this study. The SEB uncertainties determine which processes are included in the SEB estimates and which are not, complementing measurement-driven criteria. We consider insignificant (i.e. within the uncertainty) the processes with associated fluxes smaller than our 20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ uncertainty threshold. The processes considered insignificant in the context of the SEB are not necessarily insignificant at depth. In fact, in the AL, all daily-averaged sub-surface heat fluxes are within 20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (Fig. b). We apply this criterion on the nocturnal Balch ventilation (Sect. ) and the decoupling snow height ${}^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}$ (Sect. ).

Insolation differences due to the local shading effect of the coarse terrain surface, the terrain slope effects, or instrumental effects (cosine response) might lead to differences of up to 30 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (daily average) in summer (Table ) . Uncertainties in the outgoing long-wave radiation $L^{\uparrow }$ of $\pm \mathrm{10}$ % corresponding to a radiometric surface temperature difference of 7 °C might arise from patchy snow cover. The lateral advective heat transport altering the boundary layer characteristics could explain the large SEB imbalance (Fig. ) of roughly 100–150 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ during the melt-out phase. Smaller but still significant $T_{\mathrm{s}}$ differences of 2 °C can occur due to local shading (micro-topography) or uncertainty in emissivity (e.g. variable snow emissivity; ) and lead to considerable $Q_{\mathrm{H}}$ uncertainties of up to 50 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. Surface temperature differences from different emissivity values ($\mathit{\epsilon }=\mathrm{1}$ instead of $\mathrm{0.96}$ from ; Eq. ) are within $-\mathrm{1.2}$ and $+\mathrm{0.25}$ °C. The deviations tend to be largest on hot clear-sky days with little incoming long-wave radiation $L^{\downarrow }$ (Eq. ) and hence might translate into considerable uncertainties in $Q_{\mathrm{H}}$ precisely when these are largest.

Temperature and specific humidity uncertainties driving bulk fluxes (Eqs. , –) lead to considerable $Q_{\mathrm{H}}$ and $Q_{\mathrm{LE}}$ uncertainties. The parameterisation is moderately sensitive to the spatially variable wind field which might arise from the micro-topography, for example wind sheltering in the furrows. The momentum roughness length $z_{\mathrm{0}\mathrm{m}}$ is varied by a factor of ${\mathrm{10}}^{\mathrm{0.5}}$ (between 0.09 and 0.45 m), which reflects the standard deviation of the measurements (Fig. ). The roughness lengths are perhaps among the most critical parameters for the estimation of turbulent fluxes when using the bulk approach (besides the emissivity $\mathit{\epsilon }$): sensitive yet hard to constrain on rough and complex mountainous terrain . Equivalence between momentum $z_{\mathrm{0}\mathrm{m}}$ and scalar roughness lengths for heat $z_{\mathrm{0}\mathrm{h}}$ or moisture $z_{\mathrm{0}q}$ would lead to prohibitively large deviations and can be excluded. As the rough terrain turns sensor height above the surface into a somewhat vague parameter, we vary it by a typical block edge length of 0.5 m. The arising $Q_{\mathrm{H}}$ deviations of 10–20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ are similar to other parameter uncertainties. Air temperature and specific humidity measurements from the PERMOS and PERMA-XT stations, located within 50 m, differ by $\mathrm{0.6}\pm \mathrm{1.0}$ °C and $\mathrm{0.004}\pm \mathrm{0.24}$ $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}$, respectively. The calculated fluxes are similar: $Q_{\mathrm{H}}$ is within $\pm \mathrm{5}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ and $Q_{\mathrm{LE}}$ within 2 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. Maximum deviations are temporarily up to 10–15 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in winter and during the snowmelt period. The PERMOS and PERMA-XT station data yield flux estimates indistinguishable within their uncertainty.

Ground heat flux $Q_{\mathrm{G}}$ primarily reflects the rate of temperature change (RTC) of the near-surface AL and is weakly sensitive to the pyrgeometer flux measurements $Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}$. The reason behind is that $Q_{\mathrm{G}}$ (Eq. ) is dominated by the sensible heat storage changes $\mathrm{\Delta }{H}_{\mathrm{al}}^{\mathit{\theta }}$ of the 1.5 m thick blocky layer above the long-wave radiation measurement . Consequently, rather large uncertainties in $Q_{\mathrm{G}}$ (20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) come from uncertainties in the storage changes $\mathit{\delta }\left(\mathrm{\Delta }{H}_{\mathrm{al}}^{\mathit{\theta }}\right)$ from two factors, namely the porosity ${\mathit{\varphi }}_{\mathrm{al}}$ and the thermal adjustment time $\mathrm{\Delta }{t}_{\mathrm{r}}$ (time lags). First, a high porosity limits the heat storage capacity that scales with $\mathrm{1}-{\mathit{\varphi }}_{\mathrm{al}}$. We tested a plausible range of $\mathit{\delta }{\mathit{\varphi }}_{\mathrm{al}}\pm \mathrm{0.2}$ and obtained uncertainties up to 20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. Porosity might laterally vary closest to the surface, precisely where daily temperature amplitudes are largest. Second, the assumption of similar air and rock temperature profiles ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}\left(z\right)\approx T_{\mathrm{r}}\left(z\right)$ – local thermal equilibrium (LTE) assumption ) – becomes problematic in the roof of the ventilated cavity, where conditions are highly transient (convective heat transfer and effect of insolation). The entire rock mass might not adjust to the rapid temperature fluctuations of $+\mathrm{2}$ to $-\mathrm{3}$ °C from day to day. The chosen thermal relaxation time $\mathrm{\Delta }{t}_{\mathrm{r}}$ of 1 d is a minimum duration (Appendix ). We assess the influence of longer adjustment times $\mathrm{\Delta }{t}_{\mathrm{r}}$ by comparing 1 d with 5 d running averages. The differences are similar to the uncertainties related to porosity. We conclude that, due to the large and rapid heat turnover in the ventilated near-surface coarse blocky AL under highly transient conditions, the ground heat flux $Q_{\mathrm{G}}$ is arguably the least-constrained flux. This might not be a surprising finding on a landform that does not present a clearly defined surface.

Finally, we let ${}^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}=\mathrm{80}$ cm to simulate an open snow cover in early–mid-winter 2020–2021 (October–February; Fig. ) and throughout the snow-poor 2021–2022. Due to the large humidity differences across the snow cover (Fig. ), this resulted in a massive (up to ${\mathrm{10}}^{\mathrm{2}}$ fold) increase in $Q_{\mathrm{LE}}$. The closing of the snow cover is a key control factor for the SEB and ground thermal regimes (Eq. ; Sect. ) .

Table 3

Instrumental sensitivity: uncertainty due to parameter values and meteorological variables.

Parameter, variable	$\mathit{\delta }{Q}^{\ast }$	$\mathit{\delta }{Q}_{\mathrm{H}}$	$\mathit{\delta }{Q}_{\mathrm{LE}}$	$\mathit{\delta }{Q}_{\mathrm{G}}$
	[$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]	[$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]	[$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]	[$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]
$S^{\ast }\pm \mathrm{10}\mathit{\%}$ (daily total)	$\le \mathrm{30}$
$L^{\uparrow }\pm \mathrm{10}\mathit{\%}$ ($L^{\uparrow }\propto \sqrt[\mathrm{4}]{T_{\mathrm{s}}}$)	25–40	50–150		0.2
$T_{\mathrm{s}}\pm \mathrm{2}$ °C ($L^{\uparrow }\pm \mathrm{3}\mathit{\%}$)	7–12	10–50	$\le \mathrm{5}$
$\mathrm{\Delta }T\pm \mathrm{1.0}$ °C		5–25		–
$\mathrm{\Delta }q\pm \mathrm{0.5}$ $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}$			5–25	–
$u\pm \mathrm{10}\mathit{\%}$ $(\le \mathrm{0.3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}})$		10	5
${\log }_{\mathrm{10}}\mathit{\left\{}{z}_{\mathrm{0}\mathrm{m}}\mathit{\right\}}\pm \mathrm{50}\mathit{\%}$		30–80	5–20
$(z_{\mathrm{0}\mathrm{h}}/z_{\mathrm{0}\mathrm{m}})\pm \mathrm{20}\mathit{\%}$		$\le \mathrm{6}$	$\le \mathrm{2}$
$z_{\mathrm{0}\mathrm{h}},z_{\mathrm{0}q}=z_{\mathrm{0}\mathrm{m}}$		$\le \mathrm{300}$	$\le \mathrm{100}$
$z_{\mathrm{u}}$, $z_{\mathrm{T}}\pm \mathrm{0.5}$ m		10–20	$\le \mathrm{10}$
$T_{\mathrm{a}}$, $q_{\mathrm{a}}$ PERMOS/PERMA-XT		5	2
$Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}\pm \mathrm{15}\mathit{\%}$				2
${\mathit{\varphi }}_{\mathrm{al}}=\mathrm{0.5}\pm \mathrm{0.2}$				$\le \mathrm{20}$
$\mathrm{\Delta }{t}_{\mathrm{r}}$ 1–5 d				$\le \mathrm{20}$
${}^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}=\mathrm{80}$ cm			30–60

Uncertainty of daily average fluxes. $Q_{\mathrm{H}}$ and $Q_{\mathrm{LE}}$ using the modified parameterisation (cL${}^{\prime }$).

In the sensitivity analysis (Table ), input parameters and meteorological variables are varied independently from each other and based on likely maximum measurement errors. However, the meso-scale relief and local sub-surface processes like ventilation might add some systematic biases that exceed instrumental errors. We will explore parameterisation uncertainties (stability corrections) in Sect.  and uncertainties in the meteorological input variables in Sect. .

6.2 Meso-scale terrain effects on turbulent flux parameterisations

6.2.1 Katabatic wind

The interaction of the steep snow-covered wintertime terrain with a negative radiation balance induces downslope katabatic winds that govern the near-surface wind field. This, in turn, affects the calculation of the bulk turbulent fluxes that require a representative wind speed as input. Our initially calculated wintertime turbulent fluxes were “too large” by 10–35 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ on monthly average (Fig. ), especially during the overcast snow-rich winter 2020–2021, while the radiative cooling and the forcing temperature deficit $\mathrm{\Delta }T$ were weaker than in the snow-poor winter 2021–2022 (Fig. ; Appendix ; Eqs. –). This result suggests that measured wind speeds were off in relation to the snow height and not to the temperature deficit. In fact, wind tower measurements on Murtèl performed by showed strong and persistent katabatic winds in winter, with a wind speed maximum a few metres above the surface weakly correlated to snow height (Appendix ; see Fig. c). The growing snow cover “shifted” the high-velocity region of the low-level katabatic jet to the level of the wind sensor, causing an apparent wind speed increase. We compensate the wind speed $u_{\mathrm{s}}$ measured at variable height above the snow cover $(z_{\mathrm{s}}-h_{\mathrm{S}})$ for the snow height with a simple power-law relation:

23$$u=u_{\mathrm{s}}(z_{\mathrm{u}}/(z_{\mathrm{s}}-h_{\mathrm{S}}){)}^{-\widehat{m}},$$ where $\widehat{m}$ is the exponent. Note that the usual correction of the sensor height above the variable snow surface further increased rather than decreased the turbulent fluxes. Since our measurements cannot resolve the near-surface wind speed profile of the low-level katabatic jet, the intention of this ad hoc “katabatic wind correction” is not to accurately describe the wintertime wind profile but rather to pragmatically render our input data amenable to the flux–gradient parameterisations based on the Monin–Obukhov theory. We use $\widehat{m}$ as a calibration parameter (denoted by the circumflex). A literature value of $\widehat{m}=\mathrm{0.6}$ for stable conditions reduced the wind speed sufficiently to yield turbulent fluxes that close the SEB within our uncertainty threshold (Fig. ). Although the Monin–Obukhov theory is not strictly valid under such conditions , our consistent flux estimates based on a scaled wind speed are in line with and , who argue that the bulk method still provides reasonable flux estimates when measured close to the surface, as was the case here (measurements within 2 m above the variable snow cover). This illustrates the importance of accurate wind speeds for the calculation of turbulent fluxes. We suggest an alternative solution in Sect. .

6.2.2 Comparison of turbulent flux parameterisations

From the comparison of the measured eddy-covariance fluxes and the calculated Bowen and bulk fluxes (daily averages), we reach the following conclusions.

Overall, the modified parameterisation (cL${}^{\prime }$) seems the best choice to parameterise the turbulent fluxes on Murtèl for three reasons. First, as expected, the modified cL${}^{\prime }$ bulk fluxes were near-identical to the iteratively calculated Monin–Obukhov fluxes (Fig. d; $r=\mathrm{0.996}$) but at less computational cost. Second, neither filtering nor post-processing was necessary because virtually all estimates met the quality criteria. Third, these show the least deviations from the Bowen fluxes (Fig. e), which we consider the benchmark as the Bowen SEB is closed by design. The Bowen fluxes offer reasonable daily estimates based on minimal measurement requirements, namely temperature and specific humidity profiles. Wind speed is not required. The minimum time resolution is daily scale. Hourly values either are numerically unstable (small $\mathrm{\Delta }T/\mathrm{\Delta }q$, especially over a “warm” spring–early-summer snow cover) or do not account for the systematic diurnal wind speed variations. In summer, for example, hourly sensible Bowen fluxes are overestimated during the calm nights and compensated by excessive daytime fluxes.

The eddy-covariance flux systematically underestimates the sensible turbulent flux by a factor of $\sim \mathrm{2}$, leading to a large imbalance (underclosure) in the eddy SEB. The error from the lacking SND correction (Sect. ) cannot account for the eddy SEB imbalance. The ratio between the eddy sensible flux and the (sonic) buoyancy flux is between 0.93 and 1.05 (10 % and 90 % quantile, respectively; Eq. ) on a daily average. Since both the air and the ground surface are most often far from saturation, the eddy sensible flux error due to the lacking high-frequency gas analyser is less than 10 % on 617 out of 692 d (89.2 %) with valid Bowen ratios. Larger deviations occur on the remaining 75 d with $\left|Q_{\mathrm{LE}}\right|>\left|Q_{\mathrm{H}}\right|$ (Bowen ratio within $-\mathrm{0.6}$ and $\mathrm{0.5}$), typically during snowmelt periods (“warm” saturated snow surface), cloudy and rainy summer days (more frequently in summer 2021 than 2022), or during the first snowfall in autumn. Applying the SND correction showed little added value. Nonetheless, the summertime eddy-covariance fluxes correlate reasonably well with the Bowen sensible turbulent fluxes (Fig. a) and the modified bulk fluxes (Fig. b) or the driving temperature gradient (Fig. a), and therefore random instrumental or data processing errors represent an unlikely explanation. We hypothesise that some hidden systematic reason causes the flux underestimation and the large eddy SEB imbalance, e.g. secondary circulation of the anabatic afternoon winds.

The empirical Businger–Dyer (cB&D) formulation, developed over flat terrain, overestimates the turbulent fluxes and deviates more strongly for larger negative fluxes (“banana” shape; Fig. f), despite extensive filtering. Even the bulk parameterisation without stability correction (c0) outperformed cB&D for the summer fluxes under unstable atmospheric conditions (Fig. c). This finding agrees well with on the Lirung debris-covered glacier that shares many topo-climatic features with the Murtèl rock glacier (unstable atmosphere over a strongly heated debris surface with anabatic valley winds).

Finally, we describe the relation between the wintertime sensible Bowen fluxes with the driving temperature gradient $\mathrm{\Delta }T$ (“surface temperature deficit”) as quadratic (Fig. b). A quadratic functional $Q_{\mathrm{H}}$–$\mathrm{\Delta }T$ relation is unique to the katabatic model of . Note that the Bowen fluxes are independent of wind speed (Eq. ). Together with the ad hoc “katabatic correction”, our data show that the turbulent fluxes in a snow-covered cirque with katabatic winds might be better parameterised by a katabatic model than by the common Monin–Obukhov bulk method, as also found by on a sloped glacier surface and discussed by in the context of atmospheric boundary layer modelling.

Figure 11

Comparison of the turbulent sensible flux estimates with the driving temperature difference $\mathrm{\Delta }T$ (a–c) and atmospheric wind speed $u$ (d–f) motivated by Eq. () (daily and monthly averages). Measured winter eddy fluxes are weakly sensitive to local $\mathrm{\Delta }T$ (a) and $u$ (d). Sensible Bowen flux shows a seasonally differing $Q_{\mathrm{H}}$–$\mathrm{\Delta }T$ relation: quasi-linear in summer and quadratic in winter (b) (with best-fit lines). Relation to wind speed $u$ also differs seasonally (e), with a clearer increase in $Q_{\mathrm{H}}$ with wind speed in summer. Note that wind speed is not an input variable for the Bowen parameterisation. The seasonally differing $Q_{\mathrm{H}}$–$\mathrm{\Delta }T$ relation reflects the seasonally differing atmospheric stability and wind conditions: stronger dependency on $\mathrm{\Delta }T$ in the unstable summer atmosphere. The quadratic relation found for the wintertime fluxes suggests that katabatic nocturnal drainage winds govern the turbulent heat transfer . Different parameterisations of the bulk fluxes (c) show different sensitivities to $\mathrm{\Delta }T$: c0 and cB&D are over-sensitive to $\mathrm{\Delta }T$ in winter and summer, respectively.

[Figure omitted. See PDF]

6.3 Micro-scale landform effects on turbulent flux parameters

The “surface” temperature $T_{\mathrm{s}}$ and “surface” specific humidity $q_{\mathrm{s}}$ are key inputs for the Bowen ratio (Eq. ) and bulk methods to estimate turbulent fluxes (Eqs. –). However, heat and moisture can be drawn from the ventilated AL beneath the surface, provided the snow cover is sufficiently thin to allow convective exchange ($h_{\mathrm{S}}{<}^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}$). Although an “open” snow cover that allows vertical convective exchange between the AL and the atmosphere (Sect. ) shapes the wintertime ground thermal regime up to a snow height of $\sim \mathrm{70}$ cm, we will argue in Sect.  and that already at more than 20 cm of snow convective exchange across the snow cover can be ignored in the context of our SEB.

6.3.1 Snow cover and AL–atmosphere coupling

The snow cover controls the coupling between the coarse blocky AL and the atmosphere. The turbulent fluxes draw heat and moisture from the AL as long as the snow cover is “open” via snow funnels, a phenomenon widely observed on coarse blocky landforms . A thickening snow cover gradually suppresses the vertical convective coupling between the AL and the atmosphere, as more snow funnels close (sketched by the schematic envelopes in Fig. ). A snow height of $\sim \mathrm{10}$ cm begins to decouple the coarse blocky AL from the atmosphere (“semi-closed” in Fig. ), but much more snow ($\sim \mathrm{70}$ cm) is necessary to achieve the insulating effect of the snow cover on the ground thermal regime (“closed/insulating”). Then, the snow cover is thick enough to suppress convective air exchange, rapid sub-daily fluctuations are strongly attenuated, and large temperature and moisture gradients to the outside air build up , indicating that the AL–atmosphere coupling is weak. Such a high value is typical for a terrain as rough and blocky as on Murtèl, agreeing with, for example, and . The insulating effect of a thick snow cover was already observed decades ago and used to indirectly map the permafrost distribution using the bottom temperature of snow cover (BTS) method .

As regards the SEB, a snow cover as thin as ${}^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}=\mathrm{20}$ cm causes a decoupling that is strong enough for our SEB to be significant. This is shown in the snow-poor winter 2021–2022 with snow depths between $\mathrm{40}$ and 70 cm, when the ground heat flux $Q_{\mathrm{G}}$ kept fluctuating but remained below the 20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ uncertainty threshold (Table ). The convective flux across the snow cover cannot deviate strongly from $Q_{\mathrm{G}}$ because it was the single largest heat flux to supply/extract the heat for the AL sensible heat storage changes (small conductive snow heat flux, $Q_{\mathrm{S}}\pm \mathrm{2}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, $T_{\mathrm{s}}<\mathrm{0}$ °C, and no latent effects from snowmelt in that period). We take this threshold as the critical snow thickness (${}^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}$; Eq. ). We emphasise that the 20 cm threshold is an operational definition for the “decoupled” snow cover in the context of our SEB and its uncertainties.

6.3.2 “Surface” temperature and near-surface ventilation

From Sect.  it follows that for the decoupling or insulating snow cover ($h_{\mathrm{S}}{\ge }^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}$), the reference surface coincides with the snow surface. The radiometric surface temperature represents the surface temperature $T_{\mathrm{s}}$ (Eqs. , ). Whenever the snow cover is open or under snow-free conditions ($h_{\mathrm{S}}{<}^{\text{SEB}}{\widehat{h}}_S^{\text{crit}}$), the radiometric surface temperature is the input meteorological variable $T_{\mathrm{s}}$ (Eq. ). The surface is where the solar radiation is intercepted and transformed to thermal energy, hence the main source of sensible heat for $Q_{\mathrm{H}}$. The large gradient $(T_{\mathrm{a}}-T_{\mathrm{s}})/z_{\mathrm{T}}$ drives the sensible turbulent flux $Q_{\mathrm{H}}$ (Eqs. , ; Fig. a) and the wind-forced ventilation (Fig. b) . This is shown by the measured eddy-covariance flux (Fig. c) that largely follows the difference between air and radiometric surface temperatures (Fig. a, red area), as do the local (anabatic) wind in the atmosphere above the surface and the airflow speed in the near-surface AL (Fig. b). During daytime with strong heating on the low-albedo surface, the blocky surface exceeds air temperatures by $>\mathrm{15}$ °C.

However, we found evidence of near-surface ventilation that cannot be parameterised by the 2 m air temperature and the radiometric ground surface temperature: nocturnal Balch ventilation, a nighttime cooling process owed to the interplay of the air permeability, and the thermal inertia (heat retention) of the coarse blocky AL that is most apparent during fair-weather summer days with clear nights. Upwards turbulent heat export is protracted into the evening hours, long after sunset ($S^{\ast }$; Fig. c), as shown by the measured in-cavity airflow speed (WS/3; Fig. b) and the eddy-covariance flux (Fig. c). Furthermore, the measured eddy-covariance flux remained upwards-directed in the early morning hours, when the terrain surface has radiatively cooled to air temperature $T_{\mathrm{s}}\approx T_{\mathrm{a}}$ or even during reversals $T_{\mathrm{s}}<T_{\mathrm{a}}$ in clear-sky nights (that frequently occur in the Engadine: e.g. 24 out of 31 nights in August 2022). The calculated bulk fluxes driven by 2 m air temperature and radiometric ground surface temperature (rGST; red area in Fig. a; Eq. ) are then close to zero (or even downwards) and do not fully capture this nocturnal drawing of heat from the near-surface AL. The large thermal inertia of the rock mass and the protection from long-wave radiative cooling stabilise the sub-surface air temperatures and keep the air in the roof of the ventilated cavity during the nights warmer than that in the atmosphere outside and that on the ground surface (Fig. a; $\text{TK1/1}\approx \text{TK6/1}>T_{\mathrm{a}}\approx T_{\mathrm{s}}$ by $\sim \mathrm{4}$ °C). In the nighttime, the locally stably stratified air in the uppermost cavity roof is hence warmer and more unstable than that in the atmosphere (non-local static stability, ). The high permeability of the coarse blocky AL allows this air to escape upwards into the atmosphere or, equivalently, allows the colder outside air to sink into the cavity, thereby exporting heat (Balch ventilation; ). describes an occasional nocturnal sub-surface ventilation on Murtèl when $T_{\mathrm{s}}<T_{\mathrm{a}}$. Also interpreted nocturnal near-surface ventilation on the debris-covered Koxkar glacier (Xinjiang, China) from eddy-covariance and temperature data.

We assess this uncertainty by using the cavity roof temperature TK1/1 instead of the conventional 2 m air temperature as input $T_{\mathrm{a}}$ for the bulk fluxes (Eq. ) and compare the nighttime fluxes (when TK1/1 $>T_{\mathrm{s}}$). Due to the low wind speeds at night, the estimated nocturnal Balch fluxes are small despite the appreciable temperature gradients (10–20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$), maximally 10 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ larger than the “conventional” bulk flux and within our SEB 20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ uncertainty (Table ). Since the short-wave radiative forcing $S^{\ast }$ and the wind speed covary in phase, the “conventional” radiometric ground surface temperature and 2 m air temperature measurements are sufficiently adequate to parameterise the turbulent fluxes when sub-daily resolution is not required – especially considering that the chosen parameterisation (stability function) shows a much larger effect on the calculated turbulent flux (Fig. c; compare cL${}^{\prime }$ with cB&D). We reach a conclusion equivalent to that in Sect.  for the critical snow height: the nocturnal convective processes in summer are within the uncertainties of the daily-averaged SEB but might exert an important cooling effect on the sub-surface energy balance.

Figure 12

Hourly averages of (a) temperatures, (b) atmospheric wind speed, (c) near-surface AL airflow speed, and (d) sensible turbulent fluxes and net short-wave radiation $S^{\ast }$ during a summer fair-weather period (15–19 July 2022). The cavity air is stably stratified (TK1/1 $>$ TK1/3), but the air in the cavity roof is warmer and unstable compared with the outside air. During clear summer nights, the atmospheric air and the blocky surface cool down more strongly than the air in the near-surface coarse blocky AL that receives heat from the blocks (air temperature TK1/1 approaches the rock temperature TK6/1 in the night). Conversely to the calculated bulk fluxes, the measured eddy flux decays slowly in the evening–midnight (18–24 h) and remains upwards (negative) despite the small temperature gradient in the early morning (3–7 h, a). The turbulent sensible flux draws heat from the near-surface AL that is nocturnally warmer than $T_{\mathrm{a}}$ and thus unstable. (b, c) Atmospheric wind and AL airflow speeds co-vary in phase (but at $\sim$ $\mathrm{10}\times$ smaller magnitudes), follow the net radiation with some delay, and are closely related to the turbulent fluxes (d). The specific humidity shows small diurnal oscillations; the associated $Q_{\mathrm{LE}}$ varies less than $Q_{\mathrm{H}}$ in absolute terms.

[Figure omitted. See PDF]

6.3.3 “Surface” humidity and evaporation

From Sect.  it follows that for the decoupling or insulating snow cover ($h_{\mathrm{S}}{\ge }^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}$), the snow surface humidity at saturation $q_{\mathrm{ss}}^{\ast }$ represents the surface specific humidity $q_{\mathrm{s}}$ (Eqs. , ) . Whenever the snow cover is open or under snow-free conditions ($h_{\mathrm{S}}{<}^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}$), the humidity measurement in the cavity roof $q_{\mathrm{sa}}$ is the input meteorological variable $q_{\mathrm{s}}$ (Eq. ).

Perhaps contrary to the impression of a dry-looking blocky surface, the estimated summertime latent turbulent flux $Q_{\mathrm{LE}}$ is $\sim$ 30 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, corresponding to an evaporation rate of $\sim$ 1 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{d}}^{-\mathrm{1}}$, even during the severe July 2022 dry spell. This value is in the range of evaporation rates reported for debris-covered glaciers of 0.6–2.8 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{d}}^{-\mathrm{1}}$ . However, note that our evaporative flux estimate $Q_{\mathrm{LE}}$ might be an upper bound, in particular during the July 2022 dry spell. The moisture source for evaporative $Q_{\mathrm{LE}}$ is in the AL, not at the surface (except during precipitation events). However, our bulk parameterisation ($C_{\mathrm{bt}}$; Eq. ) ignores the additional resistance to vapour transport imposed by the blocky layer between the specific humidity measurement in the cavity roof and the atmosphere ($q_{\mathrm{sa}}$ is measured 0.7 m beneath the terrain surface; Fig. e). The neglected resistance to vapour transport in the coarse blocky AL during moisture-limited evaporation stages, when moisture is drawn from deeper levels and governed by vapour diffusion through the porous AL , might lead to an overestimation of $Q_{\mathrm{LE}}$. The longer the dry spell lasts, the deeper in the AL the moisture is drawn from (Sect. ), and the more important this effect of vapour transport resistance becomes. This likely explains why the specific humidity in the cavity roof is almost always (down to sub-hourly timescales) higher than that in the atmosphere despite the relentless mixing by ventilation (moisture surplus relative to the atmosphere; Sect. ), and it possibly accounts for the negative July 2022 SEB imbalance (Fig. ). The upwards vapour transfer in the comparatively large and strongly ventilated cavity during the dry spell might be overly efficient compared with that in the surrounding AL. This might have led to a non-representatively high specific humidity $q_{\mathrm{sa}}$ in the cavity roof and an overestimated evaporative flux $Q_{\mathrm{LE}}$.

6.3.4 Aerodynamic roughness lengths

Our mean (median) values of the filtered aerodynamic momentum roughness lengths $z_{\mathrm{0}\mathrm{m}}$ of 19 to 7 cm (23–8 cm) agree with the values of and of 18 and 7 cm for the snow-free and snow-covered Murtèl surface, respectively (Fig. ). The scatter range of $\sim$ 2 orders of magnitude is similar to that in other studies . The absolute values are at the upper limit of what has been typically found on debris-covered glaciers or on Juvvasshøe by (5 cm), which is plausible given the rough terrain of the Murtèl rock glacier.

The calculated roughness length decreases slightly with increasing snow height by 0.1 cm per cm of snow (Fig. b). A similar but much stronger relation between snow height and roughness lengths was found on the Haut Glacier d'Arolla, where found $z_{\mathrm{0}\mathrm{m}}$ to decrease by 2 orders of magnitude with snow heights up to 3 m. This represents a much stronger relation than that on Murtèl, where $z_{\mathrm{0}\mathrm{m}}$ decreases by $\sim \mathrm{0.5}\times$ over the observed range of snow thickness. With a maximum snow height of 120 cm (PERMA-XT measurement; Fig. ) or 200 cm (PERMOS) during the measurement period, the snow cover is thin compared with the terrain roughness (edge length of blocks $\sim$ 50 cm) and undulations. Both the comparatively high roughness length and the weak sensitivity to snow height might suggest that the aerodynamic roughness is largely controlled by the furrow-and-ridge micro-topography that is not smoothed out by the snow cover or by the few largest blocks that stick out of the snow cover rather than by the average block size. In forests, for comparison, the tallest trees of the canopy have a disproportionally large influence on the aerodynamic roughness .

For the bulk parameterisation, we linearly interpolate ${\log }_{\mathrm{10}}\mathit{\left\{}{z}_{\mathrm{0}\mathrm{m}}\mathit{\right\}}$ with snow height. Due to the lack of accurate humidity-corrected sonic temperature measurements, we cannot calculate $z_{\mathrm{0}\mathrm{h}}$ or $z_{\mathrm{0}q}$. We used the unknown ratio ${\widehat{R}}_{z\mathrm{0}}:=z_{\mathrm{0}\mathrm{h}}/z_{\mathrm{0}\mathrm{m}}=z_{\mathrm{0}q}/z_{\mathrm{0}\mathrm{m}}$ as a calibration parameter and found ${\widehat{R}}_{z\mathrm{0}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{2}}$. While the approximation $z_{\mathrm{0}\mathrm{h}}\approx z_{\mathrm{0}q}$ seems applicable, $z_{\mathrm{0}\mathrm{m}}\approx z_{\mathrm{0}\mathrm{h}}$ is not compatible with our parameterisation (Table ). This agrees with previous studies on debris-covered glaciers and hummocky ice surfaces .

6.4 Synthesis

The thick coarse blocky AL strongly insulates the underlying permafrost body. We quantified this well-known effect on Murtèl: during the thaw seasons 2021 and 2022, roughly 90 % of the net radiation $Q^{\ast }$ was exported by the turbulent fluxes and not available to melt ground ice (Fig. ). The ratio between received surface net radiation $Q^{\ast }$ to downwards-transmitted heat flux $Q_{\mathrm{r}}=Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}$ is $Q^{\ast }:Q_{\mathrm{r}}\approx \mathrm{9}:\mathrm{1}$. This follows from the ratio of the measured in-cavity net long-wave radiation $Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}$ to the surface net radiation $Q^{\ast }$ and the relation $Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}=-\mathrm{0.12}{Q}^{\ast }+\mathrm{5.6}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ ($R^{\mathrm{2}}=\mathrm{0.680}$). Our measurements corroborate the relation between surface net radiation and ground heat flux that has been found by using PERMOS data from 1997–2019. Export of the received net radiation is predominantly by sensible fluxes $Q_{\mathrm{H}}$ (50 %–70 %) and secondarily by latent fluxes $Q_{\mathrm{LE}}$ (30 %–50 %; Table ). Surface albedo (spring–early summer snow cover) and sub-surface thermal and moisture regime of the thick coarse blocky AL control the energy partitioning at/near the surface and the rock glacier's efficiency/ability to export the heat supplied by the surface net radiation $Q^{\ast }$, the primary heat input (Fig. ; Table ).

More heat was transferred into the ground (and more ground ice observed to melt) during the hot and dry summer 2022 than during the cool and wet summer 2021. Two sets of conditions enhanced the downwards heat transfer $Q_{\mathrm{G}}$ in the hot and dry summer 2022 as measured by the sub-surface pyrgeometer, in particular during dry spells and heat waves: first, less snow in winter 2021–2022 and a warm 2022 spring resulted in an early snowmelt in May, 1 month earlier than in summer 2021, and an early start of the thaw season. This exposed the dark low-albedo blocky surface to the strong insolation and the June heat wave, resulting in a rapid rise in ground temperatures and hence an increase in sensible heat storage and downwards ground heat flux (Fig. ; see , for a 20-year perspective). Second, the 11 d dry spell in July 2022 exhausted the near-surface moisture stores in the coarse blocky AL as indicated by the specific humidity deficit $(q_{\mathrm{sa}}^{\ast }-q_{\mathrm{sa}})$ in the cavity roof (Appendix Fig. ). Lack of moisture limited the evaporative cooling $Q_{\mathrm{LE}}$ and countered the rock glacier's ability to export the heat (Fig. ). The moisture storage capacity of the coarse blocky Murtèl AL is limited because the near-surface AL presents a rapid drainage and little fine material (silt) to hold water. During dry spells, the evaporation front recedes quickly to greater depths.

Figure 13

Relation between measured 10 d averaged surface net radiation $Q^{\ast }$ and measured 10 d averaged in-cavity long-wave radiation $Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}$ (of the consecutive 10 d window) summarises the relation between $Q^{\ast }$, $Q_{\mathrm{H}}$, and $Q_{\mathrm{LE}}$. The in-cavity net long-wave radiation $Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}$ represents the radiative downward heat transfer $Q_{\mathrm{r}}$ towards the ground ice table (note that upward heat transfer is primarily by convection and not represented by the radiation measurements; $Q_{\mathrm{CGR}\mathrm{3}}^{\mathrm{rad}}=Q_{\mathrm{r}}\approx \mathrm{0}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ during cooling episodes). Of the available net radiation $Q^{\ast }$, 90 % is exported into the atmosphere by $Q_{\mathrm{H}}+Q_{\mathrm{LE}}$. Only 10 % is transmitted deeper into the ground and available to warm the AL and to melt ground ice. Latent turbulent heat export is less efficient during dry spells, when the coarse blocky AL dries out and evaporation becomes increasingly moisture-limited (outliers from June, July 2022). The relations between in-cavity long-wave radiation and air temperature, surface temperature ($\propto L^{\uparrow }$), or $Q_{\mathrm{H}}$ are qualitatively similar.

[Figure omitted. See PDF]

Efficient evaporative cooling relies on frequent precipitation events to resupply, as occurred in summer 2021. With climate change, both factors – early start of the thaw season and hot and dry weather spells in summer – are projected to worsen. A shift to earlier snowmelt and an increase in the number of snow-free days has already been observed in the Swiss Alps and is projected to continue with climate change . Also, the frequency, duration, and intensity of heat waves and dry spells are likely to increase. Changes in the SEB of the thermally conditioned rock glaciers and other mountain permafrost landforms entail changes in the ground thermal regime and ground ice content first and, ultimately, morphological changes: thawing, melting of ground ice, and degradation.

7 Conclusions

We estimated the year-round surface energy balance (SEB) of the seasonally snow-covered ventilated coarse blocky active layer (AL) of the active Murtèl rock glacier that is situated in a cirque in the Upper Engadine (eastern Swiss Alps). The meso-scale landscape produces seasonally contrasting atmospheric conditions of downslope katabatic jets in winter and a strongly unstable atmosphere in summer. At landform scale, the ventilated near-surface AL acts as a buffer layer where heat and moisture transfer is coupled to the atmosphere, unless sealed by a thick snow cover. Based on a novel sensor array located above the ground surface and in the AL that expands an on-site automatic weather station from the Swiss Permafrost Monitoring Network (PERMOS), we were able to improve previous SEB calculations on Murtèl by and . The measurement period is from September 2020 to September 2022. Our monthly SEB imbalances are within 20 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ except during the snowmelt months. The main findings concern (i) the climate resilience of the Murtèl rock glacier and (ii) technical aspects of the turbulent flux parameterisations.

• i.

  Two crucial factors for climate resilience of the Murtèl rock glacier are the insulating high-albedo snow cover and a near-complete energy turnover during the snow-free thaw season. The two meteorologically contrasting years studied in this work with the cool and wet summer 2021 and the hot and dry summer 2022 were traced into the SEB and ground thermal and moisture regime. First, an early snowmelt in 2022 prolonged the thaw season, exposed the dark blocky surface to the intense July insolation, and increased AL temperature gradients and the downward heat flux. Second, about 90 % of the received surface net radiation was exported and only 10 % effectively transferred towards the ground ice table and available to melt ground ice. Heat export occurs predominantly via sensible turbulent flux and secondarily via the latent turbulent flux from evaporation. The degree of energy turnover and the turbulent flux partitioning is co-controlled by the availability of moisture for evaporation. Dry spells and heat waves counter the rock glacier's ability to export heat by limiting evaporative fluxes leaving the rapidly drying coarse blocky AL, since the moisture storage capacity is limited. Both trends – earlier snowmelt and more heat waves and dry spells – are projected to continue with climate change. This potentially renders coarse blocky landforms vulnerable to heat waves and dry spells.

• ii.

  With our in-mountain permafrost unprecedentedly comprehensive in situ measurements of eddy-covariance flux, liquid precipitation, snow height, AL temperature and humidity profiles, sub-surface airflow speeds, and sub-surface long-wave radiation, we tested different bulk parameterisations and constrained their input parameters.

  • We parameterised the year-round turbulent fluxes using the modified scheme despite seasonally contrasting atmospheric stability and wind profiles. Wintertime katabatic wind speeds needed to be scaled to close the SEB, which hints at the limits of parameterisations based on the Monin–Obukhov similarity theory in complex mountain terrain and katabatic drainage winds.

  • Sensible and latent partitioning of the turbulent fluxes using the simple Bowen ratio approach agreed with the bulk fluxes on monthly to daily resolution. Given its parsimony and independence from atmospheric stability and wind speed, the Bowen approach represents a valuable and robust tool to estimate the turbulent fluxes in complex terrain with a strongly heterogeneous wind field, provided sub-daily resolution is not required.

  • Since daily oscillations of wind speed and insolation are nearly in phase, the sensible turbulent flux $Q_{\mathrm{H}}$ is driven by the gradient between radiometric ground surface temperature $T_{\mathrm{s}}$ and air temperature $T_{\mathrm{a}}$. The heat flux $Q_{\mathrm{H}}$ from nocturnal ventilation of the permeable coarse blocky AL is small due to low wind speeds at night. Using a radiometric surface temperature is convenient for modelling $Q_{\mathrm{H}}$ with remotely sensed data.

  • During the thaw season, the evaporative turbulent flux $Q_{\mathrm{LE}}$ is 30 %–50 % of the sensible turbulent flux as a function of moisture availability in the near-surface AL. During dry spells, the near-surface moisture stores are exhausted within days, and moisture supply to the surface is limited by the convective vapour transport from the deeper AL.

  • Measured eddy-covariance fluxes were systematically too small to close the SEB. Still, the eddy-derived momentum roughness length $z_{\mathrm{0}\mathrm{m}}$ corroborates previous aerodynamic estimates from . $z_{\mathrm{0}\mathrm{m}}$ varies seasonally between 7 and 20 cm as a function of snow height that smooths the landscape.

Appendix A Turbulent flux parameterisations

Different formulations of the stability functions exist. Here, we compare the widely used Businger–Dyer relations with the formulation from and the iterative Monin–Obukhov scheme.

A1 Businger and Dyer parameterisation

In the Businger and Dyer parameterisation (denoted by “cB&D”), the bulk transfer coefficient $C_{\mathrm{bt}}$ is expressed as

A1$$C_{\mathrm{bt}}={\displaystyle \frac{k_{\ast }^{\mathrm{2}}}{\ln \left\{\frac{z_{\mathrm{u}}}{z_{\mathrm{0}\mathrm{m}}}\right\}\ln \left\{\frac{z_{T/q}}{z_{\mathrm{0}h/q}}\right\}}}({\mathit{\varphi }}_{\mathrm{m}}{\mathit{\varphi }}_{h/q}{)}^{-\mathrm{1}},$$ where the Businger–Dyer non-dimensional stability functions for heat and water vapour are expressed as a function of the bulk Richardson number ${\mathit{Ri}}_{\mathrm{b}}$ (Eq. ) , A2$$\begin{array}{rl}({\mathit{\varphi }}_{\mathrm{m}}{\mathit{\varphi }}_{\mathrm{h}}{)}^{-\mathrm{1}}& =({\mathit{\varphi }}_{\mathrm{m}}{\mathit{\varphi }}_q{)}^{-\mathrm{1}}\\ & =\left\{\begin{array}{ll}(\mathrm{1}-\mathrm{5}{\mathit{Ri}}_{\mathrm{b}}{)}^{\mathrm{2}},& \text{for}\mathrm{0.0}\le {\mathit{Ri}}_{\mathrm{b}}\le \mathrm{0.2}\text{(stable case)}\\ (\mathrm{1}-\mathrm{16}{\mathit{Ri}}_{\mathrm{b}}{)}^{\mathrm{0.75}},& \text{for}{\mathit{Ri}}_{\mathrm{b}}<\mathrm{0.0}\text{(unstable case)}\end{array}\right.\end{array}$$ that corrects for non-neutral stability of the near-surface atmosphere (atmospheric stratification). These are unity in the case of neutrally stable atmosphere (${\mathit{Ri}}_{\mathrm{b}}=\mathrm{0}$), denoted by “c0”. Used constants and parameters are the von Kármán constant $k_{\ast }$ (0.40); the roughness lengths for momentum $z_{\mathrm{0}\mathrm{m}}$, heat $z_{\mathrm{0}\mathrm{h}}$, and water vapour $z_{\mathrm{0}q}$; and the (constant) measurement heights for wind speed $z_{\mathrm{u}}$ [m], temperature $z_{\mathrm{T}}$, and humidity $z_q$ above the ground surface. In this work, the variable sensor height above the snow surface is accounted for by correcting the wind speed $u$ rather than simply subtracting the snow height from the sensor distance above the ground for site-specific reasons discussed in Sect. .

The Businger–Dyer and the modified parameterisations characterise atmospheric stability with the bulk Richardson number ${\mathit{Ri}}_{\mathrm{b}}$ defined by A3$${\mathit{Ri}}_{\mathrm{b}}:={\displaystyle \frac{g}{\stackrel{\mathrm{‾}}{T}}}{\displaystyle \frac{\frac{T_{\mathrm{a}}-T_{\mathrm{s}}}{z_{\mathrm{T}}-z_{\mathrm{0}\mathrm{h}}}}{{\left(\frac{u}{z_{\mathrm{u}}-z_{\mathrm{0}\mathrm{m}}}\right)}^{\mathrm{2}}}},$$ with the average temperature $\stackrel{\mathrm{‾}}{T}:=(T_{\mathrm{a}}+T_{\mathrm{s}})/\mathrm{2}$ [K].

However, this widely used stability function yields implausible fluxes when the atmosphere strongly deviates from near-neutral stability conditions, $\left|{\mathit{Ri}}_{\mathrm{b}}\right|>\mathrm{0.2}$. Equation () is invalid for ${\mathit{Ri}}_{\mathrm{b}}>+\mathrm{0.2}$ and tends to diverge for highly unstable atmospheres at low wind speeds where ${\mathit{Ri}}_{\mathrm{b}}\to -\mathrm{\infty }$. Filtering out these fluxes as proposed by potentially deletes a sizeable portion of the calculated fluxes at the wind-sheltered Murtèl cirque.

A2

Modified scheme

The scheme is an analytical approximation of the iterative Monin–Obukhov scheme. The bulk exchange factor is expressed as $C_{\mathrm{bt}}=f_{\mathrm{h}}{C}_{\mathrm{Hn}}$ (notation from ), where $C_{\mathrm{Hn}}=k_{\ast }^{\mathrm{2}}\ln \mathit{\{}{z}_{\mathrm{u}}/z_{\mathrm{0}}{\mathit{\}}}^{-\mathrm{2}}$ (the neutral exchange coefficient) and A4$$f_{\mathrm{h}}=\left\{\begin{array}{ll}(\mathrm{1}+\mathrm{10}{\mathit{Ri}}_{\mathrm{b}}{)}^{-\mathrm{1}},& \text{for}{\mathit{Ri}}_{\mathrm{b}}\ge \mathrm{0}\text{(stable case)}\\ \mathrm{1}-\mathrm{10}{\mathit{Ri}}_{\mathrm{b}}(\mathrm{1}+\mathrm{10}{C}_{\mathrm{Hn}}\sqrt{-{\mathit{Ri}}_{\mathrm{b}}}/f_z{)}^{-\mathrm{1}},& \text{for}{\mathit{Ri}}_{\mathrm{b}}<\mathrm{0}\text{(unstable case)}\end{array}\right.$$ with A5$$f_z={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\sqrt{{\displaystyle \frac{z_{\mathrm{0}}}{z_{\mathrm{u}}}}}.$$ However, one shortcoming of the scheme critical on rough surfaces is the assumption of equal momentum roughness length ($z_{\mathrm{0}\mathrm{m}}$) and sensible heat roughness length ($z_{\mathrm{0}\mathrm{h}}$), i.e. $z_{\mathrm{0}\mathrm{m}}=z_{\mathrm{0}\mathrm{h}}=z_{\mathrm{0}}$. Hence, we use the scheme modified by (denoted by “cL${}^{\prime }$”). The implementation is described in and .

A3 Iterative Monin–Obukhov scheme

In the Monin–Obukhov scheme (denoted by “cMO”), the atmospheric stability is characterised by the Obukhov length $L$ defined by

A6$$L={\displaystyle \frac{u_{\ast }^{\mathrm{3}}}{k_{\ast }\frac{g}{T_{\mathrm{a}}}\frac{Q_{\mathrm{H}}}{{\mathit{\rho }}_{\mathrm{a}}{C}_p}}},$$ where $u_{\ast }$ [$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$] is the friction velocity. The bulk exchange factor for heat $C_{\mathrm{bt}}=C_{\mathrm{MO}}$ is given by A7$$C_{\mathrm{MO}}={\displaystyle \frac{k_{\ast }^{\mathrm{2}}}{\left[\ln \left\{\frac{z_{\mathrm{u}}}{z_{\mathrm{0}\mathrm{m}}}\right\}-{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{u}}}{L}\right)+{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{0}\mathrm{m}}}{L}\right)\right]\left[\ln \left\{\frac{z_{\mathrm{T}}}{z_{\mathrm{0}\mathrm{h}}}\right\}-{\mathit{\psi }}_{\mathrm{h}}\left(\frac{z_{\mathrm{T}}}{L}\right)+{\mathit{\psi }}_{\mathrm{h}}\left(\frac{z_{\mathrm{0}\mathrm{h}}}{L}\right)\right]}},$$ where ${\mathit{\psi }}_{\mathrm{m}}$ and ${\mathit{\psi }}_{\mathrm{h}}$ are the vertically integrated stability functions for momentum and heat, respectively. The latent flux $Q_{\mathrm{LE}}$ is calculated analogously to scalar roughness length for water vapour $z_{\mathrm{0}q}$ and the integrated stability function for water vapour ${\mathit{\psi }}_q$, both assumed equal to the corresponding quantity and function for heat. Since $L$ depends itself on $Q_{\mathrm{H}}$, $C_{\mathrm{MO}}$ is calculated iteratively.

Different formulations are available for the integrated stability function $\mathit{\psi }$ as used by, for example, , , and . We use the following momentum ${\mathit{\psi }}_{\mathrm{m}}$ and heat ${\mathit{\psi }}_{\mathrm{h}}$ stability functions (again assuming ${\mathit{\psi }}_{\mathrm{h}}={\mathit{\psi }}_q$) : A8$$\begin{array}{rl}& {\mathit{\psi }}_{\mathrm{m}}\left(\mathit{\zeta }\right)\\ & =\left\{\begin{array}{ll}-\left[a\mathit{\zeta }+b\left(\mathit{\zeta }-\frac{c}{d}\right)\exp \left\{-d\mathit{\zeta }\right\}+\frac{bc}{d}\right],& \mathit{\zeta }>\mathrm{1}\\ & \text{(stable (Beljaars and Holtslag, 1991))}\\ \ln \left\{\left(\frac{\mathrm{1}+{\mathit{\chi }}^{\mathrm{2}}}{\mathrm{2}}\right){\left(\frac{\mathrm{1}+\mathit{\chi }}{\mathrm{2}}\right)}^{\mathrm{2}}\right\}\\ -\mathrm{2}\arctan \mathit{\chi }+\frac{\mathit{\pi }}{\mathrm{2}},& \mathit{\zeta }<\mathrm{0}\\ & \text{(unstable (Dyer, 1974))}\end{array}\right.\end{array}$$ and the heat stability function ${\mathit{\psi }}_{\mathrm{h}}$ A9$$\begin{array}{rl}& {\mathit{\psi }}_{\mathrm{h}}\left(\mathit{\zeta }\right)\\ & =\left\{\begin{array}{l}-\left[{\left(\mathrm{1}+\frac{\mathrm{2}a}{\mathrm{3}}\mathit{\zeta }\right)}^{\mathrm{1.5}}\right.\\ \left.+b\left(\mathit{\zeta }-\frac{c}{d}\right)\exp \left\{-d\mathit{\zeta }\right\}+\frac{bc}{d}-\mathrm{1}\right],& \mathit{\zeta }>\mathrm{1}\\ & \text{(stable (Beljaars and Holtslag, 1991))}\\ \mathrm{2}\ln \left(\frac{\mathrm{1}+{\mathit{\chi }}^{\mathrm{2}}}{\mathrm{2}}\right),& \mathit{\zeta }<\mathrm{0}\\ & \text{(unstable (Dyer, 1974))}\end{array}\right.\end{array}$$ with the dimensionless Obukhov parameter $\mathit{\zeta }:=z_{\mathrm{u}}/L$ (for ${\mathit{\psi }}_{\mathrm{m}}$) or $z_{\mathrm{T}}/L$ (for ${\mathit{\psi }}_{\mathrm{h}}$), $\mathit{\chi }=(\mathrm{1}-\mathrm{16}\mathit{\zeta }{)}^{\mathrm{1}/\mathrm{4}}$, $a=\mathrm{1}$, $b=\mathrm{2}/\mathrm{3}$, $c=\mathrm{5}$, and $d=\mathrm{0.35}$.

Appendix B Ground thermal and hygric regimes

The ground thermal and hygric regimes are shown in Figs.  and . These are the meteorological input variables for the calculation of the turbulent fluxes.

Figure B1

Thermal regime. (a) Radiometric ground surface temperature $T_{\mathrm{s}}$ (“rGST”) and air temperatures in the atmosphere $T_{\mathrm{a}}$ and in the instrumented cavity (cavity roof $T_{\mathrm{sa}}$ at 0.67 m and mid-cavity level at 1.95 m depth). (b) Vertical temperature gradients ${\mathrm{\nabla }}_{z}T$ between surface and 2 m air temperature (drives turbulent sensible fluxes; Eq. ), across the snow cover ($T_{\mathrm{a}}-T_{\mathrm{sa}}$), and within the cavity (indicates the stability of the in-cavity air column). The summer 2022 dry spells are referred to in the text.

[Figure omitted. See PDF]

Figure B2

Hygric regime. (a) Relative and specific humidity of the saturated snow surface $q_{\mathrm{ss}}^{\ast }$, in the atmosphere $q_{\mathrm{a}}$, and in the instrumented cavity ($q_{\mathrm{sa}}$ cavity roof at 0.67 m and mid-cavity level at 1.95 m depth). (b) Vertical specific humidity gradients (${\mathrm{\nabla }}_{z}q$) between the atmospheric air and the surface (either $q_{\mathrm{ss}}^{\ast }$ or $q_{\mathrm{sa}}$ according to Eq. ), across the snow cover ($q_{\mathrm{a}}-q_{\mathrm{sa}}$), and within the cavity (indicates the moisture transport direction). The moisture transport within the cavity is generally downwards (${\mathrm{\nabla }}_{z}q>\mathrm{0}$), and condensation rather than evaporation occurs in the deep AL. Severe dry spells that last long enough to exhaust the near-surface moisture storage ($\sim$ 5 d) strongly impact the ground moisture regime, reverse the moisture gradients, and lead to upwards moisture transport.

[Figure omitted. See PDF]

Figure B3

(a) Vertical temperature and (b) specific humidity profile (thaw-season average). Temperature and humidity are highest in the near-surface AL. The near-surface AL in contact with the atmosphere responds to dry spells and dries out within a few days after the last precipitation event. The deep AL remains close to saturation.

[Figure omitted. See PDF]

Appendix C Estimate of the thermal relaxation time

An estimate for the thermal relaxation time $\mathrm{\Delta }{t}_{\mathrm{r}}$ of a packed bed is

C1$$\mathrm{\Delta }{t}_{\mathrm{r}}\approx {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{r}}{c}_{\mathrm{r}}}{h_{\text{eff}}}}{\displaystyle \frac{V}{A_s}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{r}}{c}_{\mathrm{r}}}{h_{\text{eff}}}}{\displaystyle \frac{(d_{\mathrm{p}}/\mathrm{2})}{\mathrm{3}(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{al}})}},$$ where $h_{\text{eff}}$ is the effective heat transfer coefficient calculated via $\mathrm{1}/h_{\text{eff}}=\mathrm{1}/h+d_{\mathrm{p}}/\left(\mathrm{10}{k}_{\mathrm{r}}\right)$ . This correction accounts for the additional resistance arising from the temperature gradients within large blocks. “Large” means blocks whose Biot number exceeds $\mathrm{0.1}$ , with the Biot number defined by C2$$\mathit{Bi}:={\displaystyle \frac{h}{k_{\mathrm{r}}}}L.$$ The characteristic length $L$ is given by the ratio of the block volume to surface area: $L:=V/A_s=d_{\mathrm{p}}/\mathrm{6}$ for spheres. A minimum value for the convective heat transfer coefficient $h$ under forced convection at low wind speed ($u<\mathrm{1}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) is $\sim$ 10 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$. The value is derived from inverting Eq. () with $h:=Q_{\mathrm{H}}/(T_{\mathrm{a}}-T_{\mathrm{s}})$ (nighttime values) and agrees with an estimate in . With a thermal conductivity of the rock of $k_{\mathrm{r}}=\mathrm{2.5}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, the critical diameter is 0.15 m. With $d_{\mathrm{p}}=\mathrm{60}$ cm (typical dimension of the blocks enclosing the instrumented cavity), Eq. () yields an estimate of $\mathrm{\Delta }{t}_{\mathrm{r}}=\sim \mathrm{12}$ h. Thermal adjustment (within $>\mathrm{95}\mathit{\%}$) is attained after $\mathrm{3}\mathrm{\Delta }{t}_{\mathrm{r}}$ or $\sim \mathrm{1.5}$ d for large blocks with 60 cm edge length. Smaller blocks reach thermal equilibrium much faster, e.g. 20 cm blocks within 10 h.

Appendix D 1997–2000 wind speed profiles

Here, we justify our “katabatic correction” of wind speed measurements with wind profile data collected by and . A 10 m tower installed on the Murtèl rock glacier at the PERMOS monitoring site measured the wind velocity, air temperature, and relative humidity at 1.5, 2.0, 6.5, and 9.1 m a.g.l. between January 1997 and March 2000.

Average vertical wind speed profiles (Fig. a) show the winter-time low-level jet between December and March with maximum wind speeds close to the snow-covered surface and upwards decreasing wind speed. In summer, wind speeds increase with height and approximately show the log wind profile . The wind speed gradient is predominantly negative for more than $\sim$ 75 cm of snow measured at the PERMOS station (located on a plateau) (Fig. b), corresponding to $\sim$ 50 cm at the PERMA-XT station (located on a more wind-swept ridge). This finding justifies the compensation for snow height (Eq. ) and possibly explains why the “katabatic correction” has less effect on the SEB of the snow-poor winter 2021–2022 and a negative effect on the November 2020 SEB (Fig. ).

Figure D1

(a) Wind speed profiles for winter (December–March average) and summer (July–September) 1997–2000 . Winter storms where wind speed at 9.1 m exceeded 2.5 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ are filtered out (threshold from ). (b) The average wind speed gradient between the 2.0 and 6.5 m levels switches from a (dominantly) positive log wind profile to a negative low-level jet profile at a (PERMOS) snow height of $\sim$ 60 cm ($R^{\mathrm{2}}=\mathrm{0.271}$). Own figure based on data from .

[Figure omitted. See PDF]

Appendix E Nomenclature

Parameters and constants used in this study are tabulated in Table .

Table E1

Nomenclature: measurement variables, site-specific calibration parameters, dimensionless numbers, and constants.

Symbol	Unit	Name	Symbol	Unit	Name
$A$	1	Short-wave albedo	$t$	s, h	Time
$C_{\mathrm{bt}}$	1	Bulk turbulent heat and vapour	$u$	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$	Wind or airflow speed
		transfer coefficient	$u_{\ast }$	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$	Friction velocity
$C_p$	$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$	Isobaric specific heat capacity of moist air	$z_{\mathrm{u}}$, $z_{\mathrm{T}}$, $z_q$	m	Measurement height of wind speed,
$c_{\mathrm{w}}$, $c_{\mathrm{r}}$	$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$	Specific heat capacity of water, rock			air temperature, and humidity
$e$, $e^{\ast }$	Pa	Vapour pressure (at saturation)	$z_{\mathrm{0}\mathrm{m}}$, $z_{\mathrm{0}\mathrm{h}}$,	m	Roughness length for momentum, heat, and
$\mathrm{\Delta }{H}_{\mathrm{al}}^{\mathit{\theta }}$	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$	Sensible heat storage change	$z_{\mathrm{0}q}$		water vapour
$h_{\mathrm{al}}$, $h_{\mathrm{S}}$	m	Thickness of coarse blocky AL, snow cover	$z$	m	Vertical coordinate
${}^{\text{SEB}}{\widehat{h}}_{\mathrm{S}}^{\text{crit}}$	m	Thickness of decoupling snow cover
$K_{\mathrm{h}}$, $K_{\mathrm{w}}$	${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$	Eddy diffusivity for sensible and latent heat	$\mathit{\epsilon }$	1	Surface emissivity (snow, blocky surface)
$k_{\mathrm{S}}$	$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$	Thermal conductivity of the snow cover	$\mathit{\zeta }$	1	Dimensionless Obukhov parameter
$L$	m	Obukhov length	${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{S}}$, ${\mathit{\rho }}_{\mathrm{a}}$, ${\mathit{\rho }}_{\mathrm{r}}$	$\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$	Density of the snowpack, air, and rock
$\widehat{m}$	1	“Katabatic correction” exponent	${\mathit{\varphi }}_{\mathrm{al}}$	1	Porosity of coarse blocky AL
		(calibration parameter)	${\mathit{\varphi }}_{\mathrm{m}}$, ${\mathit{\varphi }}_{\mathrm{h}}$, ${\mathit{\varphi }}_q$	1	Stability functions for momentum, heat,
$P$	$\mathrm{Pa}$	Atmospheric pressure			and water vapour
$Q$, $S$, $L$	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$	Heat flux (general, short-wave radiation,	${\mathit{\psi }}_{\mathrm{m}}$, ${\mathit{\psi }}_{\mathrm{h}}$, ${\mathit{\psi }}_q$	1	Integrated stability functions
		long-wave radiation)
$q$, $q^{\ast }$	$\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}$	Specific humidity (at saturation)	Dimensionless numbers
$q_{\mathrm{a}}$	$\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}$	Air specific humidity (2 m a.g.l.)	$\mathit{Bo}$	1	Bowen ratio, $\mathit{Bo}:=Q_{\mathrm{H}}/Q_{\mathrm{LE}}$ (Eq. )
$q_{\mathrm{s}}$	$\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}$	Surface specific humidity	${\mathit{Ri}}_{\mathrm{b}}$	1	Bulk Richardson number (Eq. )
$q_{\mathrm{ss}}^{\ast }$	$\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}$	Saturated specific humidity
		of the snow surface
			Constants (value)
$q_{\mathrm{sa}}$	$\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}$	Specific humidity in the near-surface AL	$C_{\mathrm{pd}}$	$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$	Heat capacity of dry air (1005)
${\widehat{R}}_{z\mathrm{0}}$	1	Ratio of momentum and scalar roughness	$g$	$\mathrm{m}{\mathrm{s}}^{-\mathrm{2}}$	Gravitational acceleration (9.81)
		lengths (calibration parameter)	$k_{\ast }$	$\mathrm{1}$	Von Kármán constant (0.4)
$r$	${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$	Rainfall rate	$L_{\mathrm{m}}$	$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$	Latent heat of melting ($\mathrm{3.34}\times {\mathrm{10}}^{\mathrm{5}}$)
rH	%	Relative humidity	$L_{\mathrm{s}}$	$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$	Latent heat of sublimation ($\mathrm{2.83}\times {\mathrm{10}}^{\mathrm{6}}$)
SWE	$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$	Snow water equivalent	$L_{\mathrm{v}}$	$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$	Latent heat of vaporisation ($\mathrm{2.476}\times {\mathrm{10}}^{\mathrm{6}}$)
$T_{\mathrm{s}}$	K, °C	Surface temperature	$\mathit{\sigma }$	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{4}}$	Stefan–Boltzmann constant
		(coarse blocky AL, snow surface)			($\mathrm{5.670}\times {\mathrm{10}}^{-\mathrm{8}}$)
$T_{\mathrm{b}}$	K, °C	Temperature at base of snow cover
$T_{\mathrm{a}}$, $T_{\mathrm{wb}}$, $T_{\mathrm{v}}$	K, °C	Air temperature (dry-bulb, wet-bulb, virtual)
$T_{\mathrm{r}}$	K,°C	Temperature in blocks

Data availability

The PERMOS data can be obtained from the PERMOS network (https://doi.org/10.13093/permos-meteo-2021-01, ) and the PERMA-XT measurement data from https://www.permos.ch/doi/permos-spec-2023-1 (last access: 5 April 2024; DOI: https://doi.org/10.13093/permos-spec-2023-01, ).

Author contributions

DA performed the fieldwork, model development, and analyses for the study and wrote the manuscript. MS, MH, and BK supervised the study, provided financial and field support, and contributed to the manuscript preparation. AH and CK provided logistical support and editorial suggestions on the manuscript. HG designed the novel sensor array, regularly checked data quality, contributed to the analyses, and provided editorial suggestions on the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

This work is a collaboration between the University of Fribourg and GEOTEST. The authors wish to thank Walter Jäger (Waljag GmbH, Malans) and Thomas Sarbach (Sarbach Mechanik, St. Niklaus) for the technical support and the Corvatsch cable car company for logistical support. This publication is dedicated to Martin Scherler who laid the conceptual foundation.

Financial support

This research has been supported by Innosuisse – Schweizerische Agentur für Innovationsförderung (grant no. 36242.1 IP-EE).

Review statement

This paper was edited by Emily Collier and reviewed by two anonymous referees.