Surface energy balance closure over melting snow and ice from in situ measurements on the Greenland ice sheet

Abstract

Accurately quantifying all the components of the surface energy balance (SEB) is a prerequisite for the reliable estimation of surface melt and the surface mass balance over ice and snow. This study quantifies the SEB closure by comparing the energy available for surface melt, determined from continuous measurements of radiative fluxes and turbulent heat fluxes, to the surface ablation measured on the Greenland ice sheet between 2003 and 2023. We find that the measured daily energy available for surface melt exceeds the observed surface melt by on average 18 ± 30 W m⁻² for snow and 12 ± 54 W m⁻² for ice conditions (mean ± SD), which corresponds to 46 and 10% of the average energy available for surface melt, respectively. When the surface is not melting, the daily SEB is on average closed within 5 W m⁻². Based on the inter-comparison of different ablation sensors and radiometers installed on different stations, and on the evaluation of modelled turbulent heat fluxes, we conclude that measurement uncertainties prevent a better daily to sub-daily SEB closure. These results highlight the need and challenges in obtaining accurate long-term in situ SEB observations for the proper evaluation of climate models and for the validation of remote sensing products.

Introduction

The surface energy balance (SEB) is a fundamental physical principle stating that the sum of all the energy fluxes at the interface between the atmosphere and the subsurface must be equal to zero. The SEB is widely used to couple atmospheric models to ocean or land surface models, and over snow and ice surfaces to compute the surface melt rate when the surface is at the melting point. Hence, the assumption of a closed SEB is important for simulations of the Earth system, and in situ SEB observations are an essential part of climate model development. The SEB for an infinitesimally thin, horizontally homogeneous surface skin layer between the air and the subsurface is written as

(1)$$Q_{\rm SW} + Q_{\rm LW} + Q_{\rm H} + Q_{\rm E} + Q_{\rm G} + Q_{\rm R} = Q_{\rm M} + {\rm RES},\; $$

where Q _SW and Q _LW denote the net absorbed shortwave and longwave radiation, Q _H the turbulent sensible heat flux, Q _E the turbulent latent heat flux, Q _G the subsurface conductive heat flux, Q _R the added energy by rain and Q _M the energy used for surface melt and is defined positive. RES is the residual flux and must therefore be equal to zero if all the fluxes are included and the SEB is closed. All the fluxes are evaluated at the surface, defined positive when directed towards the surface, and expressed for an averaging time interval Δt (s) as energy flux density (W m⁻²).

Previous analyses of in situ SEB observations for various land surfaces found that the sum of measured turbulent heat fluxes is often smaller than the measured net available radiation at the surface (Wilson and others, 2002). This is commonly referred to as the ‘surface energy balance closure problem’ (Foken, 2008; Mauder and others, 2020). The apparent non-closure of the SEB over land surfaces is found to be partly explained by horizontal advection caused over a non-homogeneous surface cover which is not recorded by sensors (Stoy and others, 2013; De Roo and Mauder, 2018). Also, the subsurface fluxes are often not negligible at (sub)daily timescales, but may still be zero on average for longer timescales (Leuning and others, 2012; Grachev and others, 2020). Furthermore, instruments can underestimate the turbulent heat fluxes during stable and/or low-wind situations (De Roode and others, 2010; Helgason and Pomeroy, 2012), or when the averaging interval is not long enough to capture all the fluctuations contributing to the turbulent fluxes (Foken and others, 2006). Finally, the measured radiation over a sloping or heterogeneous surface may not be representative of the footprint of the measured turbulent heat fluxes (Georg and others, 2016). In any case, an accurate quantification of all the fluxes making up the SEB is challenging but crucial to provide a reliable benchmark dataset for model evaluation and development.

Sensors installed on the slopes of an ice sheet typically have the advantage of having a large surrounding area of a fairly homogeneous surface cover, while the persistent katabatic wind prevents very stable situations to occur. Furthermore, all the terms in the SEB related to the vegetation canopy (photosynthesis, storage) are zero, and the subsurface heat conduction is small in the ablation zone when the subsurface becomes isothermal during the melting season. Finally, the ablation of the snow/ice surface due to melt and sublimation can be directly measured and converted to an energy flux if the subsurface density is known:

(2)$${\Delta h\over \Delta t} = \left({Q_{\rm M}\over \rho_{\rm i, s}L_{\rm m}} - {Q_{\rm E}\over \rho_{\rm i, s} L_{\rm s}} \right)$$

where Δh is the measured surface lowering during a time interval Δt, and is defined positive when there is ablation, either due to melt (Q _M > 0) or due to sublimation (Q _E < 0). ρ _i,s is the snow or ice density, L _m = 3.34 × 10⁵ J kg⁻¹ the latent heat of fusion and L _s = 2.8345 × 10⁶ J kg⁻¹ the latent heat of sublimation. Hence, the measured surface lowering can be used to estimate the energy available for surface melt:

(3)$$Q_{\rm M} = \rho_{\rm i, s}L_{\rm m} {\Delta h_{\rm m}\over \Delta t},\; $$

where Δh _m is defined as the surface lowering caused by melt only, and may be estimated from the observed lowering, corrected for sublimation (Eqn (2)):

(4)$${\Delta h_{\rm m}\over \Delta t} = {\Delta h\over \Delta t} + {Q_{\rm E}\over \rho_{\rm i, s} L_{\rm s}}.$$

Besides the logistical difficulties in operating and servicing equipment for extended periods of time on glaciers and ice sheets, five main challenges of measuring the SEB of a snow/ice surface are related to: (1) the fraction of shortwave radiation that penetrates into the subsurface, (2) the evolving surrounding topography and albedo due to for example ice hummocks, crevasses, patches of snow, meltwater streams and cryoconite holes, (3) the presence of a wind maximum close to the surface due to the katabatic forcing, which potentially impacts the turbulent flux observations, (4) the fact that subsurface temperature sensors melt out during the melting season which complicates the estimation of the subsurface heat flux and (5) the partitioning of Q _E between ice sublimation and liquid water evaporation over a melting surface. Despite these challenges, both Fitzpatrick and others (2017) and Litt and others (2017) found that the measured surface lowering during an ablation season on a mountain glacier matched the cumulative sum of the observed energy fluxes, yet an RMSE of 1 cm ice ablation per day persisted in the daily SEB closure (Fitzpatrick and others, 2017), which corresponds to an average energy flux of 35 W m⁻² (Eqn (3)). No other studies have compared all the daily individual SEB fluxes over a melting ice and snow surface on an ice sheet, either due the absence of eddy covariance systems to measure the turbulent heat fluxes, or due to the absence of continuous and accurate methods to measure surface lowering. On the other hand, many studies rely on the SEB using modelled turbulent heat fluxes to quantify melt and sublimation rates using automatic weather station (AWS) data, not only on the Greenland ice sheet (e.g. Charalampidis and others, 2015; Kuipers Munneke and others, 2018; Vandecrux and others, 2020; Covi and others, 2022), but also on the Antarctic ice sheet (e.g. Buzzard and others, 2018; Jakobs and others, 2020) and mountain glaciers (e.g. Stigter and others, 2018; Brun and others, 2023).

In this study we investigate how well the SEB closes based on in situ measurements performed over melting ice and snow. We first consider five locations on the Greenland ice sheet where eddy covariance observations are available, and we then extend the analysis to 19 locations with different observed surface mass balances (SMBs). We consider various timescales, ranging from the sub-daily timescales to the entire melting season. Additionally we study different instrumental errors and discuss the possibility of missing terms in the SEB.

Measurements and post-processing

Experiments

Eddy covariance measurements were performed during selected periods at S5, S6, KAN_L, S10 and QAS_L (Table 1), which are all located in the ablation zone of the Greenland ice sheet except for S10 (Fig. 1). The combined dataset contains measurements from a CSAT3 or CSAT3B (Campbell Scientific) and LI-7500 (LI-COR) at S6 during 2003–04 (Smeets and van den Broeke, 2008), at S10 during 2012 (Lenaerts and others, 2014) and at S5 during 2005–06, 2008 and 2019–23, at QAS_L during 2019–22 and at KAN_L during 2019–22 (van Tiggelen and others, 2023). We also include AWS data from locations without eddy covariance measurements, either from the Institute for Marine and Atmospheric Research Utrecht (IMAU, Smeets and others, 2018; Kuipers Munneke and others, 2018) or from the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) AWS network (How and others, 2022, version 12). At each site, an AWS continuously recorded the four radiation components using a CNR1 or CNR4 radiometer (Kipp & Zonen), wind speed and direction, air temperature, relative humidity, air pressure and surface height change at hourly frequency (Smeets and others, 2018; Fausto and others, 2021). The longwave radiation measurements are corrected for window heating after Smeets and others (2018), and the shortwave radiation measurements are corrected for a zero-bias (Behrens, 2021) and for tilt after van den Broeke and others (2004). The considered locations and amount of available data are given in Table 1.

Table 1. Description of AWS and sonic eddy covariance (SEC) data used in this study

The coordinates of S21, S22 and S23 are at installation date. For the other stations the coordinates are from Smeets and others (2018) and Fausto and others (2021). $H_{\max }$ = maximum height of the ice hummocks used in the bulk turbulence model for the calculation of the surface aerodynamic roughness length. n.a. = not applicable since the surface is firn-covered year round.

^aTwo separate draw wire measurements are available since September 2021.

^bThese stations include in situ Q _E measurements.

^cS21 is sometimes also referred to as S20 in published literature.

^dS21 was moved 23 km northwest on 10 August 2015.

^eS22 and S23 moved by on average 0.7 and 1.7 km respectively due to ice flow.

Figure 1. Location of the AWS on the Greenland ice sheet from the IMAU and PROMICE networks used in this study. Panels contain insets for each transect. Circles denote the locations of the AWS, and triangles denote AWS locations with SEC observations. Background colour denotes the average modelled downscaled yearly SMB in the period 1958–2019 from the regional climate model RACMO2.3p2 (Noël and others, 2016). The two consecutive locations for S21 are shown.

The processing steps and corrections used to obtain Q _H and Q _E are described in van Tiggelen and others (2020). The corrections include de-spiking of the raw data using a moving median filter, 30 min averaging and rotation of the average wind vector in the mean flow direction (Kaimal and Finnigan, 1994), corrections for frequency losses (Moore, 1986), corrections of Q _H for humidity influences (Schotanus and others, 1983) and corrections of Q _E for dilatation and dilution (Webb and others, 1980). The variable instrument height due to snow accumulation and melt is taken into account using measured height changes by a sonic ranger mounted on the AWS boom. The hourly Q _H and Q _E measurements are only retained when all the following conditions are valid in order to minimise the influence of noise and non-stationary conditions: (1) air temperature does not change by more than 0.6°C h⁻¹, (2) wind speed above 1 m s⁻¹, (3) friction velocity above 0.1 m s⁻¹, (4) Q _H between −100 and 500 W m⁻², (5) Q _E between −500 and 500 W m⁻² and (6) at least 5% of valid 10 Hz data per 30 min. The latter selection removes between 15 and 25% of data at locations S5, KAN_L and S6, 33% of data at S10 and 56% of data at QAS_L.

Finally, surface lowering was measured using one of the following three methods: a draw wire (DW, Hulth, 2010), either mounted on the AWS or on a separate tripod, a sonic height ranger (SR) either mounted on the AWS boom or on a separate stake and/or a pressure transducer assembly mounted on the AWS (PTA, Fausto and others, 2012). The available type of ablation measurements per station is given in Table 1.

Model for subsurface and turbulent heat fluxes

When no direct Q _H or Q _E observations are available, which is 96 and 99% of the time when considering all the AWS station-years respectively due to limited eddy covariance data, the bulk approach is used to compute turbulent fluxes from single-level wind, air temperature and humidity observations assuming similarity theory, as described in van Tiggelen and others (2023). These data are then only used for the long-term estimation of SEB closure at all the considered AWS. The surface temperature is estimated from the measured outgoing longwave radiation, assuming an emissivity of 1, which is on the high-end of the typical values of 0.95–1.00 for snow and bare ice (Hori and others, 2006). The window heating offsets from the pyrgeometer prevent a more accurate determination of the surface emissivity and surface temperature. The bottom temperature is fixed to the multi-annual average air temperature. The roughness lengths for momentum and heat are time-varying and modelled using melt and snow depth after van Tiggelen and others (2023). The height of the ice hummocks is reset to $H_{\max }$ every year on the 1st of September. Values of $H_{\max }$ are chosen based on photographs taken during the end-of-melt season maintenance visits (Table 1).

The subsurface heat flux is computed using a vertical heat diffusion model in which the subsurface is discretised from the surface down to 25 m depth with layers increasing in thickness, ranging from 0.01 m thickness at the surface up to 2 m thickness at the bottom, and a constant time step of 10 s. The forcing is linearly interpolated between two consecutive hourly measurements. The snow density is fixed at 350 kg m⁻³, which is most likely too low for melting snow conditions (e.g. Gugerli and others, 2019) but is chosen after Fausto and others (2018) due to the lack of continuous measurements of snow density. This assumption may lead to underestimated melt fluxes in the presence of heavier snow layers, refrozen layers or slush. Ice density is set to 916 kg m⁻³. The snow/ice heat conductivity is parametrised as a function of density after Calonne and others (2019), and liquid water percolation in the snow is computed using the bucket method described in Schneider and Jansson (2004). The addition and removal of layers due to accumulation and ablation is prescribed using the AWS measurements.

The penetration of shortwave radiation below the surface is not considered due to the lack of in situ observations. Instead we assume that all the net shortwave radiation is absorbed by the skin surface layer and is used for melt when the surface is at the melting point. This assumption may lead to overestimated melt fluxes in cases when the absorbed shortwave radiation is used to warm the snow layers instead (e.g. Kuipers Munneke and others, 2009). This effect is however deemed limited for the SEB quantification, due to the corresponding reduction of the conductive heat flux in isothermal melting layers, and due to the fact that the absorbed subsurface radiation and the conductive heat flux are much smaller than the other SEB fluxes (van den Broeke and others, 2008).

Analysis structure

In the next part, the in situ observations at five locations with eddy covariance observations are used to quantify the SEB fluxes Q _SW, Q _LW, Q _H, Q _E (when available) and the surface lowering, while a model determines Q _E and Q _G. The energy available for surface melt, Q _M, is estimated using the measured surface lowering corrected for sublimation after Eqn (2). The influence of energy advected by rain Q _R is neglected due to its on average minor contribution and lack of long-term measurements, although it can become an important source of energy during extreme precipitation events (Box and others, 2023). The sum of all the SEB fluxes is used to estimate the SEB residual (RES, Eqn (1)). We consider the intra-daily, daily and seasonal timescales in order to identify the different SEB residuals and limitations in the observations.

In the final part, the long-term SEB closure and the uncertainty in melt energy using the SEB equation (Eqn (1)) is quantified at all the 19 locations with modelled turbulent heat fluxes, which are more common since eddy covariance observations are very scarce on ice sheets. The analysis is extended with up to 19 years of single-level observations at all locations shown in Figure 1. Seven stations are operated by the IMAU (Smeets and others, 2022), while the remainder are part of the PROMICE network (How and others, 2022).

Results

Seasonal SEB closure

Time series of daily averaged observed SEB fluxes, and modelled and observed lowering at QAS_L and S5 are shown in Figure 2. At site QAS_L, Q _H measurements are available for autumn and winter between 2019 and 2021, and for the entire melt season during 2022. During the ablation season, the range of measured surface lowering using the three different ablation sensors envelop the estimated lowering by closing the SEB in 2020 and 2021 (Figs 2b, c), while it is larger than the SEB sum in 2019 (Fig. 2a) and lower than the SEB sum in 2022 (Fig. 2e). During the snowmelt period in spring 2022 (Fig. 2d), the lowering measured by two sonic height rangers differ by 0.19 m (27% of the 0.7 m total lowering), yet the SEB modelled lowering, corrected for accumulation, is substantially larger: 1.77 m, that is, by a factor 2.5. This result is consistent with the result from Covi and others (2022), who also identified a persistent positive residual in the SEB fluxes during a snowmelt period. However when considering the ice ablation season of 2022 (Fig. 2e), the cumulative ice ablation measured by the pressure transducer assembly (5.27 m) and the sonic ranger (4.96 m) agree within 4% of the average of both (5.12 m). Albeit larger in absolute sense, the SEB residuals are still similar compared to previous years, which was 9% of the total energy available for surface melt in 2022. The draw wire measured 0.7 m less ablation (14%), of which about half is explained by the estimated 70–100 cm lateral displacement of the draw wire with respect to the borehole (Hulth, 2010), since the wire was attached 66 cm above the ice surface at QAS_L. The estimated lowering by summing the measured SEB fluxes yields a total ablation of 5.33 m, which is in close agreement with the pressure transducer assembly measurement. This means the SEB is closed within 6 cm ice ablation (1%) over the entire ablation season at QAS_L, and suggests that spatial variability and residuals in the SEB fluxes are of limited influence.

Figure 2. Measured surface lowering and daily SEB residual at locations QAS_L (a–e) and S5 (f–h). Each panel contains a different period during melting conditions when multiple surface lowering measurements and Q _H measurements are available. Panel (d) is a period with melting snow, while the other panels are periods with melting ice conditions. DW denotes the draw wire mounted on a separate tripod, and DW AWS the draw wire mounted on the AWS. The surface lowering Δh is computed using the SEB fluxes after Eqn (2).

At S5 (Figs 2f, g,h), Q _H measurements are available for three consecutive melt seasons from 2021 to 2023. On average, the measured cumulative lowering by the two separate draw wires differ by 0.45 m (9% of the average lowering of 5 m), with the draw wire mounted on the AWS mast recording smaller values in the first 2 years. This difference can for a large part be explained by the lateral movement of the draw wires, but also by the spatial variability in surface ablation that differ between each ablation season. During 2023, a crevasse opened under the AWS, which caused a sudden lowering of the station at the end of the melt season, thereby explaining the higher draw wire on the AWS values during 2023. During 2023, measurements from an additional sonic height ranger on a stake are available, which recorded 0.2 m less total ablation than the draw wire mounted on a separate tripod (5% of the average 4.1 m lowering). On average, the estimated lowering by closing the SEB is 0.6 m larger (≈20% of the measured ablation) than the draw wire, with large inter-annual variability, meaning that there is a positive SEB residual.

Daily SEB closure

Here we analyse the daily averaged SEB fluxes in order to quantify the daily SEB residual for snow and ice under different meteorological and surface conditions. The measurements from the five stations with eddy covariance observations are combined in three groups: (1) days with snow at the surface and without melt, defined as days when the albedo is larger than 0.6, when no ice ablation is measured and when the observed daily maximum surface temperature does not exceed −2$^\circ$C, (2) days with bare ice and melt, defined as days when ice ablation is measured by the ablatometer and when the observed daily maximum surface temperature exceeds −2$^\circ$C and (3) days at locations S6, S10 and QAS_L with melting snow, defined as days when the sonic height ranger mounted on the AWS boom records a surface lowering, when there is no measured ice ablation and when the observed daily maximum surface temperature also exceeds −2$^\circ$C. This selection results in a total of 1926  days with simultaneous AWS, Q _H and ablation measurements. The measured SEB fluxes for each case are shown in Figure 3. For all three cases, the average SEB residual is positive, and increases from 5 W m⁻² for non-melting conditions to 12 W m⁻² for melting ice conditions and to 18 W m⁻² for melting snow conditions. The SD of the residual also increases from 27 W m⁻² for non-melting conditions to 54 W m⁻² for melting ice conditions. All in all, the residual is small compared to the absolute sum of all the SEB fluxes. Although the SEB residuals for melting snow and ice conditions are similar in absolute sense, the residual is a much larger fraction of the total energy available for melt during melting snow conditions: 46% compared to 10% during melting ice. Hence, the SEB is nearly closed for melting ice, but not for melting snow conditions. The average residual of 18 W m⁻² translates to 1.3 cm d⁻¹ of snow ablation (Eqn (2) with ρ _s = 350 kg m⁻²), explaining the largely overestimated snow surface lowering at QAS_L during spring 2022 shown in Figure 2.

Figure 3. Daily averaged SEB fluxes at all locations with eddy covariance observations for non-melting conditions (a), melting snow conditions when no accumulation is observed at QAS_L, S6 and S10 (b) and for melting ice conditions at all locations with eddy covariance observations (c). The average, SD and number of SEB residual observations is shown in each panel. The fraction of the average SEB residual compared to the average energy available for melt is shown as percentage for melting conditions (b, c). Colours denote the data density. Diagonal lines represent perfect SEB closure. Selection criteria are defined in the main text.

When no melt occurs, the net radiative fluxes are smaller than the non-radiative fluxes by on average 5 ± 27 W m⁻². A possible explanation is a measurement uncertainty in the radiometer, or the omission of the penetration of shortwave radiation below the surface in the calculations. The latter may cause a reversed subsurface heat flux through warming below the surface. When selecting only days during nighttime or during the polar night when daily downward shortwave radiation is smaller than 10 W m⁻², the SEB residual persists (Fig. 4a) but remains small (4 W m⁻²), which also indicates a small measurement uncertainty in the radiometer and/or compensating errors. A different explanation may be an underestimated turbulent heat exchange during low wind conditions by the sonic anemometer, but selecting days with daily average wind speeds above 5 m s⁻¹ also yields a small but positive SEB residual (Fig. 4b). A third explanation may be the inaccurate modelling of the subsurface heat fluxes, but considering days with modelled subsurface heat fluxes smaller than 3 W m⁻² does not improve the SEB residuals either (Fig. 4c). A last explored potential explanation of the SEB residual is a combination of measurement errors. A typical uncorrelated random error of 10 W m⁻² for all the five considered SEB fluxes (Eqn (1)) would result in an expected SD of the residual of 22 W m⁻². Since the estimated SD of the observed residual for non-melting conditions ranges between 23 and 30 W m⁻², depending on the case selection (Fig. 4), we conclude that the observed residual could be mostly explained by the sum of random measurement errors.

Figure 4. Same as Figure 3 but only for non-melting conditions and for additional data selection: days when daily Q _SW < 10 W m⁻² (a), days when daily average wind speed is higher than 5 m s⁻¹ (b) and days when daily averaged |Q _G| < 3 W m⁻² (c).

Intra-daily SEB closure

Measurements during the ice ablation season from QAS_L (31 May–27 August 2022, Fig. 2e) and from S5 (2 June–20 August 2023, Fig. 2h) are averaged per time-of-day to obtain the average SEB daily cycle for bare ice conditions only (Fig. 5). At both sites, the SEB budget is dominated by the daily cycle in Q _SW, which results in a daily cycle in observed ablation. At QAS_L (Fig. 5a), the draw wire and the sonic ranger record such a daily cycle, and as seen before, somewhat less ablation compared to the SEB sum. Differences between these two sensors are due to the anchoring of the draw wire tripod several centimetres below the surface, meaning that the Q _SW-driven daytime melt is most likely overestimated when the tripod is sinking into the ice, while the Q _H-driven nighttime melt is most likely underestimated as long as it occurs above the anchoring depth. On the other hand, the pressure transducer assembly at QAS_L measured larger ablation in the morning and late afternoon, while it recorded a non-physical height increase around local noon. This may be due to malfunctioning of the pressure transducer assembly expansion bladder, which allowed the pressure of the fluid to increase due to solar heating, which is then wrongly interpreted as a height increase. Nevertheless, the pressure transducer assembly and sonic ranger still yield the best daily total ablation compared to the measured SEB budget: daily difference of 12 mm (43 W m⁻²) and 5 mm (18 W m⁻²), respectively. At S5 (Fig. 5b), the two draw wires and the sonic ranger all record a daily cycle in ablation, and the differences compared to the SEB sum are most likely explained by a combination of lateral movement of the draw wires, different tripod anchoring and heterogeneous ablation. It must be noted that, as opposed to the sonic ranger measurements, the pressure transducer assembly and the draw wire are not designed to accurately capture ablation variability at sub-daily timescales due to their sensitivity to rapid changes in melt conditions, and that the measurement uncertainty of the pressure transducer assembly as determined in the field is ~4 cm (Fausto and others, 2012).

Figure 5. Hourly evolution of measured SEB fluxes and surface lowering averaged at QAS_L for the period 31 May–16 August 2022 (a) and at S5 for the period 12 June–20 August 2023 (b). Periods are chosen based on available observations during the ice ablation season. The vertical extent of the shaded area denotes twice the standard variability per hour of day.

Discussion

Uncertainty in measured ablation

The random error in daily averaged Q _M originating from the different ablation sensors is quantified by comparing different measurements at the same location. At QAS_L, 336 days of simultaneous sonic ranger, draw wire and pressure transducer assembly measurements are available during the ablation seasons of 2019, 2020 and 2022. At S5 and S6, 282 and 60 days of simultaneous measurements are available from two separate draw wires, respectively, in the ablation seasons between 2021 and 2023. At QAS_L (Figs 6a, b) the draw wire measured on average 5 W m⁻² less daily melt energy than the pressure transducer assembly, and 19 W m⁻² less daily melt energy than the sonic ranger on a stake. The average difference is notably smaller than the RMS between all three sensors, which ranges between 89 and 109 W m⁻², corresponding to ~2.7 cm ice ablation per day. This is in agreement with previous studies that quantified the random errors from the draw wire and pressure transducer assemblies in the field (Hulth, 2010; Fausto and others, 2012). At S5 and S6 (Fig. 6c), the two draw wires differ on average by 14 and 2 W m⁻², respectively, with an RMS of 45 and 22 W m⁻². Location S5 being in very hummocky terrain, the deviations are most likely due to both lateral displacement of the sensors and the spatial variability in ablation. At S6, the terrain is much flatter, hence the smaller RMS of 22 W m⁻² corresponding to 6 mm of daily ice ablation can be interpreted as the random instrumental error for the draw wire. To summarise, the random error in daily averaged Q _M measured by the draw wire, the pressure transducer assembly and the sonic ranger is of the order of 10 W m⁻², corresponding to ~3 mm of daily ice ablation.

Figure 6. Inter-comparison of measured daily energy used for melt at QAS_L (a, b) and S5 and S6 (c) from different ablation sensors. In panels (a, b) the draw wire (DW), sonic ranger (SR) on a stake and pressure transducer assembly (PTA) are compared. In panel (c) two adjacent DW sensors are compared. The average difference and RMSE is given in each panel.

Uncertainty of modelled turbulent heat fluxes

The modelled, daily averaged Q _H and Q _E using single-level observations are compared against the available eddy covariance observations in Figure 7. At S5 and KAN_L, the modelled Q _H is on average overestimated by 10 and 14 W m⁻², respectively, while at QAS_L and S6, Q _H is accurately modelled within 5 W m⁻² and Q _E within 10 W m⁻². At S10, the modelled Q _H and Q _E are both underestimated by 32 and 9 W m⁻², although less data are available for this location. At S5 and KAN_L, the high and variable values for the roughness lengths, as well as the variable height of the surface in the footprint of the sensors (zero-referencing) are the major sources of uncertainty for modelling the turbulent heat fluxes. This is expected to be less problematic at the less hummocky locations QAS_L and S6, where the Q _H bias is indeed smaller, but where part of the measured negative Q _E values are not well represented by the model. These are days with high relative humidity (>80%) and near-zero measured net longwave radiation (Q _LW > −20 W m⁻²), suggesting either condensation/riming on the sensors or low-level clouds/fog that negatively impact the measurements. At S10, the major sources of uncertainty are the blowing snow events that influence the radiation measurements and the vertical profiles, as well larger sensitivity to measurement errors due to the smaller vertical temperature and humidity gradients typically found at these higher locations in the accumulation zone. However, longer-term measurements are required to properly quantify these effects.

Figure 7. Evaluation of modelled daily averaged Q _H (red dots) at location S5 and KAN_L (a, b) and both daily averaged Q _H (red dots) and Q _E (blue dots) at location QAS_L (c), S6 (d) and S10 (e). The value denotes the bias and the RMSE between observed and modelled daily fluxes.

Uncertainty of measured radiative fluxes

The uncertainty in measured radiative fluxes is determined by comparing the measurements from S10 and KAN_U, since both AWSs were situated at the same location between August 2010 and April 2016. After the correction for the zero bias and for the influence of tilt, the daily total shortwave measurements differ on average by 7% for the downwards, and by only 1% for the upward component (Figs S2 and S3), which is less than the 10% uncertainty given by the manufacturer. This results in a daily average difference in Q _SW of 8 W m⁻². The Q _LW observations differ by only 2 ± 6 W m⁻² SD. At the lower locations in the ablation zone, the uncertainty may be quantified by comparing 15 years of measurements from station pairs S5–KAN_L, and S6–KAN_M, due to their relative proximity and similar negative SMB. This long-term comparison yields an average daily Q _LW difference of 2 ± 7 W m⁻² (6%) and 2 ± 10 W m⁻² (6%), and a daily Q _SW difference of 7 ± 12 W m⁻² (14%) and 6 ± 53 W m⁻² (13%) between S5–KAN_L and S6–KAN_M, respectively. The differences at all the sites are possibly explained by the shading of the radiometer by the mast or solar panels and the imperfect correction for tilt, but also by spatial differences, for example, the presence of algae and impurities at the surface or by the accumulation of snow or ice on top of the sensors, which were not quantified in this study.

Long-term SEB closure from single-level observations

The 10-daily averaged SEB residuals for the seven stations located on the K-transect (Fig. 1) during the entire period of AWS observations are shown in Figure 8, and the average and SD of the 10-daily averaged SEB residuals for all the stations for different conditions are given in Table 2. The time series of the other considered AWS are given in the Supplementary materials (Figs S4–S6). Only days with 24 valid hourly measurements are considered, which reduces the number of data by between 0.4% at S5 up to 73% at S23 (Fig. S7), or 17% on average of all the data. Extending the previous analysis to include all these AWS reveals that the positive SEB residual becomes larger at the low SMB sites QAS_L, S5, KAN_L, most likely due to overestimated Q _H and Q _E (Figs 7a–c). At the higher elevation sites of the K-transect (S6–S10), the SEB residual during non-melting conditions is on average smaller than 10 W m⁻², which means that the SEB is essentially closed within measurement uncertainty, despite the fact that the turbulent heat fluxes were underestimated during the limited period with eddy covariance data at S10 (Fig. 7e). For all locations, the average SEB residual during melting conditions is smaller than the SD, but the SD is at least two times larger for melting conditions than for non-melting conditions. Hence, the validity of the SEB cannot be challenged for melting conditions, mainly due to the random error in measured ablation mentioned above. For non-melting conditions, the average and SD exceed 10 W m⁻² for nine locations, indicating that the SEB may not be closed at these sites. Exceptionally large SEB residuals occur occasionally (Fig. 8), and are mainly explained by malfunctioning of some sensors. This is for example the case at KAN_U during the strong melting season of 2019, where Q _SW was overestimated due to large AWS tilt.

Figure 8. (a) The 10-daily averaged SEB residual for each AWS on the K-transect during the period of available AWS data. Dots are days with a non-melting surface, triangles are days with a melting ice surface only and squares are days with a melting snow surface only. (b) Boxplots containing the minimum, maximum, interquartile range and median of the 10-daily residual per station and per surface type. Averages and SDs are given in Table 2.

Table 2. Average and SD of the daily SEB residuals expressed in W m⁻² at all the considered AWS stations for non-melting conditions, for melting ice conditions using the measured ablation from the draw wire (DW), pressure transducer assembly (PTA), sonic ranger (SR) or the average of all ablation measurements, and for melting snow conditions during the entire measurement period

The surface is assumed to melt if the observed daily maximum temperature exceeds 2°C.

For these long time series, the computation rather than the direct observation of turbulent heat fluxes introduces a source of uncertainty, which is of the same order of magnitude as the observed SEB residual during periods with flux observations (Fig. 7). Furthermore, changes in AWS design and placement may also affect the SEB residual over time, possibly explaining the increase in SEB residual at S6 after 2016 when a different type of AWS was installed. Nevertheless, the average daily SEB residual considering all the AWS measurements is small: 4 ± 25 W m⁻² for non-melting conditions and 14 ± 55 W m⁻² for melting ice conditions, which is very similar to the SEB residual when only considering shorter periods at selected sites when turbulent flux observations are available (Fig. 3).

Conclusion

Given the importance of the SEB to compute surface melt over snow-covered and/or glaciated surfaces, obtaining reliable in situ measurements that close the SEB is crucial for model development, calibration and evaluation and for satellite data validation. Here we use a unique dataset of continuous in situ measurements of turbulent heat fluxes, surface ablation and radiative fluxes from five locations on the Greenland ice sheet to show that the energy available for surface melt is generally higher than the actual recorded surface ablation. For bare ice conditions this is on average 12 ± 54 W m⁻², or ~10% of the average energy available for melt. Such an energy flux generates 0.28 ± 1.26 m w.e. of additional melt per year for a typical 90 day melt season when included in a SEB model. For melting snow, the residual is similar (18 ± 30 W m⁻², ~46% of the average energy available for melt), which in combination with the lower density of snow leads to large errors in snow ablation of the order of 1.5 cm d⁻¹. Assuming a total Greenland ice sheet melt flux of 600 Gt ${\rm a}^{-1}$, a representative modelled value for the period 2006–15 (van Dalum and others, 2024), and that half of this yearly melt contributes to runoff, then a 10% uncertainty in melt energy translates to 0.08 mm ${\rm a}^{-1}$ (30 Gt ${\rm a}^{-1}$) uncertainty in the global sea-level rise contribution from the Greenland ice sheet.

During non-melting conditions, the sum of the non-radiative fluxes balances the radiative fluxes within 5 ± 27 W m⁻². This small SEB residual supports the reliability of these observations, that are often taken in remote and harsh environments with limited maintenance, but are essential for model development.

A better closure of the SEB at both daily and sub-daily timescales does not only require ablation measurements with sub-millimetre accuracy and the consideration of spatially variable SEB fluxes, but also continuous measurements of the subsurface density profile and of the amount of absorbed shortwave radiation below the surface. Furthermore, measurements of both radiative and turbulent heat fluxes with an accuracy of 5 W m⁻² or better are required. Over snow, the impact of blowing snow on the measurements needs a more detailed quantification as well.

The analysis of long-term (2003–23) single-level AWS observations at 19 locations with widely ranging SMB reveals that the multi-annual observed SEB averaged over all the locations is closed within 4 ± 25 W m⁻² during non-melting conditions, and within 14 ± 55 W m⁻² during melting conditions. Hence we conclude that the skin layer SEB approach is in general an appropriate method for quantifying the energy exchange at a snow/ice surface using single-level AWS measurements.

This quality-controlled in situ dataset of SEB fluxes may serve as a benchmark for intermediate complexity or statistical methods to estimate melt rates, or may be used as independent evaluation data for climate models or for the retrieval of surface properties by remote sensing observations.