Improved processing methods for eddy covariance measurements in calculating sensible heat fluxes at glacier surfaces

Abstract

Bulk aerodynamic methods have been shown to perform poorly in computing turbulent heat fluxes at glacier surfaces during shallow katabatic winds. Katabatic surface layers have different wind shear and flux profiles to the surface layers for which the bulk methods were developed, potentially invalidating their use in these conditions. In addition, eddy covariance-derived turbulent heat fluxes are unlikely to be representative of surface conditions when eddy covariance data are collected close to the wind speed maximum (WSM). Here we utilize two months of eddy covariance and meteorological data measured at three different heights (1 m, 2 m, and 3 m) at Kaskawulsh Glacier in the Yukon, Canada, to re-examine the performance of bulk methods relative to eddy covariance-derived fluxes under different near-surface flow regimes. We propose a new set of processing methods for one-level eddy covariance data to ensure the validity of calculated fluxes during highly variable flows and low-level wind speed maxima, which leads to improved agreement between eddy covariance-derived and modelled fluxes across all flow regimes, with the best agreement (correlation >0.9) 1 m above the surface. Contrary to previous studies, these results show that adequately processed eddy covariance data collected at or above the WSM can provide valid estimates of surface heat fluxes.

1. Introduction

Turbulent heat fluxes have been observed to be important contributors to the surface energy balance of mountain glaciers (e.g. Hock and Holmgren, 1996; Greuell and Smeets, 2001; Fitzpatrick and others, 2017). Variation in the magnitude of the turbulent heat fluxes will, therefore, have a substantial influence on surface melt rates, highlighting the need for accurate estimation of these energy terms in ablation models. One of the key uncertainty sources when it comes to modelling turbulent heat fluxes is their parameterization through bulk aerodynamic methods. The performance of these bulk methods has been evaluated in relatively few glacier studies (e.g., Hay and Fitzharris, 1988; Hock, 1998; Denby and Greuell, 2000; Conway and Cullen, 2013; Radić and others, 2017), most of which highlighted a gap in our understanding of why and when these methods fail. In order to narrow uncertainties in projections of glacier melt, it is therefore necessary to narrow uncertainties in the modelling of turbulent heat fluxes.

Due to their simplicity and reliance only on standard meteorological measurements at one height (often 2 m) above the surface, the bulk methods have been the most commonly used models for deriving turbulent heat fluxes at glacier surfaces (Guo and others, 2011; MacDougall and Flowers, 2011; Conway and Cullen, 2013; Fitzpatrick and others, 2017; Steiner and others, 2018). In their foundation, the bulk methods assume flat, homogeneous surfaces with logarithmic wind speed profiles and turbulent fluxes that are near-constant in height (varying by less than 10 $\%$) within the surface boundary layer (Stull, 1988). Since the near-logarithmic wind profiles are observed only during neutral atmospheric stability conditions (Stull, 1988), corrections are commonly applied to adjust the turbulent fluxes for non-neutral stratification. The theories and empirical data used for developing these corrections were obtained from studies over non-glacierized and flat terrain (e.g. Monin and Obukhov, 1954; Dyer, 1974; Holtslag and De Bruin, 1988; Beljaars and Holtslag, 1991), and generally assume that turbulence generation will be suppressed (enhanced) in stable (unstable) conditions.

The structure of a surface boundary layer at sloping glacier surfaces can differ greatly from that of a stable surface layer over a flat surface (van der Avoird and Duynkerke, 1999). Sloping glacier surfaces under stable conditions during summer promote strong positive local air temperature gradients that drive persistent, negatively buoyant downslope winds, known as katabatic or glacier winds (Ball, 1956; Manins and Sawford, 1979). Katabatic winds are characterized by strong near-surface wind shear, large temperature gradients, and a shallow wind speed maximum (WSM) which can be below the standard measurement height (2 m) on even relatively gentle slopes (e.g. ~4$^\circ$ in Denby, 1999). Wind shear, represented in the bulk method through friction velocity ($u_\ast$), will diminish to zero as the height of a wind maximum is approached (Denby and Greuell, 2000). The closer to the surface a WSM is, the shallower the constant (variations less than 10 $\%$) momentum flux layer will be, limiting the region in which the theory will be valid (Denby and Smeets, 2000). Although the stability-based corrections in the bulk methods can correct for the effect of strong stability, they cannot account for the presence of a WSM, leading to a relatively poor performance of these methods during shallow katabatic flows (Fitzpatrick and others, 2017; Radić and others, 2017). These findings led to questioning of the validity of standard parameterizations in the bulk methods and to a development of potential alternative parameterizations (Radić and others, 2017).

To evaluate the bulk methods in simulating heat fluxes it is necessary to have reference measurements representing the true fluxes. In the absence of direct measurements of turbulent heat exchanges with the surface, the sensible (latent) heat fluxes are calculated as a covariance of high frequency measurements of wind speed and temperature (water vapour) through the eddy covariance method. Eddy covariance (EC) is a common technique in micrometeorology that requires relatively complex data processing as well as sensor maintenance to ensure that the key assumptions underpinning these techniques are being met (Foken and others, 2012). However, the installation and power requirements of the EC sensors, along with difficulties in fulfilling the necessary measurement assumptions, have limited the use of EC systems at glacier surfaces and the length of usable datasets where measurements exist.

The EC community has developed a set of best practices to improve the robustness of EC-derived fluxes and ensure consistency between various studies (e.g. Lee and others, 2004; Aubinet and others, 2012). The validity of the EC method is based on the assumption that the flow is fully turbulent and that measured fluctuations are solely attributed to eddy motion (Foken and others, 2012). The temporal averaging window from which the covariance is calculated must be short enough to avoid contamination by non-turbulent motions, but also long enough to capture motions of large eddies. In many applications, the spectral gap separating turbulent motions from changes in mean flow is assumed to be approximately 30 min (Stull, 1988). Thus 30min is the most commonly chosen interval length for the covariance calculations. However, in the presence of strong stratification and low wind speeds, as is often the case at glacier surfaces, the spectral gap can be on the order of minutes (Vickers and Mahrt, 2003; Mott and others, 2020; Nicholson and Stiperski, 2020). Furthermore, the optimal interval length is shown to be highly dependent on flow characteristics (Sun and others, 2018), while covariances assessed from a highly variable flow show strong sensitivity to the choice of the interval length, even at a scale of minutes or seconds (Mahrt and others, 2015).

The established best practices for EC data processing, such as the use of constant 30 min interval lengths over the observational period, have been generally adopted as a standard at glacier surfaces (e.g. Conway and Cullen, 2013; Fitzpatrick and others, 2017). However, flow conditions and turbulence characteristics at glacier surfaces can vary substantially over a melting season, with mean near-surface wind speeds spanning an order of magnitude. Furthermore, constant interval lengths cannot account for conditions that change during the 30 min interval, such as a sudden burst of cross-glacier wind or a shallow WSM moving past the sensor height. Additionally, because the constant flux layer is suppressed during shallow katabatic winds, EC-derived turbulent heat fluxes are unlikely to be representative of surface conditions when measurements are collected close to the WSM, which can be at or below the standard measurement height of 2 m. These observations highlight the need for further improvements of processing methods for EC data at glacier sites.

A key motivation for this study is to address the understudied role of EC processing methods in deriving the turbulent heat fluxes at glacier surfaces. Our first objective is to improve the processing methods to ensure the validity of calculated fluxes for conditions such as highly variable flow and low-level wind speed maxima. Our second objective is to re-examine the performance of the most commonly used bulk methods relative to the EC-derived fluxes under different near-surface flow regimes. To address these objectives, we utilize a two month EC and meteorological dataset measured at three different heights (1 m, 2 m, and 3 m) at a glacier site in the Yukon, Canada. The improved EC data processing methods are aimed to ensure that: (1) the covariances are derived from the interval lengths optimized as a function of flow characteristics displayed throughout the observational period, and (2) the EC-derived fluxes are representative of surface conditions, i.e. those well below the wind speed maxima. In developing these methods, we prioritize their applicability to one-level EC measurements, thus making the methods independent of multi-level EC measurements or any type of atmospheric profiling of wind and temperature.

2. Study site

The Kaskawulsh Glacier is a large, ~50 km long, temperate mountain glacier in Kluane National Park that drains from the St. Elias Icefields in the Yukon, Canada. The St. Elias Icefields are the largest non-polar icefields, and melting in this region is responsible for roughly 9 $\%$ of observed sea level rise in the latter half of the twentieth century (Arendt and others, 2002). The St. Elias Mountains are characterized by substantial topographic variations, with Canada's second tallest mountain located less than 20 km from sea level. This pronounced orography generates strong environmental gradients of temperature and precipitation from the Gulf of Alaska to the Yukon interior. Kaskawulsh Glacier accounts for roughly 9 $\%$ of glacier ice volume in the Yukon (Farinotti and others, 2019). The glacier has an estimated geodetic mass balance of -0.46 m water equivalent per year between 2007–2018 (Young and others, 2021), which is in line with the average rate of glacier thinning for this region (Berthier and others, 2010). Between 1956-2007, the glacier's terminus retreated by 655 m (Foy and others, 2011). Its meltwater contributed to Łhú’áán Män (Kluane Lake) through Ä’äy Chú (Slims River) until 2016 when its retreating terminus rerouted runoff into the Gulf of Alaska (Shugar and others, 2017). This rerouting has caused the Yukon's largest lake to drop by multiple metres and has increased the dust output from the now-empty Slims River, agitated by the persistent down-glacier katabatic winds (Bachelder and others, 2020). The size of Kaskawulsh Glacier facilitates the study of katabatic winds, as large glaciers produce relatively frequent and strong winds due to their larger fetch and consequent resistance to disturbances from other wind systems (Ohata, 1989).

3. Data

3.1 Measurements

An automated weather station was installed near the confluence of the north and central arms of Kaskawulsh Glacier (60$^\circ$45.517’N, 139$^\circ$07.513’W) at an altitude of 1666 m above sea level (Fig. 1). The automated weather station was comprised of two quadpods separated by approximately 5 m perpendicular to the glacier flow-line. Continuous measurements using a set of meteorological sensors at multiple heights above the glacier surface were made from June 30 to August 27, 2019. The local slope angle is 1.4°, and our microtopography surveys indicate that the surface roughness is approximately homogeneous in all directions within 500 m of the site. The measurements allowed for the estimate of surface melting through a surface energy balance (SEB) model accounting for all relevant fluxes: net shortwave and longwave radiation, turbulent heat fluxes, flux into the snow/ice and flux due to rain.

Figure 1. (Left) Map of confluence of the north and central arms on Kaskawulsh Glacier with the regional map in the bottom right corner. The automated weather station (AWS) is indicated by a white circle and the primary direction of glacier flow is indicated by two arrows. (Right) Setup of the field installation including Main I and Main II and the data-logger structure with solar panels.

One quadpod (Main I) recorded meteorological variables, and the other (Main II) recorded 20 Hz EC measurements (Fig. 1). Main I was equipped with aspirated temperature and humidity sensors installed at 1 m, 2 m, and 3 m, wind sensors at 2 m and 3 m, a four-component radiometer, a rain gauge, and a thermistor array drilled to an initial depth of 4 m, with thermistor spacing of 25 cm in the upper 1 m, and 1 m in the lower 3 m. Thermistor measurements were made every 30 s and averaged over 30 min. All other meteorological measurements were made every second and averaged every minute. Main II had three sonic anemometors sensors installed at 1 m, 2 m, and 3 m, all operating at a frequency of 20 Hz. An IRGASON anemometer, which also has an open path sensor for detecting fluctuations in water vapour, was installed at a height of 1 m, oriented directly up-glacier to reduce flow interference by the quadpod and minimize flow distortion effects observed by Horst and others (2016). Two Gill R3-50 anemometers were installed at heights of 2 m and 3 m. Prior to their installation in the field, the three sensors were tested for any biases in their measurements when simultaneously operating at the same site and at the same height. At Main II, the IRGASON was aligned parallel with the primary axis of the glacier, and the Gill anemometers were each placed at a 60$^\circ$ offset to minimize the interference between the sensors and crossarms. A final structure was equipped with two SR50 ultrasonic distance sensors to measure surface lowering. Campbell Scientific CR3000 dataloggers for both stations were installed approximately 5 m down-glacier from the two quadpods on a separate infrastructure to minimize any interference in measurements. A camera was installed in the vicinity of the station to record the conditions at the site every three hours throughout the observational period.

All the sensors were installed over relatively homogeneous terrain. Bare ice was exposed under the sensors for the entire campaign. The station operated autonomously over the observational period, with no need for manual readjustments or maintenance. Because of the robustness of the quadpods, the surface melting caused minimal effects on the alignment or tilt of the sensors as detected by the inclinometer (total change of 3$^\circ$ over two months). At the infrastructure with the dataloggers (Fig. 1), the shading from the logger boxes caused inhomogeneous melting resulting in a localized hummocky surface.

4. Methods

Our methods can be summarized as follows: first, we perform a clustering analysis to establish the most prevalent near-surface flow regimes, based on the multi-level measurements of wind and temperature at the study site. Second, we investigate the impact of EC data processing on the calculated turbulent heat fluxes. In particular, we investigate the use of different methods in determining the optimal interval length for calculating covariances, and we propose a filtering method for detecting EC data representative of surface conditions, i.e. those well below the WSM. Thirdly, we model the turbulent heat fluxes using the most commonly utilized aerodynamic bulk methods at glacier surfaces. We aim to quantify the effect that processing and filtering of EC data has on the EC-derived fluxes, as well as on the evaluation of modelled fluxes. This analysis is performed over the whole observational period and for each of the identified near-surface flow regimes. Throughout this paper, fluxes are defined in accordance with glaciological convention: positive (negative) fluxes denote a flux directed toward (outward) the surface.

4.1 Identification of near-surface flow regimes

We aim to identify near-surface flow regimes as characterized by mean profiles of wind speed and temperature in the first 3 m above the surface. To do so, we cluster the standard (30 min) measurements of wind speed and temperature from the three measurement heights. Prior to the clustering, the dataset is ‘compressed’ to variables that carry the bulk of the variance over the whole observational period. This ‘compression’ is achieved through principal component analysis (PCA), a standard method for dimensionality reduction and identification of dominant modes of variability within a dataset (Hsieh, 2009). PCA is applied to the whole dataset consisting of 30 min averages of: downslope wind, cross-slope wind, and temperature at 1 m, 2 m, and 3 m, as well as the gradients (differences) of downslope wind, cross-slope wind, and temperature between 3 m and 2 m, and 2 m and 1 m (e.g. u ₃ − u ₂ and u ₂ − u ₁). This yields 15 total variables. Prior to PCA, each of the 15 variables is standardized to give zero mean and unitary standard deviation. Once the dominant modes of variability are identified, each represented by an eigenvector and principal components (PCs), we focus only on the first few modes that collectively carry the bulk (>90%) of the variance in the data. The PCs of these selected modes are then clustered using agglomerative hierarchical clustering with Ward's method (Ward, 1963). The method recursively clusters data points by grouping the points with the highest similarity (smallest Euclidean distances) into bigger clusters while limiting the increase in inter-cluster variance at each step. This sequential procedure of merging smaller clusters into larger ones is represented by an ‘inverted tree’, or dendrogram. The initial large number of clusters (bottom of the inverted tree) yields smaller, more specific clusters, while the merged bigger clusters (top of the inverted tree) are more generic, ultimately leading to one cluster that contains all data points. While there is no objective way to determine the optimal number of clusters for the given dataset, a visual inspection of the dendrogram allows for an informed guess of the optimal number of clusters (Hsieh, 2009).

Once the clusters are identified, we assign a name to each cluster or regime. Each name is associated with a potential driver (e.g. katabatic) or a key characteristic of each flow (e.g. downslope). To support the analysis of potential drivers, in addition to our meteorological observations we also look into synoptic sea level pressure maps reconstructed from ERA5 reanalysis (Hersbach and others, 2020) for this region. Note that in the absence of observed wind speed and temperature profiles above 3 m, we are limited in our analysis to provide a more in-depth flow regime characterization.

4.2 EC data processing: interval length

EC data are used to compute sensible heat flux (Q _H) through covariance of high frequency vertical wind speed (w) and potential temperature (θ),

(1)$$Q_H = -\rho_ac_p\overline{w'\theta'},\; $$

where the prime denotes the turbulent component (as a deviation from the temporal mean) and the overbar denotes a temporal mean. ρ _a is air density and c _p = 1005 Jkg⁻¹K⁻¹ is the specific heat capacity of air. Note that the negative sign is added so that the sign of the flux agrees with the glaciological convention (positive flux into the surface) while still defining positive w′ as directed away from the surface. The raw (20 Hz) data underwent a series of preprocessing steps using the EddyPro data package (Fratini and Mauder, 2014) as outlined in Fitzpatrick and others (2017). The primary components of this preprocessing are spike removal (Vickers and Mahrt, 1997), planar fit coordinate rotation (Wilczak and others, 2001), and the Schotanus correction to account for effects of humidity (Schotanus et al., 1983). We applied additional preprocessing to account for high frequency and low frequency flux loss, following the methods of Ibrom and others (2007) and Moncrieff and others (2004), respectively.

First, we focus on the choice of interval length and temporal averaging window over which the covariance is calculated. A generalized procedure to accommodate for a spectral gap separating turbulent motions from changes in mean flow is to split 30-minute periods with n samples (observations) into m subsamples, split at locations in the timeseries τ _{1 : m} = (τ ₁,  …,  τ _m) and calculate the block-average covariance as (Howell and Mahrt, 1997),

(2)$$\overline{\zeta'\mu'} = {1\over n}\sum_{i = 1}^{m + 1}( \tau_{i}-\tau_{i-1}) \overline{\zeta'\mu'}_{( \tau_{i-1}\, \colon \, \tau_i) },\; $$

where ζ and μ can be any of u, v and w (three components of the wind speed vector), or θ. The subscript (τ _i−1 : τ _i) denotes a covariance calculated between τ _i−1 and τ _i. The most common approach is to calculate fluxes over a fixed window length, i.e. duration Δτ over which the measurements are taken (Δτ = τ _i+1 − τ _i), often set to 30 minutes. Here, we test three approaches for setting the subinterval length or the averaging window, and we refer to these methods as: (1) 30 min intervals, (2) Multiresolution-Flux Decomposition (MRD), and (3) Changepoint Detection (CPD):

30 min intervals: Covariances are calculated using $\Delta \tau = {30} {\min }$. Most calculations of turbulent fluxes from EC data on glaciers employ this method (e.g., Cassano and others, 2001; Conway and Cullen, 2013; Litt and others, 2017; Radić and others, 2017).

MRD: This method was originally introduced by Howell and Mahrt (1997) as a data analysis tool to assess time scales that are dominant contributors to the flux. Rather than setting a fixed 30 min interval, the MRD method determines the optimal average length scale Δτ from the time-scale-dependent contributions to covariance measurements. Following Vickers and Mahrt (2003), a record of 2^M EC data (e.g. w and θ) points is partitioned into averages containing 1,  2,  4,  ...,  2^M consecutive data points. We truncate our 30 min record to 27.3 min, containing 2¹⁵ = 32768 20 Hz data points. First, the lowest-order average, containing all 2^M data points, is subtracted from the record. The next-lowest-order mode comprised of two averages of 2^M−1 data points is then removed. The process is repeated for each mode and can be interpreted as a series of successive high pass filters. At each stage, the record is split into 2^M/2^M−m segments, where m = 0,  1,  2,  3...,  M. The filtered data are averaged over each of these segments, leaving a record with 2^m data points. As an example, if applied to 20 Hz w and θ data, taking the covariance of the filtered records with 2⁸ data points will yield the contribution of the 2⁸ (1/20) s=12.8 s timescale to the calculated $\overline {w'\theta '}$. The iteration over m = 0,  1,  ...,  M generates estimates of covariance as a function of averaging timescale. In the ‘covariance versus timescale’ plot, the zero-crossing of the covariance curve indicates the optimal gap scale Δτ for calculating the covariances for the entire observational period. The assumption is that covariances calculated over timescales smaller than Δτ are the result of turbulent motions while covariances calculated over timescales larger than Δτ are the result of non-turbulent motions and are thus omitted. More succinctly, MRD can be viewed as successive applications of Haar wavelets (Howell and Mahrt, 1997).

CPD: This method was originally used as a optimization technique in order to identify where statistical properties of a time series change (Scott and Knott, 1974). We adopt it here in order to account for potentially varying optimal interval length throughout the observational period. We apply CPD on the high-frequency time-series of u, v, w and T to automatically isolate turbulent motions from non-constant flow structures, such as turbulent rolls, breaking non-linear mountain waves aloft, cross slope winds, or the shallow WSM crossing the measurement height. In CPD, a set of candidate changepoints are tested by evaluating four-dimensional (u,  v,  w,  T) distributions before and after introducing the changepoint (for schematic illustration see Fig. 2). If the distributions on either side of the candidate changepoint are sufficiently similar (according to a kernel-based cost function), the changepoint is rejected, as it likely does not represent a change in a physical process. If the distributions on either side of the candidate changepoint are sufficiently dissimilar (based on a threshold for the kernel-based cost function), the changepoint is accepted. Following the notation of Killick and others (2012), we assume our data to be an array of the form y _{1 : n} = (y ₁,  …,  y _n) where y _k = (u _k,  v _k,  w _k,  t _k)^T, containing m changepoints at locations τ _{1 : m} = (τ ₁,  …,  τ _m) that are used to compute the covariance in Eqn. (2). The i ^th changepoint corresponds to the slice of data $y_{( \tau _{i-1} + 1) \, \colon \, \tau _i}$. The principle of the method is to select the number of changepoints and their locations in order to minimize

(3)$$\beta f( m) + \sum_{i = 1}^{m + 1}{\cal C}( y_{( \tau_{i-1} + 1) \, \colon \, \tau_i}) ,\; $$

where $\beta f( m) = \beta \log m$ is a penalty to prevent over-fitting. ${\cal C}$ is a cost function whose value is informed by the expected data distributions. For example, if the distributions of data are expected to be normal with a changing mean, then ${\cal C}$ = L ²-norm is sufficient to detect the changepoints. However, as our data do not adhere to an a priori distribution, we employ non-parametric kernel-based detection (Garreau and Arlot, 2016; Truong and others, 2020). The original y is mapped by features ϕ onto a reproducing kernel Hilbert space ${\cal H}$ implicitly through the Gaussian radial basis kernel k,

(4)$$k( p,\; q) = \exp( -\xi\Vert p-q\Vert _2^2) ,\; $$

where ξ > 0 is the bandwidth parameter, set by the inverse median of all pairwise distances between parameters in 12-dimensional space (our four variables measured at three heights). We select the Gaussian kernel due to its smoothness and popularity in machine-learning applications when little is known about the data distribution a priori (Truong and others, 2020). Regardless of the kernel chosen, the cost function in a kernel-based detection on an interval I = τ _i : τ _i+1 is:

(5)$${\cal C}_I = \sum_{\,j\in I} k\left(\phi( y_j) -\overline{\phi( y_I) },\; \phi( y_j) -\overline{\phi( y_I) }\right).$$

The changepoint locations are found using the pruned exact linear time algorithm (Killick and others, 2012). Practically, a direct implementation of the pruned exact linear time algorithm is computationally expensive because the kernel k contains over one billion elements for each 30 min period. We find that coarse-graining the columns and rows of k by a factor of ten (summing over 10x10 blocks within the kernel to reduce its size by two orders of magnitude) speeds up the pruned exact linear time algorithm by 5–10 times compared to current implementations (Truong and others, 2020) with no change in the method's performance. We perform CPD on the 12-dimensional input data comprised of u, v, w, and T at three heights for ease of comparison between heights, but this technique is also applicable to EC measurements made at only a single height.

Figure 2. Simplified schematic example of changepoint detection applied on a set of two variables (u, T). Here we assume one changepoint has already been established (top panel) and test two candidate changepoints. For visual clarity, only the two-dimensional (u, T) distributions are shown, while in the study we use the four-dimensional data (u, v, w, and T) at each of the three heights.

4.3 Intercomparison of the EC-processing methods

To analyse differences in EC-derived Q _H due to the three processing methods, we investigate characteristic sensible heat flux profiles, i.e. profiles along the three measurement heights, across the whole observational period and for each of the identified near-surface flow regimes. In particular, we want to examine how the occurrence of different characteristic profiles varies across the processing methods. The characteristic profiles are identified using self-organizing maps (SOMs), an unsupervised machine learning method that clusters the data on a two-dimensional map (Kohonen, 1982). A general feature of such a map is that more similar patterns are placed closer together on the map while more dissimilar patterns are placed further apart. The observed 30 min flux profiles, used as input data to the SOM algorithm, are normalized by dividing each profile by the measured Q _H at 1 m such that the normalized Q _H (z = 1 m) is 1. Observations where the measured Q _H at 1 m is less that 5 Wm⁻² are discarded to not skew the results toward the profiles that have a division with a small number. For consistency, if an observation is discarded for one processing technique, then the corresponding observation is also discarded for the other two processing techniques.

4.4 EC data filtering: detection of a WSM from one-level measurements

As mentioned earlier, the presence of a WSM close to the measurement height, as is often the case during shallow katabatic flow, leads to a misrepresentation of surface fluxes by EC-derived fluxes (van der Avoird and Duynkerke, 1999; Denby and Smeets, 2000; Finnigan, 2008). Here we introduce a method that detects the presence of a WSM from one-level EC measurements, so that these data segments can be omitted or filtered out from the calculation of covariances. The idea behind this filtering builds upon the methods of Grachev and others (2016) and is based on expected scatter plots of temperature versus wind speed at different heights relative to a WSM (below, at, and above a WSM) that is present during stable atmospheric conditions under which the katabatic flow develops and persists (see Fig. 3 for schematic illustration). We consider the following three theoretical cases for a stably stratified flow over a melting glacier surface (i.e., at 0$^\circ$C) in summer where air temperature increases with height. In all the cases presented, the displacement of an air parcel is considered over a small vertical distance.

1. 1. Below the WSM (Fig. 3A): Here, wind speed, like temperature, increases with height. A parcel of air displaced upward and away from the glacier surface (w′ < 0) will be colder (T′ < 0) and slower (u′ < 0) than its new surroundings, so $\overline {u'T'}> 0$. Similarly, a parcel of air displaced downward and toward the glacier (w′ > 0) will be warmer (T′ > 0) and faster (u′ > 0) than its new surroundings, so again $\overline {u'T'}> 0$. The positive covariance between u′ and T′ implies that the outline of the u′ − T′ scatter cloud can be approximated by an ellipse that has a positive angle ($\vartheta$) between its semi-major axis and the (T′) axis. As the vertical gradient of wind speed and temperature increases, so does the ratio (η) between the ellipse's semi-major and semi-minor axis (Fig. 3A relative to 3E). Thus, η is expected to approach 1 as the vertical gradients vanish (3D).

2. 2. Above the WSM (Fig. 3C): Here, wind speed decreases with height while temperature increases with height. A parcel of air moving upward and away from the glacier surface will be colder (T′ < 0) and faster (u′ > 0) than its new surroundings, so $\overline {u'T'}< 0$. A parcel of air moving downward and toward the glacier will be warmer (T′ > 0) and slower (u′ < 0) than its surroundings, so again $\overline {u'T'}< 0$. The outline of the scatter cloud of u′ − T′ measurements can be approximated by an ellipse with $\vartheta < 0^\circ$ and η > 1. Here, we implicitly assume temperature stratification above the jet is not strong enough to suppress turbulence.

3. 3. Near the WSM (Fig. 3B): Here, a parcel of air displaced upward and away from the glacier surface will be colder than its new surroundings (T′ < 0), but its horizontal speed will experience a negligible change (u′ ≈ 0) because ∂u/∂z = 0 at the WSM. Similarly, a parcel of air displaced downward and toward the glacier surface will display T′ > 0 and u′ ≈ 0 relative to its new surroundings. In both cases, $\overline {u'T'}\approx 0$, implying an ellipse with $\vartheta \approx 0$. However, the presence of a temperature gradient will produce an ellipse with η > 1.

Figure 3. Schematic example of scatter plots for assumed downslope wind speed (u) versus assumed temperature (T) for the case where the measurements are taken at a height below, at and above the wind speed maximum (WSM; left panel), as well as for the case where the WSM, if present, is well above the measurement height (right panel). Schematic profiles of wind speed (solid line) and temperature (dashed line) are also shown for the two cases. (A), (B), and (C) show the assumed measurements below, at, and above the WSM, respectively. (E) shows the assumed measurements close to the surface where gradients of T and u are relatively large, and (D) shows the assumed measurements far from the surface where gradients are low.

We hypothesize that the measurements representative of surface conditions are those taken below the WSM (case 1 above) and therefore should be detected by the following characteristics: a positive covariance ($\overline {u'T'}> 0$), positive angle $\vartheta$ between the ellipse's semi-major axis and the horizontal axis ($\vartheta > 0$), and a semi-major to semi-minor axis ratio greater than unity (η > 1). We introduce thresholds on positive values of $\vartheta$ and η, or in other words:

(6)$$ \eqalign{ & \matrix{0^\circ< \vartheta_{{\rm low}}< \vartheta< \vartheta_{{\rm high}}< 90^\circ,\; }\cr & \,1< \eta_{{\rm low}}< \eta,\; }$$

where $\vartheta _{{\rm low}} = 25^\circ$, $\vartheta _{{\rm high}} = 65^\circ$, and η _low = 1.3 are selected after testing a range of values on our data. The selection of these threshold values is based on striking a balance between data quality and data retention. Selecting more strict criteria (e.g., narrowing the range of acceptable $\vartheta$ and η) further improves data quality but reduces data quantity.

4.5 Evaluation of bulk methods

Our objective here is to evaluate the most commonly-used aerodynamic bulk methods in their estimates of turbulent heat fluxes using the EC-derived fluxes as our reference data. At its core, a bulk aerodynamic method is rooted in gradient transport theory or K theory, in which the turbulent fluxes of momentum and sensible heat (Q _H) are proportional to the time-averaged vertical gradients of wind speed (u) and temperature (T), respectively (Stull, 1988). The multi-level meteorological measurements, as collected in this study, would allow for the application of a profile ‘bulk’ method that relies on differences between two-level measurements. This profile ‘bulk’ method, however, is known for large errors (Denby and Smeets, 2000; Hock, 2005), a result that is also corroborated by our data (not shown). Thus we focus only on the bulk methods based on one-level meteorological measurements. Although there are many variants of the bulk method, mainly related to the stability corrections used, the three most often employed on glacier surfaces are those tested by Fitzpatrick and others (2017): the bulk method without any stability corrections, the bulk method with stability corrections using the bulk Richardson number (hereafter the ‘bulk Richardson method’), and the bulk method with stability corrections using the Obukhov length. Initially, we evaluated all three methods against the EC-derived fluxes, but for the brevity of the paper we choose to focus on the method that overall performed the best, which is the bulk Richardson method. The bulk Richardson correction has the additional advantage of relying only on mean meteorological variables and not Obukhov length L, which has been criticized for use on glaciers beyond serving as a proxy for local stability through z/L (e.g. Grisogono and others, 2007; Monti and others, 2014). We note that all conclusions based on the results with this method also hold for the other two methods.

The bulk method for deriving Q _H at glacier surfaces is based on the mixing-length theory by Prandtl (1935), which assumes that friction velocity ($u_\ast$) and wind speed (u) at a given measurement height z are linearly related by a dimensionless exchange coefficient C _v,

(7)$$u_\ast = C_v( z) u( z) ,\; $$

while the expression for sensible heat flux Q _H is derived as:

(8)$$Q_H = \rho_a c_P{1\over {\rm Pr}}C_t( z) u_\ast ( T( z) -T_0) .$$

Here, Pr is the Prandtl number (Pr = 0.7, from Pope, 2000), C _t is the dimensionless exchange coefficient for temperature, T ₀ is the temperature at the glacier surface (often set to 0$^\circ$C), and T(z) is the temperature at measurement height z. $u_\ast$ is the modelled friction velocity from Eqn. (7). The dimensionless exchange coefficients for momentum (C _i = C _v) and temperature (C _i = C _t) are modelled as

(9)$$C_i = \kappa/\ln( z/z_{0, i}) ,\; $$

where κ = 0.4 is the von Kármán constant, z _0,v is the roughness length for momentum, and z _0,T is the roughness length for temperature. The roughness lengths are derived from our EC data following Radić and others (2017). We derive a separate roughness length for each measurement height and for each EC processing and filtering technique, as the roughness lengths are fitting parameters that represent the dynamic effects of the surface on momentum and heat transfer, and thus are not necessarily directly related to the true surface roughness (Sun and others, 2020). z _0,v is calculated from Eqns. (7) and (9) as

(10)$$z_{0, v} = z \exp( -\kappa u/u_\ast ) ,\; $$

where $u_\ast$ is calculated from the EC data as

(11)$$u_\ast = \sqrt{\overline{u'w'}^2 + \overline{v'w'}^2}.$$

Similarly, z _0,T is calculated from Eqn. (8), incorporating the expression for Q _H as derived from the EC data (Eqn. (1)) and using the EC-derived $u_\ast$:

(12)$$z_{0, T} = z \exp( \kappa u_\ast ( T( z) -T_0) /\overline{w'\theta'}) .$$

The above expressions for z _0,v and z _0,T theoretically hold for neutral stability, and thus the EC data are filtered to ensure that $u_\ast$ and Q _H are only considered during near-neutral conditions. The stability is assessed through the EC-derived Obukhov length L, and we omit strongly stable and unstable conditions by restricting measurements to |z/L| < 0.1. In addition to this filter, a series of other filtering steps is applied to ensure high quality data. For completeness, we list the filters briefly here, but refer the reader to Radić and others (2017) for a more detailed explanation of the filters. The filtering steps employed are:

1. 1. Wind direction filter: Restrict incident wind direction to ±45$^\circ$ of the central axis of the EC sensor.

2. 2. Temperature filter: Restrict measurements to $T( z) > {1}^\circ$C as errors in deriving roughness lengths are comparatively large for small temperature gradients.

3. 3. Realistic value filter: restrict z _0,v and z _0,T to between 10⁻⁷ m and 1 m.

These filtering steps are applied to all EC-derived fluxes, regardless of the processing method used (30 min, MRD, or CPD) prior to calculating the roughness lengths.

To account for the suppression of turbulence and reduction of flux due to the prevalent strong near-surface stratification, we employ the bulk Richardson correction when calculating Eqn. (8). Following Webb and others (1980), the exchange coefficients of Eqn. (9) become

(13)$$C_i = \kappa/\ln( z/z_{0, i}) ( 1-5{\rm Ri}_b) ,\; $$

where Ri_b is the bulk Richardson number, calculated as:

(14)$${\rm Ri}_b = {gz( T( z) -T_0) \over T( z) u( z) ^2},\; $$

with temperature expressed in Kelvin, and gravitational acceleration g = 9.81ms⁻². To evaluate the performance of the modelled Q _H relative to EC-derived Q _H, we compute a root-mean-square error (RMSE), correlation (r), and mean bias error (MBE) for each of the three heights. The EC-derived Q _H and $u_\ast$, as well as the roughness lengths, are calculated using the three processing methods (30 min averages, 1 min averages, and CPD) as well as the ellipse filtering method for the presence of the WSM. All fluxes are averaged to 30 min for ease of comparison with previous studies, regardless of the processing method used.

To assess uncertainties due to measurement errors in the calculations of roughness lengths and sensible heat flux we follow the standard methods for error propagation of multivariate functions (Bevington and others, 1993). To quantify error in the roughness lengths, we assume constant measurement errors in u and T of δu = 0.3m s⁻¹ and $\delta T = {0.1}^\circ$C, respectively (Table 1). From Laubach and Kelliher (2004), we assume relative errors in measured $u_\ast$ and $\overline {w'\theta '}$ of 5 $\%$. Once we quantify the relative errors in roughness lengths for each 30 min interval, we determine the mean relative error in z _0,v and z _0,T for the whole observational period. These errors, together with the errors in u and T, are then propagated in the calculation for Q _H (Eqn. (8)).

Table 1. Instrumentation used in this study and their manufacturer-stated accuracy

Sensors installed on Main I collected low frequency measurements (1 Hz and slower) and those installed on Main II collected high frequency measurements (20 Hz).

5. Results

5.1 Clusters of flow regimes

For a summary of the meteorological conditions observed over the study period, we refer to fig. S1 and accompanying text. Here we first present the results of our clustering algorithm, as the evaluation of measured and modelled fluxes will be performed across these clusters, as well as for the entire observational period. The first four of the 15 modes of PCA are found to explain 97 $\%$ of the variance in the mean data, enabling a significant reduction of dimensionality (Fig. S2). The first mode (explaining 58 $\%$ of variance) is mainly represented by the variability in the downslope wind speed (u) and temperature (T) at all three measurement heights. The first mode displays positively correlated downslope wind speed and wind shear, temperature, and temperature gradient. The second mode (19 $\%$ of variance) shows downslope wind speed and wind shear anticorrelated with temperature and temperature gradient. The third mode (11 $\%$ of variance) shows correlated cross-slope wind and wind shear, and the fourth mode (8 $\%$ of variance) shows anticorrelated temperature and temperature gradient. Although we perform PCA on u, v, and T at three heights and two finite differences per variable, we find similar results – in terms of wind speed and temperature carrying the bulk of the variance – when only inputting into PCA one measurement height and one finite difference per variable. The variance percentages differ slightly (up to 5% for the first mode), but the key features of each eigenvector are the same. Similar results are obtained when the same analysis is performed on mean values in each subinterval established through MRD and CPD. We initially included slope-normal velocity w as an input variable, but after analysing the results, found that variations in w between regimes were small and not significant until higher-order modes, which explained little variance. Therefore, we omit w in the interest of simplicity.

Performing hierarchical clustering on the data in the four-dimensional principal component space produces six clusters, where the number of optimal clusters is determined from the dendrogram (fig. S3). For each regime (cluster) we plot the cluster-averaged wind profile and temperature profile, as well as the distribution of each cluster's occurrence within a day (Fig. 4). Downslope flow primarily originates from the northern arm of the glacier (Fig. 1). Below we list the six regimes with their assigned names and briefly describe their key characteristics. The clusters are listed in descending order according to their frequency of occurrence over the observational period.

‘Downslope’ regime: the most frequent regime (34 $\%$ of data points associated with this cluster) is characterized by persistent downslope winds (with mean 2 m wind speed of 4.0ms⁻¹) with moderate near-surface temperature gradients (mean gradient of 2.5$^\circ$C between 1 m and the surface). This regime occurs primarily (57% of the time) between midnight and noon (Fig. 4).

‘Strong downslope’ regime: the second most frequent regime (22 $\%$ of data) is most prominent between noon and midnight during clear-sky conditions. The regime displays strong downslope winds (mean 2 m wind speed of 5.7ms⁻¹) and strong near-surface temperature gradients (mean gradient of 5.4$^\circ$C between 1 m and the surface) with small temperature gradients above 1 m (Fig. 4). According to our ellipse filtering method applied to the data from this and the ‘downslope’ regime, the presence of a WSM is not detected within the first 3 m above the surface. However, this does not mean that the two regimes are not associated with deeper katabatic flow.

‘Katabatic’ regime: the third-most prevalent regime shows occasional ellipse-flattening at 3 m ($\vartheta < 25^\circ$ observed in 41 $\%$ of this regime), implying the presence of a nearby WSM not far above 3 m. This regime exhibits moderate temperature gradients (mean gradient of 3.4$^\circ$C between 1 m and the surface) with low wind speeds (mean 2 m wind speed of 2.0ms⁻¹; Fig. 4). Due to the ellipse flattening indicating the presence of a WSM above 3 m, relatively weak winds, and strong temperature gradient, we call this regime katabatic.

‘Cold synoptic’ regime: this regime (11 $\%$ of data) occurs during episodes associated with storm conditions (rainfall and snowfall) that took place in the middle of August. This regime is characterized by low wind speeds (mean 2 m wind speed of 1.8ms⁻¹) and a constant-with-height air temperature of approximately 0$^\circ$C (Fig. 4). We label this regime as ‘cold synoptic’ since it coincides with the passage of cold fronts according to the synoptic pressure maps for this region (not shown).

‘Shallow katabatic’ regime: this regime (9 $\%$ of data) is characterized by a WSM below 3 m and strong temperature gradients across the WSM (mean gradient of 3.4$^\circ$C between 1 m and 2 m; Fig. 4). Two thirds of all shallow katabatics are observed between noon and midnight. An example of how CPD and ellipse flattening is used to observe the presence of a WSM in this regime is shown in fig. S4.

‘Upslope’ regime: upslope flow accounts for the remainder of the data (5 $\%$ of data) and occurs most often late at night (Fig. 4). This regime is only observed on seven days, with wind direction exclusively up-glacier, strong near-surface wind gradients and moderate temperature gradients. Note that because of the alignment of the IRGASON sensor to measure downslope wind, the EC measurements at 1 m are likely not valid for this regime.

Figure 4. Mean vertical profiles of wind speed (top panel) and temperature (middle panel) for each of the six flow regimes (clusters). Frequency of occurrence (f) of each regime over the observational record is above each column. Shaded regions show the standard deviation derived from the measurements associated with each cluster. Number of times (counts) each regime is observed as a function of time of day (bottom panel). Black lines are the raw counts and coloured lines show the smoothed curves (running averages). Time of day is given in local time (Mountain Standard Time, UTC -7 h).

We also test the use of different numbers of clusters according to the same dendrogram (fig. S3). Adding a seventh and eighth cluster splits the cold synoptic regime into three different regimes that only vary slightly in incident wind direction and temperature profile. These clusters occur infrequently and do not meet the ellipse filtering conditions or the general data quality filters of Radić and others (2017). Collapsing to five regimes instead of six aggregates upslope flow and cold synoptic regimes, despite one having neutral and the other stable stratification. Further collapsing to four clusters combines downslope and katabatic regimes, even though the latter shows the presence of a WSM near 3 m, so important information is lost by selecting fewer than six regimes. Thus, six is the optimal number of clusters required to capture the dominant flow regimes in our measurements.

5.2 Fluxes from processed EC data

We process the EC data using the three methods with different interval lengths for covariance calculation (30 min, MRD and CPD) in order to derive sensible heat fluxes. MRD finds an optimal interval length of 1 min in our measurements (fig. S5). In CPD, 30 min records are split, on average, into 10 subintervals, with two records that are split into 20 subintervals and two records that are not subdivided at all. The shortest subinterval is 12.5 s long, while the average subinterval is 3 min long (see fig. S6 for the distribution of subintervals). The temporal variability in the optimal interval length is most pronounced in the ‘katabatic’ and ‘shallow katabatic’ flow regimes, with averages of 14.0 and 12.7 subintervals per 30 min record, respectively. The variability is least pronounced in the ‘downslope’ and ‘strong downslope’ flow regimes, with averages of 6.9 and 7.4 subintervals per 30 min record, respectively. In the case of the ‘cold synoptic’ regime, the standard 30 min method yields similar 30 min fluxes as MRD (1 min interval length) and CPD. Applying MRD to each flow regime provides similar results (fig. S8). The ‘katabatic’ and ‘shallow katabatic’ flow regimes exhibit height-dependent gap scales and suggest a shorter mean averaging window, while the gap scales from the ‘downslope’ and ‘strong downslope’ regimes are longer and show little height dependence.

Next, we examine how the characteristic profiles of measured sensible heat flux depend on the flux processing method using SOMs. After testing different numbers of clusters (size of the SOM), we settle on a 2 × 3 SOM, i.e. a map showing six characteristic flux profiles determined from all three processing methods (Fig. 5). We expect theoretically physical profiles of sensible heat flux in the surface layer to show either a very small dependence on height, or a monotonic decrease with increasing height. Clusters (nodes) #2 and #6 of the SOM fit these criteria of a physical flux profile, with node #2 showing a small dependence on height and node #6 showing a monotonic decrease (Fig. 5). The remaining profiles are theoretically unphysical, with nodes #1 and #3 showing non-monotonic flux profiles, and nodes #4 and #5 showing fluxes that increase significantly as a function of height. Thus, we hypothesize that these profiles are likely observed due to inadequate EC data processing when determining EC-derived Q _H. As each observation is associated with one cluster (characteristic flux profile), we calculate the frequency of occurrence of each cluster across the dataset from each processing method separately. For 30 min averages, only 33 $\%$ of all observations are associated with the theoretically physical profiles (22 $\%$ with node #2 and 11 $\%$ with node #6). For the MRD method, 46 $\%$ of all observations are associated with the theoretically physical profiles (34 $\%$ with node #2 and 12 $\%$ with node #6). Finally, when processing fluxes with CPD, 76 $\%$ of all observations are associated with the theoretically physical profiles (53 $\%$ with node #2 and 23 $\%$ with node #6).

Figure 5. Six clusters (#1 to #6), presented as a 2 × 3 self organizing map (SOM), of sensible heat flux profiles computed from the data with all three eddy covariance (EC) processing methods: 30 min, multiresolution flux decomposition (MRD), and changepoint detection (CPD). The SOM is calculated using profiles from all three processing techniques, but the frequency of occurrence of each cluster is calculated separately for each processing technique. The three percentages above each cluster present the frequency of occurrence of that cluster for 30 min, MRD, and CPD processing, when read from left-to-right. The profiles with shaded grey backgrounds are those deemed theoretically unphysical as the flux increases with increasing measurement height, either from 1 m to 2 m, or from 2 m to 3 m.

To further look into the differences in profiles across the three processing methods, we analyse differences in the daily running mean of Q _H between a pair of heights, i.e. Q _H(3 m) − Q _H(2 m), and Q _H(2 m) − Q _H(1 m), over the whole observational period (Fig. 6). Calculating fluxes with the 30 min method, the positive gradient of Q _H between 1 m and 2 m has a maximum of 31.2Wm⁻². In comparison, using MRD gives a maximum positive gradient of 15.2Wm⁻², and processing with CPD gives a maximum positive gradient of 4.6Wm⁻². The percentage increase of Q _H between 1 m and 2 m exceeds 10 $\%$ for 36.0 $\%$ of the record when processing with 30 min averages, 28.6 $\%$ of the record when processing with 1 min averages, and 5.7 $\%$ of the record when CPD is used.

Figure 6. Differences in EC-derived sensible heat fluxes between 3 m and 2 m (black) and between 2 m and 1 m (red). EC data are processed with 30 min method (top), MRD 1 min interval length (middle), and CPD (bottom). Fluxes are smoothed with a 1-day moving average. Grey shading indicates periods where the flux at 2 m exceeds the flux at 1 m by more than 10 $\%$, provided the absolute value of the flux at 1 m exceeds 5Wm⁻².

We also compare the percentage change of sensible heat flux with height when fluxes are averaged across each of the six identified flow regimes (Table 2). Looking across all the regimes, the 30 min method yields heat fluxes that, on average, increase between 1 m and 2 m, but decrease between 2 m and 3 m. A similar pattern with a maximum flux at 2 m is observed when EC data are processed with the MRD method (1 min interval length). When data are processed using CPD, fluxes are found to decrease monotonically with height in the downslope, strong downslope, katabatic, and shallow katabatic regimes. According to the CPD method, the differences in Q _H between heights is shown to be statistically significant (to a significance level of 0.05) in each of these four regimes except between 1 m and 2 m in the katabatic regime. The cold synoptic regime does not show a statistically significant difference in Q _H between heights. However, we note that the mean sensible heat flux in the cold synoptic regime is approximately 3Wm⁻², so the small absolute differences in Q _H present as large relative differences. The upslope regime shows a statistically significant flux increase between 1 m and 2 m for all processing methods, but the flow in the upslope regime at 1 m is obstructed by the quadpod, making the EC-derived Q _H at 1 m likely erroneous.

Table 2. Median percentage change of sensible heat flux between heights, as calculated with the three flux processing methods (30 min, MRD, CPD) for six identified flow regimes

Negative (positive) percentages denote a decrease (increase) with increasing height. Shaded cells indicate that the difference between the compared fluxes is statistically significant at a significance level of 0.05, as assessed by the two-sample t-test for equality of the means.

5.3 Fluxes from filtered EC data

Here, we present the results of our ellipse filtering method that ensures the EC-derived fluxes are representative of surface fluxes: The ellipse filtering is computed with both MRD and CPD at each measurement height. Applying ellipse filtering criteria on EC data processed with $\vartheta _{{\rm low}} = 25^\circ$, $\vartheta _{{\rm high}} = 65^\circ$, η _low = 1.3 retains 52 $\%$, 42 $\%$, and 34 $\%$ of the high frequency data at 1 m, 2 m, and 3 m, respectively, for both CPD and MRD. Here, the ellipse filtering omits the following data scatters: scatters with negative ellipse angle (corresponding to Fig. 3C), flat ellipse angle (Fig. 3B), or ambiguous ellipse angle because $\eta \ \approx 1$ (Fig. 3.D). This filter additionally omits data with vertical ellipse orientation (wind speed varies while temperature does not), and more ambiguous wind-temperature scatters that do not resemble an ellipse (e.g., η = 1, $\vartheta = 0$). Examples of the characteristic u′ − T′ scatters from our EC measurements are presented in fig. S7 in the supplementary material, where T′ and u′ are computed as deviations from the 30 min means. Although the percentage of total data that pass the ellipse filtering criteria is similar between MRD and CPD, the two methods do not agree on which 30 min records pass the filtering criteria. For example, the overlapping 30 min records for which at least 25 min of the record pass the ellipse filtering criteria for both CPD and MRD occurs 82 $\%$ of the time at 1 m, 88% at 2 m and 92% at 3 m. We also analyse how well the identified periods with WSM height below 3 m, according to the ellipse filtering method, agree with those from standard wind speed measurements at three heights. We find that for 97% of those time segments, according to the ellipse filtering method, the measured wind profiles also indicate a WSM below 3 m.

As the goal of ellipse filtering is to ensure measurements reflect surface conditions, we expect measurements that pass ellipse filtering to have fluxes which are roughly constant in height (variations less than 10 $\%$), which is one of the conditions defining a surface boundary layer (Stull, 1988). Here, we compare the EC-derived Q _H, and its variability with height, as derived with and without the ellipse filtering. Using MRD (1 min interval length), the relative mean bias error between Q _H at 2 m and 3 m without the data filtering is 22.5 $\%$, while after the filtering it is decreased to 10.2 $\%$ (Fig. 7). For the CPD method, the same relative MBE is decreased from 28.4 $\%$ to 3.4 $\%$, implying that the ellipse filtering is more effective when applied with CPD than with the MRD method. When applied to measured Q _H between 1 m and 2 m, both MRD and CPD produce relative mean bias error of less than 10 $\%$, although MRD yields fluxes that reach maximum values among the three heights at 2 m, while for CPD the maximum fluxes are observed at 1 m, as is theoretically expected. These findings imply that CPD with ellipse filtering is the most successful among the methods in identifying the data whose fluxes vary by less than 10 $\%$ in the first 3 m above the surface. In the following analysis, we restrict the use of ellipse filtering to CPD intervals only.

Figure 7. Comparison of EC-derived sensible heat fluxes between 2 m and 3 m (top) and 1 m and 2 m (bottom) using 1 min MRD-derived interval length (left) and variable CPD interval lengths (right). Grey dots show all 30 min records and pink dots show 30 min records that pass the ellipse filtering criteria. Statistical metrics (RMSE in W m⁻², mean relative bias error (MRBE), and correlation coefficient r) are shown for both cases.

Over the observational period, the difference between the mean energy available for surface melting, as assessed by a surface energy balance model at our site (see Supplementary Material for details), increases by 0.2% when the EC-derived Q _H is used with CPD method at 1 m relative to the EC-derived Q _H with standard 30 min method at 2 m. At daily scales, however, the difference in EC data processing can lead to a difference in estimated melt energy by up to 5%, while at hourly scales the difference can be up to 20%. We also compared the modelled melt to the observed melt as inferred from the surface lowering measured by the SR50 sonic rangers. Relative to the standard (30 min) EC-derived Q _H at 2 m, we find that the improved (CPD with ellipse filtering) estimate of Q _H at 2 m can reduce the bias between modelled and observed melt by 10–25% at sub-daily scales.

5.4 Modelled versus EC-derived fluxes

In this section, we show the results of the bulk method evaluation, first performed over the whole dataset and then across the six flow regimes. We start, however, by analysing the relationship between measured wind speed (u) and EC-derived friction velocity ($u_\ast$) because the bulk method for assessing the momentum flux is grounded in this relationship. According to Eqn. (7), there should be a linear relationship between the two variables, with the slope equal to the dimensionless exchange coefficient C _v. We assess when the $u-u_\ast$ scatter resembles this linear relationship depending on the processing method used (30 min, MRD, CPD, and CPD with ellipse filtering) at each measurement height (Fig. 8). The results show that the linear fit to $u-u_\ast$ scatter is performs worse as the measurement height increases from 1 m to 3 m (Fig. 8). At all heights, however, the linear fit substantially improves if ellipse filtering is applied to the CPD-processed data (R ² improving by at least 20 $\%$ relative to 30 min averages).

Figure 8. Scatter plots of 30 min averaged u versus EC-derived $u_\ast$, each at 1 m (bottom), 2 m (middle), and 3 m (top) for the four EC processing techniques: 30 min, MRD (1 min), CPD, and ellipse-filtered CPD. Grey points indicate all data and pink points denote the data that pass the filtering criteria of Radić and others (2017, percentage given by n _p). In the fourth column, n _p is the percentage of data that pass the filtering of Radić and others (2017) and the ellipse filtering criteria. The dashed black line shows the average $u_\ast$ for each bin interval of u, with a bin width of Δu = 0.5ms⁻¹. The red line shows the trendline derived from a linear regression on the pink points, while a coefficient of determination (R ²) for the fit is indicated in the bottom-right corner of each plot.

Next, we compute the roughness lengths z _0,v (Eqn. (10)) and z _0,T (Eqn. (12)) from the data that pass the filters listed in the Methods section. At 1 m, the mean logarithmic z _0,v ranges between 10^−2.8 m and 10^−3.1 m depending on processing method used (Fig. 9). Above 1 m, estimates of z _0,v vary more, ranging between 10^−2.5 m and 10^−3.5 m, depending on the EC processing technique and height. The mean z _0,T varies between 10^−4.8 m and 10^−5.2 m among the three heights. The scatter (standard deviation from the logarithmic mean) in momentum and temperature roughness length increases with height for all EC processing techniques. When testing the performance of the bulk method in simulating sensible heat fluxes, for each height and each EC processing method we use the mean estimates of $\log {z_{0, v}}$ and $\log {z_{0, T}}$ as derived in Figure 9.

Figure 9. Probability density function (PDF) of EC-derived momentum (z _0,v; red) and temperature (z _0,T, blue) roughness lengths for four EC processing techniques at 1 m (bottom), 2 m (middle), and 3 m (top). The vertical dashed line denotes the mean in log-space and temporal variability is given by ± one standard deviation.

We calculate the mean relative error, i.e. the ratio in the error of roughness length to the roughness length, δz ₀/z ₀, where δz ₀ is derived through the propagation of errors as explained in the Methods section. For momentum, the mean δz _0,v/z _0,v varies between 0.60 and 0.75 depending on measurement height and processing technique selected, with no height nor processing technique providing a systematic advantage over any other. A standard deviation in the relative error, as assessed over the whole observational period, varies between 0.1 and 0.25. A mean momentum roughness length of z _0,v = 0.001 m with a mean relative error of 0.75 can be expressed as 10^−3.0±0.7 m. The mean δz _0,T/z _0,T varies between 0.45 and 0.6, with a standard deviation in the relative error ranging between 0.05 and 0.15. We note that the magnitude of the errors in z _0,v and z _0,T is of the same order of magnitude as the temporal variability (one standard deviation) in our EC-derived roughness lengths (Fig. 9). We use the mean relative error of 0.69 for z _0,v and 0.51 for z _0,T in the assessment of errors in modelled Q _H by the bulk method.

The evaluation of the bulk Richardson method in simulating Q _H over the whole observational period (Eqn. (13)) shows the worst performance when the standard 30 min covariances are used at each height (Fig. 10). Q _H is overestimated at each height, with the largest overestimation (MBE =13.0Wm⁻²) at 3 m, followed by 2 m (MBE =7.8Wm⁻²), and then by 1 m (MBE =4.3Wm⁻²). As EC-processing complexity increases (from 30 min to MRD to CPD to ellipse-filtered CPD), the overestimation in Q _H decreases at all heights, to a minimum of 7.1Wm⁻², 3.0Wm⁻², and 0.6Wm⁻² at 3 m, 2 m, and 1 m respectively when applying ellipse-filtered CPD. A similar trend is observed in RMSE: the error decreases closer to the surface (e.g., from 27.1Wm⁻² at 3 m to 14.1Wm⁻² at 1 m using 30 min averages) and as EC-processing complexity increases (e.g., from 14.1Wm⁻² using 30 min averages to 4.8Wm⁻² using ellipse-filtered CPD at 1 m). Similarly, the correlation coefficient (r) also increases closer to the surface and with increasing EC-processing complexity, with the smallest values of r = 0.76 at 3 m for the 30 min method, and the highest values of r = 0.98 at 1 m for CPD and ellipse-filtered CPD.

Figure 10. Modelled versus observed (EC-derived) 30 min sensible heat flux at 1 m (bottom), 2 m (middle), and 3 m (top). Solid lines are the bin-averaged Q _H, calculated by averaging the modelled fluxes that fall within each 5 Wm⁻² bin of observed fluxes. Dashed lines are the 1:1 lines. Light grey vertical lines show propagated measurement error. Root-mean-square error (RMSE, W m⁻²), mean bias error (MBE, in W m⁻²), and correlation (r) are shown for each case.

In each of the regimes, the correlation between EC-derived and modelled Q _H increases as the measurement height drops from 3 m to 1 m (Fig. 11). In the ‘downslope’ and ‘strong downslope’ regimes, both RMSE and MBE between EC-derived and modelled Q _H decrease as the height drops from 3 m to 1 m, while in the ‘katabatic’ and ‘shallow katabatic’ regimes, only RMSE decreases as the height drops, although we note that there are very few records that pass the ellipse filtering criteria in the shallow katabatic regime at 2 m and 3 m. Although the absolute MBE increases slightly closer to the surface in the katabatic regime (0.4Wm⁻² at 2 m to –1.3Wm⁻² at 1 m), the correlation improves significantly (from 0.34 at 2 m to 0.78 at 1 m).

Figure 11. Modelled versus observed (EC-derived) 30 min sensible heat fluxes in each of the six flow regimes at all three heights. Ellipse-filtered CPD is used for the EC-derived fluxes. Dashed line shows the 1:1 line. Light grey vertical lines plotted for each point show estimated uncertainty in the modelled Q _H. Statistical metrics (RMSE in W m⁻², MBE in W m⁻², and r) are shown for each case.

6. Discussion

6.1 EC data processing and filtering

In the comparison of three methods for covariance calculations (30 min, MRD, and CPD), we find that CPD best identifies appropriate flux window lengths for the whole observational period, as well as for the near-surface flow regimes. As our results demonstrate, the CPD processing yields profiles of heat fluxes that best agree with the theoretically expected profiles in the surface boundary layer. According to the clustering analysis with the SOMs (Fig. 5), CPD yields the largest number of observations that resemble the theoretically expected characteristic profiles of Q _H. This result is also corroborated by analysing the differences (and their statistical significance) in Q _H across the three measurement heights for each of the flow regimes. The same analysis also reveals that the conventional 30 min method for covariance calculation consistently performs the worst, i.e. yielding the smallest number of, and the least similarity with, the theoretically expected characteristic profiles of Q _H. While MRD performs better than the 30 min method, corroborating previous findings at glacier surfaces (e.g., Nicholson and Stiperski, 2020; Mott and others, 2020), CPD is still the preferred method for several reasons, listed below.

Firstly, in the presence of non-logarithmic wind speed profiles, (unfiltered) EC-derived fluxes are expected to be highest near the surface and decrease as a function of height closer to the WSM. This theoretically expected profile is most consistently observed with CPD, i.e. highest fluxes derived at 1 m and lowest fluxes measured at 3 m. On the other hand, MRD consistently calculates a flux that is maximized at 2 m. This finding is consistent across all flow regimes, except ‘shallow katabatic’. Although there have been reported differences in EC-derived fluxes between IRGASON and Gill sonic anemometers (e.g., Wang and others, 2016), we did not find any systematic bias during our sensor calibration and testing. Applying the correction from the previously reported differences between the two sensors (Wang and others, 2016) to the MRD data does not remove the unphysical flux maximum at 2 m. If these unphysical profiles are indeed an artefact due to measurements from different sensor manufacturers, then CPD may be a promising approach to circumvent this type of bias, but more targeted research is needed on this subject.

Secondly, CPD gives varying optimal lengths throughout the observational period, while MRD provides only one optimal length. As MRD calculates an optimal interval length in the frequency domain and CPD calculates an optimal interval length in the time domain, CPD is better suited to operations that act in the time domain, such as computing time-varying fluxes. To enable MRD to detect more than one optimal length, we test the application of the method separately to each of the six flow regimes and derive regime-specific optimal averaging lengths. By doing so, however, we encounter some practical challenges in determining the optimal length from the MRD curves due to the absence of their ‘zero-crossings’ (fig. S8). In addition, some of the MRD-determined optimal interval lengths are counter-intuitive. For example, MRD determines the optimal interval length at 1 m to be 1 min for the cold synoptic regime, and 15 min for the strong downslope regime. However, in terms of reducing both MBE and RMSE in Q _H between heights, the strong downslope regime benefits more than the cold synoptic regime if covariances are calculated with 1 min intervals (instead of 30 min).

Another question which has not yet been directly addressed in previous studies is the representation of true surface heat fluxes by EC-derived fluxes at glacier surfaces. While some studies have deployed EC sensors at roughly 1 m above the glacier surface (Munro, 1989) assuming that measurements closer to the surface better represent true surface fluxes, many used the standard height of 2 m above the surface. Our results demonstrate that 1 m EC-derived fluxes are indeed more representative of surface fluxes than 2 m EC-derived fluxes due to the presence of shallow WSM at this site. We also conclude that 3 m EC-derived fluxes are the least representative of surface conditions. However, we advise caution in installing an EC sensor at 1 m above the surface, especially in conditions where: the surface is particularly rough and the sensor path length is large (Burba, 2013), the sensor cannot be positioned to avoid flow distortions (e.g. Geissbühler and others, 2000; Horst and others, 2016), or the post-processing steps to avoid frequency loss substantially reduce the amount of valid data (e.g. Moncrieff and others, 2004; Ibrom and others, 2007). Nevertheless, as we demonstrate that with our ellipse filtering method, it is possible to process the EC data, collected at or above the WSM, to obtain EC-derived fluxes that are representative of surface fluxes. The ellipse filtering method works best in combination with the CPD averaging method, as determined by the similarity in EC-derived fluxes among the three heights.

We find that the visual guide of ellipses of u′ − T′ scatter, as proposed by our filtering method, more clearly identifies WSM presence when compared to only looking at $\overline {u'T'}$ as in Grachev and others (2016), especially when using EC data from only one measurement height. Since Grachev and others (2016) suggested that both multilevel measurements of $\overline {u'w'}$ and $\overline {u'T'}$ could be used to determine the height of the WSM, we additionally test the use of the ellipse filtering applied to u′ − w′ scatter. We find that the results from u′ − w′, in terms of the detection of a WSM presence, do not always agree with the results from u′ − T′, and more importantly do not always agree with wind profile measurements from standard wind sensors at our site. In comparison, these instances of disagreement between the WSM (below 3 m) detected through ellipse filtering on u′ − T′ and through the measured wind profiles are rare ($< 3\%$ occurence). Recent high-resolution simulations of sub- and super-critical jets by Salinas and others (2021) show zones of negative shear production above and below the WSM in which $\overline {u'w'}$ and ∂u/∂z are anticorrelated. These simulations may serve as an explanation for the observed poorer performance of filtering on u′ − w′ scatter relative to u′ − T′ scatter, although further analysis is required to confirm if negative shear production is being detected in our data.

While the ellipse filtering method ensures that EC-derived fluxes are representative of the surface fluxes, the method can substantially reduce the amount of high-quality data obtained. The reduction in data is especially striking for the shallow katabatic regime: only 22% of EC-derived fluxes at 1 m are representative of surface fluxes. For the same conditions, 2 m (3 m) fluxes are representative of the surface only 4% (2%) of the time. On the other hand, in the strong downslope flow regime, 76%, 68%, and 57% of EC-derived fluxes are representative of surface fluxes, at 1 m, 2 m, and 3 m, respectively. By its design, ellipse filtering will retain more observations with higher fluxes during strong wind regimes and discard low fluxes during weak winds with near-surface WSM. Overall, the use of CPD with ellipse filtering, relative to the standard 30 min method with no filtering, can lead to a difference in estimated melt energy by up to 20% on sub-daily scales. This difference is similar to the reported error in the SEB closure that used the standard 30 min covariance calculation in deriving Q _H at several mid-latitude glaciers (Fitzpatrick and others, 2019). Thus, at sub-daily timescales, the use of CPD with ellipse filtering can improve the assessment of sensible heat fluxes and simulation of surface melt.

6.2 Bulk method performance

The correct representation of surface fluxes by adequately processed EC data enables some new insights in the performance of the bulk methods. A few previous studies found that the standard bulk method overestimates $u_\ast$ during shallow katabatics (e.g., Radić and others, 2017). This result is not surprising, as EC-derived $u_\ast$ approaches zero at the WSM. We show here that the EC-derived $u_\ast$ close to the WSM is not representative of the surface $u_\ast$, which features in the bulk method estimate of Q _H. The bulk method is designed to represent $u_\ast$ near the surface and not $u_\ast$ near the WSM. Thus, the systematic overestimation of sensible heat fluxes observed in previous studies is likely attributed to their reference data not being representative of the surface conditions, and not a failure in the bulk method as was previously proposed (e.g., Conway and Cullen, 2013; Radić and others, 2017). Although we derived these results using the bulk Richardson method, our key conclusions do not change when other commonly-used bulk methods are applied: Regardless of the flow regime, the bulk method performance in simulating Q _H is better if the measurements of temperature and wind speed are taken at 1 m above the glacier surface instead of at the standard 2 m height. More importantly, even for the measurement heights above 1 m and in the vicinity of a WSM, the bulk method is shown to perform well as long as the reference EC-derived fluxes are adequately processed, i.e. with the use of CPD and ellipse filtering.

While improved processing of EC-derived heat fluxes leads to a better match with the modelled fluxes for each of the flow regimes, some biases between modelled and observed fluxes still remain. In particular, the bulk method overestimates Q _H during the ‘downslope’ and ‘strong downslope’ regimes, with the largest overestimation for measurements taken at 3 m. Likely, the overestimation stems from the misalignment of the observed scatter trend between u and $u_\ast$ and the linear relationship assumed by Eqn. (7). The observed linear trendline in the $u-u_\ast$ scatter does not have zero intercept as expected by Eqn. (7) and the mixing-length theory. Instead, the $u-u_\ast$ scatter would be better suited to a piece-wise linear, or ‘hockey stick’, fit where $u_\ast$ is roughly constant as u increases from zero until some velocity threshold is reached and then linearly increases with u beyond this velocity threshold. The ‘hockey-stick’ fit in $u-u_\ast$ scatter has been observed previously over non-glaciated surfaces (often above 10 m) and is attributed to the suppression of turbulence generation due to strong stratification (e.g., Sun and others, 2012; Sunand others, 2015; Freundorfer and others, 2019; Grisogono and others, 2020; Sun and others, 2020). Since the linear fit in the $u-u_\ast$ scatter determines the value of C _v in Eqn. (7) which then features in the calculation of Q _H in Eqn. (8), it is possible to empirically adjust this fit to a piecewise-linear fit that better represents the ‘hockey stick’ pattern and thus improve the calculations of C _v and Q _H. We attempted this bias correction, which led to a decrease in MBE for Q _H to <3Wm⁻² from the original 13Wm⁻² observed at 3 m using 30 min mithout without ellipse filtering. However, as the bias correction is empirical and likely site-specific, for ease of comparison with previous studies and with little understanding of what physically drives this bias, we refrain from proposing these corrections to glacier sites in general.

The roughness lengths at ice-exposed glacier surfaces have been shown to vary substantially from one location to the next (e.g. van den Broeke and others, 2005; Brock and others, 2010). As an important control on turbulent flux generation, the roughness lengths require accurate representation in the bulk method. Our EC-derived z _0,v = 10^−2.8±0.7 m and z _0,T = 10^−5.0±0.7 m, using the 30 min method for processing EC data at 2 m, are in line with previous findings over glaciers with bare ice exposed (e.g., Munro, 1989; Brock and others, 2006; Fitzpatrick and others, 2019). The ratio of z _0,v/z _0,T ≈ 100 is also consistent with previous findings (Hock, 2005). Nevertheless, the estimates of roughness lengths are shown to vary by up to an order of magnitude depending on the measurement height and the EC processing method used. For example, using CPD with ellipse filtering applied to data at 2 m yields z _0,v = 10^−3.4±0.6 m and z _0,T = 10^−5.2±0.6 m. We note that the temporal variability (one standard deviation) in these estimates across the whole observational period is of the same order or larger than the uncertainty of individual 30 min estimates of z _0,v and z _0,T. The relatively large temporal variability, especially for z _0,v, is likely due to the fact that z _0,v reflects the total dynamic effect of the surface on momentum transfer, and thus may not represent solely the physical surface roughness that is relatively constant over the observational period (Sun and others, 2020). As EC-derived roughness lengths are the most commonly used reference values when evaluating other techniques for deriving z _0,v, such as those developed from photogrammetry and remote sensing (Fitzpatrick and others, 2019), the relatively large sensitivity in z _0,v to the choice of EC data processing method and measurement height should be taken into consideration.

7. Conclusions

The primary objectives of this study are to: (1) improve the EC data processing methods, targeted for one-level measurements, to ensure the validity of calculated fluxes for conditions such as highly variable flow and low-level wind speed maxima, and (2) evaluate the most commonly used bulk methods relative to the EC-derived fluxes under different near-surface flow regimes. To that end, standard meteorological and EC measurements were collected at three different heights (1 m, 2 m, and 3 m) at a site on the Kaskawulsh Glacier in the Yukon over a two-month period in summer 2019. We summarize our key findings as follows:

The length of the time window over which covariances between wind speed and temperature are computed from EC data has a substantial impact on the EC-derived fluxes at hourly and daily scales. By intercomparing the three methods – standard 30 min method, 1 min interval length as derived by Multiresolution-Flux Decomposition (MRD), and our proposed Changepoint Detection (CPD) method – we find that the CPD method best determines the optimal averaging window that varies throughout the observational period and produces physically realistic near-surface profiles of sensible heat flux. Although the difference between MRD and CPD is small for observations taken at 1 m, the differences are larger at 2 m and 3 m. As most previous studies on glaciers have installed sonic anemometers at or above 2 m, CPD may be able to improve the flux measurements of these previous studies.

We propose a filtering method applied to one-level EC data to ensure EC-derived fluxes are computed from measurements representative of surface conditions. The filtering method can successfully remove the data ‘contaminated’ by the presence of the WSM at or in the vicinity of the measurement height.

With the CPD and filtering methods applied to the EC data, the agreement between modelled and EC-derived sensible heat fluxes is substantially improved relative to standard processing methods, at each measurement height. This agreement also holds during the shallow katabatic flow regime, directly contradicting previous findings which highlighted the failure of the bulk method and asked for improved theory. We show that the standard theory works provided EC data are adequately processed, and provide a processing procedure to ascertain when the bulk method can be relied upon, even in the presence of highly variable wind speed and a maximum wind speed at or below the measurement height.

EC measurements taken at 1 m above the surface more frequently pass our filtering criteria than those at 2 m or 3 m, implying that measurements at 1 m are more representative of surface conditions. Relative to the measurement heights above, the bulk method at 1 m shows the least scatter, the best correlation, and the smallest bias from the reference EC-derived fluxes for the whole observational period, as well as for different flow regimes.

Our results highlight that future assessments of turbulent heat fluxes on glaciers should prioritize adequate EC data processing that goes beyond the standard practices initially established for non-glacierized flat terrain. Although the results presented in this study show an improvement in both deriving turbulent heat fluxes from EC data and establishing a better agreement between EC-derived fluxes and those modelled through bulk methods, it remains to be seen how transferable these findings are to other glaciers. As Kaskawulsh is a very large mountain glacier, we suspect the observed near-surface flow regimes to differ from those at smaller mountain glaciers, and we suggest a similar analysis be performed prior to a long-term installation of EC systems at glacier surfaces.