A Vertical Propeller Eddy-Covariance Method and Its Application to Long-term Monitoring of Surface Turbulent Fluxes on the Greenland Ice Sheet

On the Greenland ice sheet, the sensible heat flux is the second largest source of energy for surface melt. Yet in atmospheric models, the surface turbulent heat fluxes are always indirectly estimated using a bulk turbulence parametrization, which needs to be constrained by long-term and continuous observations. Unfortunately, such observations are challenging to obtain in remote polar environments, especially over ablating ice surfaces. We therefore test a classical eddy-covariance method, based on propeller anemometers and thermocouple measurements, to estimate the momentum and sensible heat fluxes on the Greenland ice sheet. To correct for the high-frequency attenuation, we experimentally derive the sensor frequency-response characteristics and evaluate the universal turbulence spectra on the ice sheet. We show that the corrected fluxes are accurate and that the sampling interval can be reduced to 4 s to increase the system’s autonomy. To illustrate its potential, we apply the correction to one year of vertical propeller eddy-covariance measurements in the western ablation area of the ice sheet, and quantify the seasonal variability of the sensible heat flux and of the aerodynamic roughness length.


Introduction
The total mass balance of the Greenland ice sheet, defined as the integrated surface mass balance minus the calving of ice at marine-terminating glaciers, is a primary component of the global sea-level budget. Between 2012 and 2016, the ice sheet lost on average 247 Gt yr −1 of mass (≈ 0.7 mm yr −1 sea-level equivalent), which accounts for 37% of all the land-ice contribution to global sea-level rise (Bamber et al. 2018). This recent strong mass imbalance of the ice sheet has been linked to a significant increase in surface melt , which is either measured in-situ or calculated by closing the surface energy balance, B Maurice van Tiggelen m.vantiggelen@uu.nl where M is the surface melt, R net is the net absorbed radiation by the surface, H is the sensible heat flux, L E is the latent heat flux, and G is the ground heat flux. Here we define H , L E, and G positive upwards and express them in W m −2 . On the ice sheet, a positive R net drives most of the surface melt, while L E and G are an order of magnitude smaller (Kuipers Munneke et al. 2018b). The sensible heat flux H , however, is an important source of energy for the melt of seasonal snow in mountain regions (Mott et al. 2011) and in the Arctic tundra (Pohl et al. 2006), but also for the melt of Arctic sea ice (Tjernström et al. 2015) and at the surface of Antarctic ice shelves (Kuipers Munneke et al. 2018a). Despite efforts to measure the various components of the surface energy balance on the ice sheet Van As et al. 2011;Kuipers Munneke et al. 2018b), direct measurements of turbulent heat fluxes are still limited. Instead, an indirect bulk method is typically used to estimate the turbulent surface fluxes from measurements made by singlelevel automatic weather stations. This has revealed that the sensible heat flux is also an important source of energy for surface melt on the ice sheet, both in the western ablation area (Kuipers Munneke et al. 2018b) and the southern ablation area (Fausto et al. 2016). However, the modelled sensible heat flux using these methods can be highly uncertain, either due to underlying assumptions (Radić et al. 2017) or due to physical parameters that are not well constrained in time and space, such as the aerodynamic roughness length (Smeets and Van den Broeke 2008).
One way to directly measure the turbulent surface fluxes uses the sonic eddy-covariance (SEC) method, based on fast measurements of the three-dimensional wind vector and temperature acquired with sonic anemometers. Such instruments are costly, require a continuous and significant power supply, and do not function under drifting snow conditions or precipitation. This makes them less than practical for long-term experiments in remote polar areas. Yet several experimental campaigns have successfully measured the turbulent surface fluxes with sonic anemometers on the ice sheet (Henneken et al. 1994;Forrer and Rotach 1997;Box and Steffen 2001;Smeets and Van den Broeke 2008;Miller et al. 2017;Madsen et al. 2019). Unfortunately these datasets rarely span more than several weeks, and are not always representative of areas with the highest surface melt rates.
A more feasible alternative to measuring turbulent fluxes on the ice sheet is the vertical propeller eddy-covariance (VPEC) method, which relies on propeller anemometers and thermocouples (Blanford and Gay 1992). One then faces two practical obstacles: the icing of the instruments, and their limited frequency response (Horst 1997). Conveniently, icing is not a frequently occurring problem in the katabatic wind zone, where the air is usually undersaturated (Smeets et al. 2018). This paper aims to provide a solution to the second challenge: the high-frequency attenuation of the measured fluxes due to the limited frequency response of the propeller anemometers and thermocouples.
Although the high-frequency attenuation of the measured fluxes has received quite some attention, the latest developments mainly focus on experimental set-ups that use a slowresponse sensor in combination with a fast sensor to measure fluxes of atmospheric gases, such as carbon dioxide and water vapour (Ibrom et al. 2007), or methane (Peltola et al. 2013). We focus on the high-frequency attenuation of momentum and sensible heat fluxes caused by applying a combination of two slow-response sensors. Our aim is to accurately model this attenuation, in particular under the stable conditions commonly observed in the atmospheric boundary layer over the ice sheet. . The set-up of Experiment 3 is not shown here but can be found in Lenaerts et al. (2014). The numbers indicate different instruments: (1) Fine-wire thermocouples (also indicated by arrows), (2) Gill vertical propeller anemometer, (3) Young horizontal propeller anemometer, (4) CSAT3 sonic anemometer, (5) Intelligent weather station and logger, (6) CNR4 net radiation sensor, (7) Ablation draw-wire sensor, (8) CR1000 datalogger The paper is organized as follows: in Sect. 2 we give an overview of the field experiments and the instruments, in Sect. 3 we give a detailed description of the model that is used to correct for the high-frequency attenuation of the fluxes. In Sect. 4 we then evaluate the model and the corrected VPEC fluxes against SEC fluxes, and we quantify the influence of a reduced sampling rate. Finally, we apply the high-frequency correction to one year of VPEC measurements made in the western ablation area of the ice sheet in Sect. 5, and quantify the temporal variability of the fluxes.

Instrumental Set-up
In the following experiments, we test the VPEC method. For the horizontal wind measurement we use a Young wind vane anemometer (model 05103-L, R.M. Young Company, Traverse City, Michigan, USA) fitted either with polypropylene or carbon fiber thermoplastic blades. For the measurement of the vertical wind speed we use a Gill vertical propeller fitted with expanded polystyrene blades (model 27106, R.M. Young Company, Traverse City, Michigan, USA). An alternative instrument for the VPEC method is the K-Gill propeller vane (Ataktürk and Katsaros 1989), which has a higher sensitivity to wind fluctuations. Unfortunately, this sensor requires sensitive material for both the propellers, which makes it unsuitable for long-term studies in polar conditions. Besides, it is more convenient to install a vertical Gill propeller next to the Young wind vanes that are already used on many existing weather stations.
The temperature is measured by a 76.2 µm diameter, type E, fine-wire thermocouple (model FW3, Campbell Scientific, Logan, USA). The propellers are positioned at the same height above the surface, and the centre of the blades are separated by 0.50 m. The fine-wire thermocouple referred to hereafter is fitted just below the blades of the vertical propeller.
As a reference, we apply the SEC method with a non-orthogonal sonic anemometer (CSAT3, Campbell Scientific, Logan, USA) sampling at a rate of 10 Hz. The sonic anemometer is positioned at 0.50 m distance from the Young wind vane anemometer and at the same height as the centre of the propeller blades.
When mentioned, net absorbed radiation by the surface is measured by a net radiometer (model CNR4, Kipp & Zonen, Delft, the Netherlands) and ice ablation is measured by a drawwire (model FD115 -15000, Althen, Leidschendam, the Netherlands). For the measurements presented below, the draw-wire sensor is positioned on the same mast as the other instruments, at a height of 3.08 m (see Fig. 1). A weight is attached to the tip of the draw-wire, which is drilled 10 m into the ice. The cumulative ablation of the ice surface causes the wire to roll around a spring-loaded spool, and the wire's linear extension is then measured every 30 min with a potentiometer (Hulth 2010).

Description of the Experiments
We use measurements from four separate field experiments, two of which were performed on the Greenland ice sheet. The characteristics of each experiment are summarized in Table  1 and are further detailed below.

Experiments 1 and 2: CESAR Site in 2011 and 2019
The full set-up consisting of both of the VPEC and SEC instruments was tested twice at the Cabauw Experimental Site for Atmospheric Research (CESAR, e.g., Monna and Bosveld 2013), located on a flat grassland in the Netherlands (51.970 • N, 4.927 • W,−0.8 m a.s.l). The first experiment took place during August and September 2011, while the second experiment took place during February 2019.
In the first experiment, all the time series were sampled at 10 Hz. In the second experiment, however, the sampling rate of the VPEC system was reduced to 2 Hz in order to test a different sampling of the horizontal wind speed. Differences between the two experiments also comprise the height and the orientation of the sensors, as well as the material of the horizontal propeller. All information can be found in Table 1.

Experiment 3: Site S10 in 2012
This experiment took place at site S10 of the K-transect (67.00 • N, 47.02 • W, 1850 m a.s.l), which is a transect of automatic weather stations and mass balance observations located in the ablation area of the western Greenland ice sheet. It spans from the ice edge up to 1850 m elevation and has been operated by the Institute for Marine and Atmospheric research Utrecht since 1993. Further details about the K-transect can be found in Smeets et al. (2018) and Kuipers Munneke et al. (2018b). The measurements for this specific experiment are documented by Lenaerts et al. (2014) and are used to test the validity of the Kaimal et al. (1972) turbulence spectra on the ice sheet. We use the raw SEC time series, which were Table 1 Instrumental set-up during the four different field experiments, performed at the Cabauw Experimental Site for Atmospheric Research (CESAR) located in the Netherlands (NL), and at two sites on the Greenland ice sheet (GrIS). z m denotes the measurement height above the local surface, d is the displacement height and h is the height of the local surface above the surrounding average topography; L u is the response length of the horizontal anemometer; A w is the calibration constant of the vertical anemometer; and τ T is the response time of the fine-wire thermocouple, which were used to correct for high-frequency attenuation and are defined in the body of the paper. Their measurement uncertainty is given in the brackets. Propeller types used are identified as follows: PP: polypropylene, CFT: carbon fibre thermoplastic, EPS: expanded polystyrene recorded four hours per day between August and October 2012. The comparison between SEC and VPEC fluxes cannot be done with this experiment, as there was no Young anemometer adjacent to the Gill vertical propeller. Furthermore the sensible heat fluxes are small, with an average of 20 W m −2 , and the frequency of both riming and blowing-snow events reduces the amount of valid data for comparison. The snow surface at this site is very homogeneous, and slopes downward from east to west with a slope angle of about 0.4 degrees. This gives rise to a south-easterly katabatic flow more than 70 percent of the time (Smeets et al. 2018). The height of the SEC instruments was recorded every 30 min by a sonic height ranger, and decreased from 4.9 m in August to 4.2 m in the end of October due to snowfall.

Experiment 4: Site S5 in 2016-2017
To illustrate the potential of the method, the VPEC instruments without a sonic anemometer have also been operated on the K-transect for a longer period of time. The measurements were made at site S5 of the K-transect (67.09 • N , 50.06 • W , 540 m a.s.l) between September 2016 and August 2017. All the instruments are connected to a low power logger which continuously records the time series at an interval of 4 s.
The local ice surface is composed of rough hummocks and domes, interlaced by meltwater streams. The station is located on the top of an ice dome, and the dome itself is located approximately 1.5 m above the average surrounding topography, denoted h. The effective measurement height is thus, where z m is the height of the sensors above the local surface, which was recorded every 30 min by a sonic height ranger and ranges between 3.7 m and 4.0 m during the measurement period. The displacement height d and the height of the local surface above the average surrounding topography h are assumed constant and equal to 0.5 m and 1.5 m, respectively, after Smeets and Van den Broeke (2008).

Notations
We work with the measured time series of the horizontal wind vector U = (u, v), vertical wind speed w, and air temperature T , where u and v are the along-wind and cross-wind components of U, respectively, and w is the wind component normal to the local surface slope. The 30-min average of x is written as x, while the 30-min covariance between the fluctuations of x and y is written as x y . The temperature flux is thus w T , and the along-wind momentum flux is u w . Hereafter we write the 30-min average of U as U .

Preliminary Corrections
The raw time series of each experiment are first screened for non-physical values, and the measured vertical wind speed from the vertical propeller anemometer is multiplied by a factor 1.25, in order to account for the non-cosine response of propeller anemometers at high angles of attack (Gill 1975). Then multiple iterations of a median absolute difference threshold filter (Mauder et al. 2013) are performed to remove individual spikes in the raw time series. The latter filter proved unnecessary for the propeller and thermocouple measurements, due to the very small amount of spikes in these time series. The raw time series are block-averaged in 30-min windows, linearly de-trended, and windows with more than 5% missing data are flagged. Furthermore, block-averaged time windows with wind directions that are suspected to contain flow distortion are also flagged. For both the VPEC and the SEC instruments, a yaw rotation followed by a pitch rotation, or double rotation (Kaimal and Finnigan 1994), is used to rotate the raw measured wind vector in the local horizontal reference frame, thereby correcting for the flux biases induced by tilted sensors. This rotation method was chosen over the planar fit method (Wilczak et al. 2001), as the tilt angles of the weather station in the ablation area of the ice sheet change over time due to melt. Besides, only a narrow band of downslope wind directions are available for analysis in the katabatic zone of the ice sheet. The final step of the preliminary processing involves the calculation of the raw turbulence (co)spectra, which are smoothed with an averaging window that exponentially expands with frequency (Kaimal and Finnigan 1994).
During Experiments 1 and 2, the buoyancy flux measured by the sonic anemometer is corrected for humidity influences according to Schotanus et al. (1983), using the latent heat flux measured at a nearby location by the Royal Netherlands Meteorological Institute (KNMI) at the CESAR observatory. 1 Finally the SEC fluxes are corrected for path averaging after Moore (1986), using a path length of 0.12 m.
The downward shortwave radiative flux measured during Experiment 4 is corrected for the pitch and roll of the net radiometer. For this we use the tilt angles measured by the inclinometer located in the station logger, and the geometrical equations found in, e.g., Wang et al. (2016). The tilt angle of the weather station varied between 4 • and 6 • towards the west over one year. This tilt orientation means that the correction mostly shifts the phase of the measured shortwave incoming radiation, which affects the daily averaged radiation by less than 2 W m −2 . The ablation draw wire measurement is not corrected for the movement of the station, as the horizontal displacement of the station with respect to the borehole was less than 0.5 m after one year. This offset results in an error in yearly ice ablation of less than 0.04 m, which is also the difference with the manual ablation measurement at a nearby stake and ≈ 1 % of the total yearly ablation. The ice ablation is converted to an energy flux using a latent heat of melting of 334 × 10 3 J kg −1 and a constant ice density of 916 kg m −3 .

Fluxes
When analyzing the fluxes measured during all four experiments, we minimize the influence of flow obstruction and propeller stalling at low wind speeds by applying the following data selection criteria: where WD is the wind direction of the 30 min-averaged wind vector and 'Flag' is the quality flag of the preliminary flux corrections. The interval of unobstructed wind directions is determined a priori by the design of the mast and the relative location of each instrument. It is then adjusted iteratively until an optimal trade-off between VPEC and SEC flux agreement and data quantity is found. The final chosen intervals for each experiment are given in Table  1. We do not include any filter related to the error in the cross-momentum flux v w , as we assume that the same error is present in both the SEC and the VPEC fluxes.

Variance Spectra
For the analysis of the variance spectra measured at the CESAR site during Experiments 1 and 2, we remove the ill-defined spectra in terms of signal-to-noise ratio, and thus extend the previous data selection with the following criteria:

Covariance Spectra
Finally, when analyzing the turbulence cospectra measured at site S10 on the ice sheet during Experiment 3, we limit the influence of drifting snow and apply a very strict near-neutral stability range: The resulting selected fluxes and spectra for each experiment are summarized in Table 1.

Methods: High-Frequency Attenuation Correction
When used to measure turbulent fluxes, propeller anemometers have the following limitations: (1) a limited frequency response, (2) a non-linear directional sensitivity, and (3) a threshold starting wind speed (Wyngaard 1981). In this section we provide a spectral correction for the limited frequency response. The data selection criteria are used to mitigate for limitations (2) and (3). The measured (co)variance between atmospheric quantities x and y, denoted x y , is a fraction A xy of the true (co)variance. We calculate the attenuation coefficient A xy after Moore (1986), where x y m is the measured (attenuated) (co)variance between atmospheric variables x and y, and f is the frequency (Hz). With this definition, the attenuation of the flux is equal to 1 − A xy . This method thus requires an expression for the reference turbulence (co)spectrum S xy and the total transfer function of the system T xy .

Model
The sensor transfer functions T xy are the product of both the low-pass and the high-pass filters, where T D denotes the high-pass filter caused by the block-averaging procedure, T p is the low-pass filter associated to the averaging of the x measurement along a path length p x . The latter filter is not used for the VPEC system. In the above, G x is the frequency response function of the x sensor, T s denotes the low-pass filter caused by the spatial separation of the x and y sensors by a distance s xy , and G xy is the frequency response of the covariance x y . For T s and T p we use the exponential expressions from Moore (1986). For T D we use the analytical expression from Moncrieff et al. (2005).
The frequency response of a propeller anemometer or a thermocouple is approximated as a first-order gain function, with a response time of τ x (Wyngaard 1981), With this definition of τ x , the cut-off frequency is f c = 1/(2πτ x ), such that G 2 x ( f c ) = 1/2. The a priori response time of a sensor is not known, but we will assume that the response time of a horizontal propeller anemometer is inversely proportional to the horizontal wind speed (MacCready and Jex 1964), where we define L u as a response length (or distance constant), which only depends on the physical sensor characteristics. We assume that the response time of the vertical propeller also depends on the angle of attack, which is defined as the angle of the instantaneous wind vector with the plane of rotation of the propeller (Fichtl and Kumar 1974), where A w is an empirical calibration constant depending on the propeller type, and σ w is the standard deviation of the vertical wind speed. Assuming that the phase difference between x and y is small and independent of frequency, the transfer function of the covariance is then written as (Horst 1997),

Experimental Derivation of Sensor Response Times
We derive the response times of the propeller anemometers and the thermocouples using the measurements acquired at the CESAR site during Experiments 1 and 2. This is done by fitting the square of a first-order gain function, G 2 x ( f ), to the experimental transfer functions using the non-linear least-squares Levenberg-Marquardt algorithm. The experimental transfer functions are defined as, with where S V P EC x x is the variance spectrum measured by the slower instruments and S SEC x x is the simultaneous spectrum measured by the sonic anemometer. We include a normalization coefficient N to force the lower-frequency part of the transfer function to be equal to one. The (1/(2πτ x )) = 0.5. The red and blue lines denote the whole range of fitted first-order functions. The '-4' ('+ 4' ) dashed lines denote the 6 dB per octave decrease (increase) with increasing frequency related to first-order roll-off and instrumental noise, respectively. Bottom panels: corresponding estimated response times of the vertical propeller anemometer (d), the horizontal propeller anemometer, (e) and the fine-wire thermocouple (f). The linear regression of the response time characteristics is denoted by the solid line, and the chosen uncertainty interval by the dashed lines limit frequency f is chosen such that the normalization coefficient is not affected by highfrequency attenuation (Aubinet et al. 1999). To reduce the influence of noise and aliasing, we .5] (Hz) frequency range. Finally, we reject half the estimated response times that yield the poorest fit residuals, with the aim of rejecting spectra with poor signal-to-noise ratios that result in non-physical response times. The resulting estimated response times are shown in Fig. 2. The associated error bar is taken differently for each instrument. For the vertical propeller the error bar is taken as the interval containing 80% of the estimated response times, while for the horizontal propeller it is taken as the interval containing 80% of the estimated values for wind speeds below 5 m s −1 . The thermocouple response time is assumed constant and taken as the average estimated response time after three weeks of operation in the field.
The heavier material of the propeller blades increases the response length of the horizontal anemometer from L u = 1.78±0.2 m for carbon fibre thermoplastic to L u = 3.15±0.2 m for polypropylene. These values differ significantly from the values reported by the manufacturer (2.2 m and 2.7 m, respectively). During Experiment 1, the horizontal wind speed was sampled at 10 Hz by counting the amount of revolutions of the horizontal propeller every 100 ms. This sampling method results in violet noise in the S uu spectrum, or noise increasing as f 2 , which is visible in Fig. 2b at frequencies above 1 Hz. During Experiment 2, the horizontal wind speed was sampled differently from Experiment 1 by averaging the time between all the propeller revolutions every 200 ms. This removes the noise and results in more well-defined first-order spectra. However, with the latter method, wind speeds less than one propeller revolution within a 200-ms interval cannot be measured. For the type of propellers used in this study this corresponds to wind speeds of less than 2 m s −1 .
For the vertical anemometer, only expanded polystyrene blades were tested, and the average response constant is equal to A w = 0.45 ± 0.1 m. In the surface layer, σ w /U typically ranges between 0.02 and 0.2, depending on stability and surface roughness. This means that the effective response distance of the vertical anemometer L w ≡ τ w U can range from 1.3 m to 6.1 m (see Eq. 7), depending on the average angle of attack. This is a known result (Fichtl and Kumar 1974;Garratt 1975), and partly explains why different wind-tunnel (McBean 1972;Hicks 1972;Gill 1975) and field experiments (Horst 1973) find a different value for L w . The drawback of applying Eq. 7 to estimate the τ w in time is that it requires a priori knowledge of σ w , which is underestimated by a propeller anemometer. As such, we simply assume that σ w = 2σ w,V P EC in Eq. 7, where σ w,V P EC is the uncorrected standard deviation of the vertical wind speed, measured by the vertical propeller anemometer. Conveniently, the spectra of the vertical anemometer do not contain high-frequency noise, as the vertical wind speed is sampled by measuring a voltage that is directly proportional to the propeller revolution speed.
The fine-wire thermocouple response could only be estimated during the first experiment because of the higher sampling rate. The response time does not show any significant variation with wind speed. It does however show an increase in time, from 0.08 s to 0.14 s after 6 weeks of operation in the field, which we attribute to accretion of material on the fine wires due to air pollution and rain. In the remainder of the paper we will assume τ T to be constant and equal to 0.13 ± 0.04 s.

Model
We assume that the normalized turbulence spectra follow the functions experimentally derived by Kaimal et al. (1972). Under stable stratification, i.e., z/L O > 0, these functions are written as, where n = f z/U is the normalized frequency. A more general relation is given by Horst et al. (2004) and by Massman and Clement (2005),  Kaimal et al. (1972) and the solid lines are the optimal fit of Eq. 11 on the average spectra. The solid (resp. dashed) lines in (d) denote the linear regression of the measured (resp. Kaimal) spectra. The red dashed line is the simplified function given by Horst (1997) where A 0 is a normalization coefficient, μ a broadness parameter and m = 3/4 for the cospectra. The spectral peak frequency f m is then parametrized as an increasing function of the atmospheric stability z/L 0 . Equation 11 is used to experimentally estimate an expression for f S xy ( f ), which we then compare to the Kaimal et al. (1972) spectra (Eq. 10).

Experimental Verification of Turbulence Cospectra Models
Virtually all high-frequency attenuation corrections are based on Moore (1986), Horst (1997), or Massman (2000, which assume that the Kaimal et al. (1972) spectra are valid (Eq. 10). This assumption does not have a notable effect on the high-frequency correction, as long as the response time of the sensors falls within the well-defined inertial subrange. For a VPEC system, the response time is of the order of ≈ 1 s. This means that the system also misses a small fraction of the flux in the lower frequency part of the turbulence spectrum, which is not necessarily well defined, as demonstrated by Smeets et al. (1998) for katabatic flow due to the influence of large-scale flow. As such, we test the validity of the Kaimal et al. (1972) spectra in a katabatic flow regime with measurements from Experiment 3. We first estimate a fixed shape parameter μ by fitting the average turbulence (co)spectra to Eq. 11. We then estimate the spectral peak frequency f m by fitting each individual (co)spectrum. The results are plotted as function of atmospheric stability in Fig. 3. The observed peak frequency shows significant scatter, but its increase with stability is well visible and roughly follows the relation derived by Kaimal et al. (1972). Furthermore, the observed averaged spectra appear wider than the Kaimal et al. (1972) spectra. Finally, our data suggest that the vertical wind speed spectra are slightly shifted to lower frequencies, although the limited sampling rate makes it difficult to fit these spectra due to aliasing. It must also be noted that the relation by Horst (1997) is a reasonable approximation for the momentum cospectra under near-neutral circumstances, i.e., z/L 0 < 0.1 (see Fig. 3).

Analytical Model
Several analytical models of Eq. 3 have been derived (Horst 1997;Massman 2000). These present an analytical function for both the normalized turbulence cospectra S xy and the sensor transfer functions T xy such that the attenuation coefficient in Eq. 3 can be integrated analytically. In the general case of two slow-response sensors, and assuming that the phase difference between x and y is small and independent of frequency, Horst (1997) writes the attenuation coefficient of the covariance x y as, where the following model for the cospectral peak frequency is used, However, both our measured cospectra during Experiment 3 and the peak frequency of the spectra derived by Kaimal et al. (1972) suggest that the following models are more accurate in the z/L 0 ∈ [0; 0.2] range, momentum flux f S uw : sensible heat flux f S wT : The modelled attenuation coefficient based on Eq. 12 and Eq. 3 using the experimentally derived expressions for τ x , τ y and f m agree within 1%. This difference is mostly due to the high-pass filter T D that is not taken into account in Eq. 12. Hereafter we will numerically calculate the integral in Eq. 3 to estimate the attenuation factor.

Summary: Model Parameters
The model described in the previous section used to estimate the attenuation of the (co)variance x y is thus entirely described by the following parameters:

Results: Evaluation of the Correction for the High-Frequency Attenuation (Experiments 1 and 2)
We apply the correction for the high-frequency attenuation to the measured VPEC fluxes. We then compare these corrected fluxes to the simultaneously measured SEC fluxes during Experiments 1 and 2. The sensible heat flux is converted to an energy flux according to H = ρ a C p w T , where the air density ρ a and air heat capacity C p are calculated using the 2-m air temperature, the 2-m air specific humidity and the surface pressure measured by the KNMI at the CESAR tower.

Accuracy of the High-Frequency Correction
The comparison of VPEC and SEC fluxes is shown in Fig. 4, where we have used the response times derived in Sect. 3 and written in Table 1. Both the corrected VPEC momentum and sensible heat fluxes show a small bias and root-mean-square error (RMSE) when compared to the SEC fluxes. The bias, which we define as the average of the difference between the two time series, is for the friction velocity equal or smaller than 0.01 m s −1 for both experiments. The RMSE value is 0.03 m s −1 and 0.05 m s −1 for Experiments 1 and 2, respectively. The sensible heat flux is slightly overestimated by the VPEC system compared to the SEC system: 6.3 W m −2 during Experiment 1 and 4.4 W m −2 during Experiment 2. This small bias is also present when comparing the sensible heat fluxes obtained with the sonic temperature and with the thermocouple attached to the SEC system. Hence the bias is not related to the vertical propeller but to the sensitivity of the thermocouples. The RMSE value of the sensible heat flux is 12.4 W m −2 and 11.9 W m −2 during the two experiments. The difference between the VPEC fluxes and SEC fluxes is of similar magnitude as the difference obtained when measuring fluxes with two adjacent sonic anemometers (Mauder and Zeeman 2018).
The only adjustable parameters in the correction for the high-frequency attenuation are the response times of the instruments. In order to quantify the sensitivity of the correction to these input parameters, the same comparison as presented above is done using perturbed response times. The optimal response parameters of the vertical anemometer, horizontal anemometer, and thermocouple are perturbed by ± A w , ± L u and ± τ T respectively. These error bars are derived from the calibration procedure (see Fig. 2 and Table 2), and result in three combinations of sensor response times that we denote τ − x , τ re f x , and τ + x . The bias and RMSE value of the corrected fluxes for each parameter set are presented in Table  2. The comparison with the SEC fluxes shows that changing the input parameters within the defined range does not affect the RMSE value of the fluxes, neither does it significantly affect the bias of the correct friction velocity. It does, however, slightly (< 1 W m −2 ) affect the bias of the corrected sensible heat flux.

Influence of the Sampling Rate
The power supply and memory usage are the most important limiting factors when considering the sampling rate of an automatic weather station, denoted F s , deployed in remote polar areas. As such we test whether the sampling rate can be decreased without significantly increasing the bias and RMSE value of the measured fluxes, as presented by Bosveld and Beljaars (2001). For this purpose, the measured VPEC time series from Experiments 1 and 2 are artificially downsampled from the initial sampling interval to sampling intervals of 1 s and 4 s. This is done by taking the first sample in each sampling interval. We then apply the processing steps and high-frequency correction from Sect. 3 to the downsampled time series, and compare the corrected and VPEC fluxes to the SEC fluxes. The resulting bias and RMSE value of the fluxes calculated with the downsampled time series are shown in Table  2.
Increasing the sampling interval from 0.1 s to 4 s increases the RMSE value of the sensible heat flux and the friction velocity by ≈ 1 W m −2 and ≈ 0.01 m s −1 , respectively. This is of the same order of magnitude as the uncertainty related to the response time of the sensors shown Table 2 Bias and root-mean-square error (RMSE) of corrected VPEC fluxes compared to SEC fluxes during Experiments 1 and 2. The experiments are described in detail in Sect. 2 and in Table 1. The bias is written in bold and the RMSE value is given in the brackets. Three different sampling rates F s of the VPEC system are compared. For each experiment, three different parameter sets used to model the instrument response times τ x have been applied, based on the uncertainty when deriving these parameters experimentally. The values for the optimal set τ re f x and the two perturbed sets τ − x and τ + x are given in Table 1 Bias (RMSE) F s 10 Hz 1 Hz 0.25 Hz Exp.1 5.9 6.3 6.6 6.1 6.4 6.8 5.9 6.3 6.6 in the previous paragraph. Furthermore, the downsampling does not significantly affect the bias, as expected from Bosveld and Beljaars (2001). We conclude from these comparisons that the overall precision and accuracy of the corrected VPEC fluxes are neither dominated by the high-frequency correction, nor by the sampling rate. After correcting for their limited dynamical response, VPEC instruments sampling at an interval of 4 s are able to measure the sensible heat flux during the first two experiments with a typical bias of 6.3 W m −2 and a RMSE value of 13.4 W m −2 (Table 2). For the friction velocity the bias is less than 0.01 m s −1 , and the RMSE value of the order of 0.05 m s −1 . Overall, this means that a corrected VPEC system measures the turbulent surface fluxes with an accuracy similar to the deviation between two adjacent sonic anemometers (Frank et al. 2015;Mauder and Zeeman 2018). Yet the main advantage of the VPEC system is its simplicity, which allows for longer periods of unattended operation in remote polar areas.

Results: One Year of Turbulent Fluxes in the Western Ablation Area of the Greenland Ice Sheet (Experiment 4)
To demonstrate its potential, we apply the correction for the high-frequency attenuation of the fluxes to one year of VPEC measurements at site S5 on the Greenland ice sheet (Experiment 4). The values for A w , L u , and τ T are found in Table 1 and are based on the calibration of identical sensors during the first two experiments.

Corrected Turbulent Fluxes
The time series of both the corrected turbulent fluxes and of the modelled attenuation coefficients are presented in Fig. 5. Over the course of one year, the attenuation factor A xy of both fluxes remains in the [0.6 − 0.8] range, which is the same interval of modelled corrections during the first two experiments. At this location, katabatic flows continuously generate turbulent mixing despite the stable stratification. This means that near-neutral conditions are nearly always observed close to the surface (i.e., z/L 0 < 0.2), thereby keeping the spectral peak to lower frequencies (Fig. 3). On the other hand, the high wind speeds shift the spectral peak to higher frequencies, but this effect is compensated by the simultaneous decrease in the response time of the propeller anemometers (see Eqs. 6 and 7 and Fig. 2). The highest flux attenuation factor of 0.5 is modelled in winter, during short periods with low wind speeds and increasing stability (Fig. 5). During such periods, the propeller anemometers become slower while the most energetic turbulent fluctuations become smaller, which results in high flux attenuation. Only one period with probable propeller freezing was identified in January (Fig.  5c1), which results in a gap of several days in the measurements. The smallest attenuation factor of 0.8 is modelled in summer, during periods with high wind speed and near-neutral stability. These are also the periods when the highest fluxes are measured: sensible heat fluxes up to −300 W m −2 and friction velocities up to 1 m s −1 .

Error Caused by the Uncertainty in Response Times
For the first two experiments, the uncertainty in the response times propagates in an uncertainty of ± 0.4 W m −2 in the corrected VPEC sensible heat flux and less than 0.01 m s −1 for the corrected friction velocity (Table 2). However the measured fluxes during the first two experiments were smaller, so these error intervals are not representative of the fluxes measured during Experiment 4. As such, the same response time perturbation exercise as in Sect. 4 is performed for Experiment 4. The different corrected fluxes show a maximal deviation of ± 10% for the sensible heat flux and ± 1.5% for the friction velocity (not shown in Fig. 5). The uncertainty of the corrected flux depends on the flux, and reaches up to 30 W m −2 for measured sensible heat fluxes of −300 W m −2 . These intervals must be interpreted as the widest interval of all possible corrected fluxes, using the response times derived in Fig. 2. It is often smaller, for instance for higher wind speeds and more neutral conditions. This maximum difference interval can only be reduced further if the response times are known with greater accuracy.

Contribution of the Sensible Heat Flux to Surface Ablation
In Fig. 5 we compare the measured sensible heat flux with the other measured components of the surface energy balance. During winter, the surface cools due to net emission of longwave radiation, which is on average compensated by the downwards sensible heat flux. The net imbalance between the two fluxes then contributes to either warming or cooling of the surface, which rarely exceeds the melting point during winter. During summer, both the net absorbed radiation by the surface and a downward sensible heat flux warm the surface. The excess energy supplied when the surface is already at the melting point results in enhanced melt.
During the first half of the melting season (May-June), the measured ice ablation closely follows R net . Then, after several consecutive warm events in July, the daily ice ablation is on average 100 W m −2 larger than R net (Fig. 5d1). This additional energy flux can for the most part be explained by the measured sensible heat flux (Fig. 5a1, d1). It must be noted that we did not measure the latent heat flux and the ground heat flux, which are also an important part of the surface energy balance. In fact during winter, important temperature gradients in the ice will compensate for the difference between the emitted radiation and the sensible heat flux. During summer, the ice is mostly at the constant melting temperature, which makes Time series of selected measurements from the VPEC system during Experiment 4 at site S5, located in the western ablation area of the Greenland ice sheet. From top to bottom: (a1) friction velocity (u * ) and sensible heat flux (H ), (b1) modelled attenuation coefficients A xy , (c1) standard deviation of the vertical wind speed σ w divided by the horizontal wind speed (U ), (d1) daily averaged ice ablation, net absorbed radiation by the surface (R net ) and sensible heat flux (H ). The panels on the right are the probability histograms for each variable calculated for the whole measurement period. Three notable periods are shaded as follows: (I) a period with a rimed or stalling vertical propeller, (II) a period in the middle of the melting season without observed ablation, and (III) an extreme warm summer event when the strongest ablation is observed. A daily average of 100 W m −2 in panel (d1) corresponds to a cumulative daily ice ablation of 2.82 cm, assuming a latent heat of melting of 334 × 10 3 J kg −1 and a constant ice density of 916 kg m −3 the ground heat flux negligible, but will enhance latent heat fluxes due to important vertical gradients in specific humidity. Furthermore, small snowfall events during summer (as is the case during event (II) in Fig. 5) will also absorb a large part of the total melt energy during the consecutive days. Finally, the footprint of the sensible heat flux is not necessarily representative of the area in direct vicinity of the weather station. Especially after several warm events, the melt water will accumulate in the surrounding narrow channels and melt ponds, which remain invisible to the ablation and radiation sensors as they are located on top of an ice hummock.

Variability of the Aerodynamic Roughness Length
When the sensible heat flux and the friction velocity are calculated using a bulk turbulence model, the aerodynamic roughness length z 0m is often unknown and becomes an adjustable parameter. Assuming Monin-Obukhov similarity, the aerodynamic roughness length is defined as the height at which the logarithmic wind profile extrapolates to zero. Here we evaluate z 0m by extrapolating the measured wind speed to the surface using the measured momentum flux at the same height, according to where we use the expression of Holtslag and De Bruin (1988) for the integrated stability functions for momentum Ψ m . We only select the measurements when z/L O < 0.2, and we assume that the last term on the right-hand side of Eq. 15 is negligible and thus set it to zero. The estimated aerodynamic roughness length is then shown as function of both time of year and wind direction in Fig. 6. The aerodynamic roughness length z 0m shows a very significant variability over the course of one year, and ranges between 10 −4 m in winter to nearly 10 −1 m in summer, which is consistent with the two-level wind profile and sonic anemometer measurements by Smeets and Van den Broeke (2008) at the same location. The reduction of z 0m from September to February (Fig. 6a) is mainly attributed to the accumulation of snow that gradually reduces the size of the ice hummocks. Furthermore, winter time sublimation smoothens the top of the hummocks, which explains that the smallest values of z 0m are observed in March and April. From May onwards, the melting of the snow increases the amplitude of the ice hummocks. The resulting rapid increase of z 0m from March to August (Fig 6b) is further enhanced by the differential melting of the ice hummocks.
The roughness length z 0m shows a remarkable dependency on the wind direction as well, with minimal values found in summer at 95 • and in winter at 115 • . The maximum values of z 0m are found for the most southerly wind directions, independently of the season. We attribute this strong directional dependency to changes in the effective obstacle area, which is most likely a direct consequence of the complex geometry and spacing of the ice hummocks (Miles et al. 2017;Fitzpatrick et al. 2019). A detailed topographical survey is required to property quantify this effect at this location.
Remarkably, the southerly wind directions are also the directions when the warmer and more turbulent air masses generate the highest sensible heat fluxes (Fig. 6). The southerly wind directions are most likely caused by the interaction between katabatic winds and barrier winds (Van den Broeke and Gallée 1996), although this interaction remains to be investigated in more detail.
The results from this experiment confirm even more that using a constant value for z 0m over snow and ice surfaces is not recommended, as shown by similar experiments at other locations (e.g., Miles et al. 2017;Vignon et al. 2017;Fitzpatrick et al. 2019). This is especially the case at this location on the western ablation area of the Greenland ice sheet, where a shift in wind direction is often associated with fast changes in air mass properties. This raises the question of the parametrization of turbulent heat fluxes in climate models, and of the possible feedback between surface roughness and surface ablation through turbulent heat exchange.

Conclusions
Motivated by the important contribution of turbulent heat fluxes to surface ablation on the Greenland ice sheet, we tested a vertical propeller eddy-covariance (VPEC) system, which is capable of continuously measuring the surface turbulent fluxes for longer periods without regular maintenance. By comparing the VPEC system to a sonic eddy-covariance (SEC) system, we found that the frequency responses of propeller anemometers and thermocouples may accurately be approximated as first-order functions, with typical response times of less than 1 s. We have shown that the resulting flux attenuation can be accurately modelled, as long as the normalized turbulence cospectra are known. Furthermore, the sampling interval can be reduced to 4 s to increase the system's autonomy even further in terms of power supply and data storage.
We presented one year of measured VPEC turbulent sensible heat and momentum fluxes at site S5 of the K-transect, located in the western ablation area of the ice sheet. Near the margins of the ice sheet, persistent density-driven katabatic flows are the main source for near-surface turbulence. Such a forcing results in quasi-continuous stable, but near-neutral, conditions. These are very favourable conditions for a VPEC system and maintain the attenuation factor of the fluxes above 0.6 (e.g., the attenuation below 40%). This long-term and continuous dataset of turbulent fluxes is invaluable for the evaluation of atmospheric numerical models, but also for the fundamental understanding of processes and drivers of surface ablation. At this location, downward sensible heat fluxes as large as 300 W m −2 have been measured, both during winter and during summer. Such values are similar or even more important than the surface net absorbed radiation, which makes them an essential component of the surface energy balance. Furthermore, we have shown that the aerodynamic roughness length is very variable in time and space, and that the highest value of nearly 10 −1 m is estimated when the sensible heat fluxes are also at their maximum.