Interactive comment on “ A study of turbulent fluxes and their measurement errors for different wind regimes over the tropical Zongo glacier ( 16 S ) during the dry season ”

Abstract. Over glaciers in the outer tropics, during the dry winter season, turbulent fluxes are an important sink of melt energy due to high sublimation rates, but measurements in stable surface layers in remote and complex terrains remain challenging. Eddy-covariance (EC) and bulk-aerodynamic (BA) methods were used to estimate surface turbulent heat fluxes of sensible (H) and latent heat (LE) in the ablation zone of the tropical Zongo Glacier, Bolivia (16° S, 5080 m a.s.l.), from 22 July to 1 September 2007. We studied the turbulent fluxes and their associated random and systematic measurement errors under the three most frequent wind regimes. For nightly, density-driven katabatic flows, and for strong downslope flows related to large-scale forcing, H generally heats the surface (i.e. is positive), while LE cools it down (i.e. is negative). On average, both fluxes exhibit similar magnitudes and cancel each other out. Most energy losses through turbulence occur for daytime upslope flows, when H is weak due to small temperature gradients and LE is strongly negative due to very dry air. Mean random errors of the BA method (6 % on net H + LE fluxes) originated mainly from large uncertainties in roughness lengths. For EC fluxes, mean random errors were due mainly to poor statistical sampling of large-scale outer-layer eddies (12 %). The BA method is highly sensitive to the method used to derive surface temperature from longwave radiation measurements and underestimates fluxes due to vertical flux divergence at low heights and nonstationarity of turbulent flow. The EC method also probably underestimates the fluxes, albeit to a lesser extent, due to underestimation of vertical wind speed and to vertical flux divergence. For both methods, when H and LE compensate each other in downslope fluxes, biases tend to cancel each other out or remain small. When the net turbulent fluxes (H + LE) are the largest in upslope flows, nonstationarity effects and underestimations of the vertical wind speed do not compensate, and surface temperature errors are important, so that large biases on H + LE are expected when using both the EC and the BA method.


Introduction
Surface turbulent heat fluxes play a significant role in the energy balance of mountain glaciers (Sicart et al., 2008;Gillett and Cullen, 2011), a crucial calculation for assessing melt and mass balance using meteorological variables.Uncertainty assessment is necessary for sensibility studies in distributed energy balance models dedicated to glacier melt estimation from meteorological variables (Conway and Cullen, 2013).Generally, energy balance studies precisely assess the radiation components of the balance, as well as the related errors, but turbulent fluxes and associated errors remain poorly understood.Commonly, turbulent fluxes are estimated with the bulk-aerodynamic (BA) method, and only a few studies have relied on direct measurements with the eddy-Published by Copernicus Publications on behalf of the European Geosciences Union.covariance (EC) method (e.g.Cullen et al., 2006;Van den Broeke et al., 2009;Conway and Cullen, 2013;Litt et al., 2015).Over mountain glaciers, several characteristics of turbulent flows challenge the measurement of turbulent fluxes, such as flux divergence between the surface and the instruments under frequent katabatic winds (Grisogono and Oerlemans, 2002), non-equilibrium surface layers related to largescale orographic disturbances (Smeets et al., 1999), temperature overestimates due to strong shortwave radiation and large surface albedo (Huwald et al., 2009), or intermittency in weak wind and strong stable stratification (Mahrt, 2007).These characteristics can cause significant errors in estimates of turbulent fluxes and thus in the energy balance.
On the Alpine Pasterze Glacier, Smeets et al. (1999) found evidence for the influence of outer-layer turbulent structures on the surface layer that impacted turbulent fluxes, and Denby and Greuell (2000) showed that under katabatic winds the BA method underestimated surface fluxes, albeit less so than when applying the profile method between two measurement levels in the air.Conway and Cullen (2013) studied the influence of uncertainties in roughness lengths, surface temperature and the choice of stability corrections on the calculation of surface energy balance on a glacier in the Southern Alps of New Zealand.Sicart et al. (2014a) studied impacts of uncertainties in roughness lengths and the illdefined zero-reference level on turbulent fluxes on the tropical Zongo Glacier in Bolivia.On the same glacier, Litt et al. (2015) compared results of the BA method to those of the EC method and attributed large underestimates of turbulent sensible heat flux (H ) by the BA method to the influence of katabatic flow oscillations or to the influence of outer-layer, large-scale structures.However, a comprehensive study of errors is necessary to interpret such a potential underestimate correctly.
Based on similarity theory, the BA method assumes that the vertical turbulent fluxes of momentum and sensible and latent heat scale with mean vertical gradients of wind speed, temperature, and specific humidity inside the surface layer (Monin and Obukhov, 1954).Divergence from similarity assumptions (Mahrt, 2007) or uncertainties in stability functions (Berkowicz and Prahm, 1982) can lead to large systematic errors.Random errors mainly result from random uncertainties in measurements or from poor estimates of roughness lengths (Smeets and Van den Broeke, 2008).With the EC method, turbulent fluxes are derived from measurements of air temperature, specific humidity and three-dimensional fluctuations in the wind speed.Systematic errors can arise from uncertainties in post-field datatreatment methods, spectral losses at high frequency (Massman, 2000), underestimating the vertical velocity component in non-orthogonal sonic anemometers (Nakai et al., 2006;Frank et al., 2013) and nonstationarity of the flow (Vickers and Mahrt, 1997).Random errors result from poor statistical sampling of the largest eddies of the flow (Mann and Lenschow, 1994;Vickers and Mahrt, 1997;Finkelstein and Sims, 2001;Hollinger and Richardson, 2005) rather than from the quite small random instrumental noise around measurements (Billesbach, 2011).Finally both methods can be systematically affected by flux divergence, advection and storage between the surface and the height of the sensors (Aubinet et al., 2012).
On high-altitude Andean tropical glaciers during the dry season, turbulent latent heat fluxes (LE) are generally an important sink of energy since they mainly consist of sublimation, which is favoured by the low atmospheric pressure.The magnitude of H can change significantly between night and day due to changes in thermal stratification (e.g.Wagnon et al., 2003;Sicart et al., 2008).Synoptic forcing generally remains weak, and thermally driven winds dominate wind circulation (Litt et al., 2015).At night and during the morning, a marked temperature inversion at low height above the surface (2-3 m) favours the development of katabatic flows (Fedorovitch and Shapiro, 2009).For weak synoptic forcing, a wind-speed maximum is frequently observed at low height.Strong synoptic forcing causes strong winds at the glacier surface, while complex outer-layer interactions with the surface layer occur (Litt et al., 2015).During afternoons of the dry season, anabatic valley winds are frequently observed (Fedorovitch and Shapiro, 2009), while the temperature inversion is less marked.Net turbulent surface heat flux (H + LE) and associated errors may strongly depend on these wind-regime characteristics.A comprehensive review and quantification of sources of uncertainties, and comparison of these errors with net turbulent fluxes under different wind regimes, is required to improve uncertainty assessment in energy balance studies.
This study is one of few concerning error assessment of turbulent fluxes measurements on glaciers (e.g.Box and Steffen, 2001;Conway and Cullen, 2013;Sicart et al., 2014a).It focuses on some of the most important errors expected in flux estimates by the BA and EC methods over the tropical Zongo Glacier.We describe measurements from a micrometeorological field campaign performed in the ablation area during the 2007 dry season from July to September and associated data analysis.Hourly mean incident and reflected shortwave and longwave radiation and 2 m temperature, wind speed and humidity were measured, as were 2 m EC data and temperature and wind-speed vertical profiles with a 6 m mast.For each of the three wind regimes observed -i.e.purekatabatic flows, strong downslope flows and anabatic flowswe derived H and LE using both the BA and EC methods and calculated their main systematic and random measurement errors.We then identified the error sources that most impacted estimates of net turbulent flux and those of secondary importance.Finally, for each wind regime, we discuss changes in its net turbulent fluxes and the influence of errors. 2 Location and data

Site and measurements
The site and field campaign are extensively described in Sicart et al. (2014a) and Litt et al. (2015).A short description is given here for consistency.The Zongo Glacier is a small (approximately 1.8 km 2 ), high-altitude (4800-6000 m a.s.l.) Andean Glacier located in the outer tropics of Bolivia (16 • S, 68 • W) in the Cordillera Real.A micro-meteorological field campaign was conducted in the ablation area at 5080 m a.s.l.from 22 July to 1 September of the 2007 dry season, during the austral winter.Temperature and wind-speed vertical profiles were measured using a 6 m high tower with 12 levels for temperature and 8 levels for wind speed (Table 1 and Fig. 1).Data were sampled every 10 s, and 15 min means were stored on a data logger.Temperature sensors were shielded and continuously mechanically aspirated.Two EC systems, measuring the three wind-speed components, temperature and humidity at 20 Hz, were installed at roughly 2 m above the ground, on two additional masts.The three masts were separated from each other by approximately 20 m and aligned perpendicularly to the glacier flow and to the main wind direction.An automatic weather station (AWS) recorded halfhourly averages of all radiation components SW inc , SW out , LW inc and LW out (SW: shortwave; LW: longwave; subscripts inc and out for incident and outgoing terms, respectively), wind speed and direction, mechanically aspirated relative humidity (RH) and air temperature.
The surface below the masts remained homogeneous and relatively smooth throughout the campaign.It was covered with snow at the beginning of the campaign.Low melt rates were observed (mean daily melt <1cm SWE).The snow became slightly rougher over time and completely disappeared from the ablation zone by the end of the campaign.Typical height changes of the snow or ice surface were 10-20 cm over horizontal distances of ∼ 10 m (Sicart et al., 2014a).The local slope on the glacial plateau around the measurement site is approximately ∼ 5 • .

Data processing
All data were split into 1 h runs.High-frequency data from the EC systems were checked for quality (Vickers and Mahrt, 1997), and low-quality runs (∼ 16 %) were discarded.A further ∼ 17 % of the remaining runs were removed because wind direction was ill-defined or wind blew through the mast structures.The remaining good-quality runs (GQRs, ∼ 70 % of all runs recorded) were despiked (Vickers and Mahrt, 1997).Planar-fit rotation (Wilczak et al., 2001) was applied in ∼ 10-day periods (between two field visits) to derive longitudinal (u), lateral (v) and vertical (w) wind-speed compo- nents.Use of constant rotation angles over a melting surface may not be advisable.Since the daily melt rates were low (measured surface height changes were 1.6 cm day −1 ) and no noticeable evolution was reported in the rotation angles calculated on a daily basis, we assumed the surface changes were sufficiently slow in order to apply the planar-fit method between two field visits, in ∼ 10-day periods.Sonic air temperature was corrected for the influence of water vapour content (Schotanus et al., 1983).Comparison of aspirated air temperature measurements from the profile mast with the temperature derived from the nearest sonic anemometer revealed residual influence from solar radiative heating, which was corrected according to Sicart et al. (2014a).Surface temperature was derived from LW out assuming a surface emissivity of 0.99 (Dozier and Warren, 1982).Since the emissivity of snow and ice can vary from 0.98 to 1.00 (Wiscombe and Warren, 1980), we accounted for this uncertainty when deriving random errors of surface temperature (see Appendix A).

Meteorological conditions and wind regimes
The austral winter dry season in the Cordillera Real in Bolivia, from May to August, is characterized by clear skies, with a mean of 20 clear-sky days per month (Sicart et al., 2014b).The air is dry at this high-altitude site since more water is needed to reach the saturation pressure.Above the ablation zone of Zongo Glacier, katabatic flows are regularly observed during the night and late morning.They are associated with a strong thermal inversion in the first few metres above the surface and a wind-speed maximum at low heights, between 2 and 3 m.During the field campaign for these conditions, we observed a mean wind speed of 1.8 m s −1 (max-imum 3.9 m s −1 ), and the surface layer (the "constant-flux" layer) was poorly defined (Sicart et al., 2014a;Litt et al., 2015).Synoptic forcing is sometimes strong during the dry season and is associated with a westerly wind in the Cordillera Real.This flow roughly aligns with that of the Zongo Glacier, generating moderate to strong downslope winds, mostly during the night.Under these conditions, we observed a mean wind speed of 3.5 m s −1 during the campaign (maximum 9.5 m s −1 ).Temperature stratification in the first few metres above the ground was also pronounced, but wind-speed maxima were not observed below 5 m, and the surface layer was more developed (Litt et al., 2015).
Anabatic flows (referred to as "upslope" herein) generally occur around midday and in the afternoon, when the surface is melting or near melting due to large daytime radiative fluxes directed towards the surface.Due to the high elevation, air above the glacier is rarely warmer than a few degrees above 0 • C, which results in a nearly neutral thermal stratification in the first few metres above the surface.During the experiment, wind speed was moderate under these conditions: a mean of 2.1 m s −1 , with a maximum of 6.8 m s −1 .
For these three wind regimes, low-frequency perturbations affect the surface-layer flow (Litt et al., 2015).The mean cospectra of vertical wind speed, w, with potential temperature, θ , and those of w with specific humidity, q, are plotted in Fig. 2 against normalized frequency n = f z/u, where f is frequency, u is the horizontal wind speed, z is the measurement height, and the overline indicates a temporal average over a 1 h run.These cospectra reveal that the lowfrequency perturbations observed affect the turbulent fluxes since they show stronger contributions at low frequency (below n = 10 −1 ) than the reference curve of Kaimal et al.  10 -3  10 -2 10 -1 10 0 10 1 x x n -5/7 0.1 0 HF peak LF peak "downslope" and "upslope" inertial subrange high frequency losses Kaimal et al. (1972) "purekatabatic" Figure 2. Cospectra of w with x = θ and of w with x = q measured in the surface layer of Zongo Glacier during the 2007 field campaign with one of the eddy-covariance systems (adapted from Litt et al., 2015).Cospectra were calculated over 1 h runs and then averaged over each wind-regime subset.Mean cospectra for the "upslope" and "downslope" subsets (dotted black line) and the lowwind-speed "pure-katabatic" subset (solid black line) are presented with the inertial subrange slope (solid green curve) and the reference Kaimal curve (Kaimal et al., 1972) (dashed red curve) for comparison.The peaks used to calculate the integral length scale of the Mann and Lenschow (1994) method are identified by the vertical blue lines; LF and HF stand for high and low frequency, respectively.The high-frequency losses assessed in Sect.4.1.2are shaded in grey.(1972), which was measured under ideal undisturbed conditions.At low wind speeds, in pure-katabatic flows, the low-frequency perturbations likely result from oscillations in the katabatic flow (McNider, 1982).When outer-layer forcing is strong, these oscillations probably originate from interactions between outer-layer eddies and the surface layer (Högström et al., 2002;McNaughton and Laubach, 2000).Under these conditions, the surface layer is probably out of equilibrium; i.e. local production of turbulent kinetic energy (TKE) does not balance dissipation since TKE transport might not be zero (e.g.Smeets et al., 1999).

Run selection according to wind regimes
We classified the GQRs into three subsets corresponding to the main wind regimes (Sect.2.3) to study characteristics of net turbulent flux for each regime during the campaign.We tried to consider as many runs as possible in each wind regime without selecting specific turbulent conditions.Hence, selection was based only on wind direction and the detection of a wind-speed maximum below the highest profile measurement.
Runs for which wind direction was downslope (between 260 and 360 • ) and for which a wind-speed maximum was observed below 5 m were classified as "pure-katabatic" (42 % of the GQRs).The GQRs for which wind blew downslope and no maximum was observed below 5 m were classified as "downslope" (26 % of the GQRs).Runs with wind direction of 45-180 • were classified as upslope (32 % of the GQRs).Throughout the text, temporal averages of a variable over the runs from one subset are noted .

Eddy-covariance fluxes
We measured the turbulent fluxes with the two EC systems, following Eqs.( 1) and ( 2), where downward (upward) fluxes were set positive (negative).The ec subscript indicates flux estimates derived from the EC method, ρ is the air density (kg m −3 ) and c p is the specific heat of humid air (J kg −1 K −1 ).The value of L e was set equal to the latent heat of sublimation of the ice (283×10 4 J kg −1 ) when the surface was below freezing and to the latent heat of vaporization (250 × 10 4 J kg −1 ) when the surface was melting.The primes denote fluctuations of the variables around their 1 h means.The Webb-Pearman-Leuning term (WPL) (Webb et al., 1980) generally remained below 2 W m −2 throughout the campaign.

Potential random errors in eddy-covariance fluxes
For the EC method, the principal source of random error is poor statistical sampling of the main transporting eddies (Vickers and Mahrt, 1997).To characterize this error, we applied the familiar and straightforward Mann and Lenschow (1994) (ML) and Hollinger and Richardson (2005) (HR) methods.Results are compared with those from studies applying other methods (i.e.Vickers et al., 2010;Salesky et al., 2012).The ML method assumes that temperature and vertical velocity fluctuations are joint-normally distributed.Under stable conditions, it expresses random error σ F in the flux F as where τ f is the integral timescale (s) of the flux, which is related to the timescale of the largest transporting eddies of the flow; P is the temporal averaging period; and r wx is the correlation coefficient between w and the scalar variable x associated with the flux.This method requires determining τ f , and several methods have been proposed to calculate it (Wyngaard, 1973;Mann and Lenschow, 1994;Finkelstein and Sims, 2001).We calculated τ f using the cospectra of w and x; the integral timescale is given by the inverse of the frequency of the maximum in the cospectra (Wyngaard, 1973).
The mean cospectra, calculated over 1 h runs and then averaged over each wind-regime subset (Fig. 2), exhibited two distinct peaks in the downslope and upslope subsets but only one in the pure-katabatic subset (Litt et al., 2015).The highfrequency peak is associated with small-scale fast turbulent eddies, whereas the low-frequency peak is associated with large-scale slow eddies.Within a 1 h run, fewer slow eddies are sampled than fast ones, which leads to larger random errors.Thus we allowed τ f to change according to wind conditions: when a katabatic wind-speed maximum was found and only a single peak was observed in heat fluxes' cospectra (Fig. 2), we used the frequency of this single peak.When no maximum was detected below 5 m and two peaks were observed in the cospectra, we used the frequency of the peak found at lowest frequency, associated with the largest eddies, which control the random sampling error.
The Vickers et al. (2010) method is based on the calculation of sub-record fluxes, i.e. fluxes computed with timescales shorter than 1 h.It evaluates the random error from the within-run variance of the sub-record fluxes.The Salesky et al. (2012) method is based on similar assumptions to the ML method.A filter of changing time constant is applied to the 1 h flux.Filtered fluxes with increasing time constant converge to the 1 h averaged flux following a power law.Salesky et al. (2012) show that the parameters of this power law can be related to the random error in the flux.Both methods have the advantage of not requiring an estimation of the integral timescale.
In contrast, the HR method assumes that two EC systems installed in separated locations -but over similar terrain, wind conditions and fetch -measure the same flux independently.This assumption holds if they are located far enough from each other so that the largest eddies do not influence both sensors at the same time.Let F 1 be hourly fluxes measured with the first EC system, F 2 be those measured by the second one, and the true value of the flux be F .If the conditions described above are fulfilled, the mean value of F 1 −F 2 during the campaign is 0. We have (i = 1 or 2): where δF i is the random measurement error in each flux, which is a random variable with mean 0 and standard deviation (SD) σ (δF i ).Since the mean And then, since the δF i are assumed to be independent, we have where δF = δF i .The SD of δF provides an estimate of the random error in the flux.This method requires collecting enough pairs of flux samples representing each reported wind condition to calculate a typical random error for each.We measured 1 √ 2 σ (F 1 − F 2 ) over equally spaced bins of wind speed over the 3 weeks of the campaign when both EC systems were available.
The effect of random instrumental noise on fluxes is expected to be small because the C-SAT3 and the LICOR7500 are accurate sensors with small measurement errors (Table 1).The ML and HR methods implicitly include this source of error (Businger, 1986;Billesbach, 2011), and we do not calculate its individual contribution.

Potential systematic errors in eddy-covariance fluxes
Non-orthogonal sonic anemometers may underestimate vertical wind velocity (Kochendorfer et al., 2012a;Nakai et al., 2006).Underestimates are larger when the wind-attack angle, relative to the horizontal plane of the sonic anemometer, increases, which can lead to underestimating fluxes by as much as 15 %.Corrections have been proposed for several sonic anemometers (Van der Molen et al., 2004;Nakai et al., 2006;Kochendorfer et al., 2012a;Nakai and Shimoyama, 2012).To evaluate the degree to which this error could affect our measurements, we selected runs during specific events for which wind conditions were representative of the three main wind regimes.For each event we compared the covariances of w with θ and q, obtained without any correction related to the attack angle, and covariances obtained after correcting the wind components for transducer shadowing following the method proposed by Horst et al. (2015).
The relative difference between these two covariances, plotted against ω, the attack angle relative to the u − v plane of the instrument, helps us to evaluate the degree of flux underestimation in each wind regime.Burns et al. (2012) also showed that sensible heat fluxes at high wind speeds might be underestimated in some EC systems (above 8 m s −1 ) due to firmware issues.This error may not be of concern on Zongo Glacier, where wind speed remains moderate (Sect.2.3), but may be of importance over other glaciers where wind speed is higher.
Flux losses at high frequency may occur due to inability of the EC system to sample the fastest, small-scale eddies because of sensor path averaging or because of sensor separation when the variable considered for the flux is measured separately from wind-speed fluctuations (Moore, 1986).This error is usually corrected using frequency-response corrections and transfer functions (e.g.Massman, 2000).The method assumes a theoretical shape for the cospectra.Since the cospectra measured on Zongo Glacier were influenced by low-frequency perturbations and did not follow commonly expected forms (Fig. 2), we did not use this method.Instead, we estimated these losses by adjusting, for the mean cospectra for each wind regime, an inertial subrange curve that follows a n −5/7 slope (Kaimal and Finnigan, 1994).The cospectra were adjusted over the frequency range in which the expected inertial subrange slope was observed, i.e. from their peak at n = 5 × 10 −1 up to the frequency at which the cospectra deviated from the inertial subrange slope (Fig. 2).Between this frequency and the high-frequency end at n = 10 1 , deviations from the n −5/7 slope were attributed to high-frequency losses (Fig. 2).
Insufficient averaging time can lead to incomplete sampling of the largest eddies and to flux losses at low frequency and thus to systematic underestimates of flux.Nevertheless, long averaging times would include spurious fluxes due to changes in the state of the turbulent flow under the effect of changes in external forcing, i.e. nonstationarity.A good choice for the averaging timescale is a trade-off between these two constraints.To check if the 1 h averaging timescale was appropriate, we calculated multi-resolution decomposition (MRD; Vickers and Mahrt, 2003) cospectra of w with θ and of w with q.This decomposition isolates the contributions to H and LE coming from increasing timescales as if the fluxes were estimated using increasing averaging periods.The sampling must be long enough so that the MRD cospectra fall to zero at the longest timescales considered.Another application of this method is identification of a gap scale between the timescales of fast turbulent eddies and the nonstationary fluctuations induced by the slowest motions.

Bulk-aerodynamic method
Sensible heat fluxes H and latent heat fluxes LE were estimated by applying the BA profile method between the surface and each of the eight measurement levels of wind speed and the corresponding temperature levels from the profile mast (Table 1).The subscript ba on H and LE refers to fluxes estimated by this method.Specific humidity was calculated from relative humidity measured at the AWS and ventilated temperature profiles (see Appendix B).We have LE ba,j = ρL e k 2 u(q j − q s ) ln with the same sign convention as for the EC method.The subscripts j and s refer to the j th wind-speed measurement level and the surface level, respectively.The von Karman constant k was set to 0.4.The length scale L * is the Obukhov length, defined as with u * being the friction velocity and the subscript v referring to virtual temperature.Stability corrections for momentum, sensible heat and humidity transfers (respectively, ψ m , ψ h and ψ q ) were taken from Brutsaert (1982).The different forms of these corrections are given by the following equations, where the variable x is defined as x = (1−16z/L * ) 1/4 .If z/L * <0, we have and if z/L * > 1, we have The length scales z 0 , z t and z q are the momentum, thermal and humidity roughness lengths (m), respectively.They were derived from least-square iterative fitting of the mean temperature and wind-speed vertical profiles (Sicart et al., 2014a) and assuming z q = z t (Andreas, 2002).We use their median values, z 0 = 2.07 × 10 −3 m and z t = 0.09 × 10 −3 m, since they did not vary significantly during the campaign.
The bulk formulation is designed for surfaces normal to the gravity vector and might not be suitable for sloping surfaces under the influence of slope flows (e.g.Denby, 1999;Grachev et al., 2015).In that case, the formulation of the Obukhov length in Eq. ( 9) may have to be modified to include scalar fluxes tilted to the gravity vector (Horst et al., 1988;Grachev et al., 2015).We did not attempt to derive such a formulation in the present study.Rather, we resorted to quantifying the bias of turbulent flux measurements over glaciers resulting from the use of the bulk method in its conventional, horizontal-surface form, which is typical of most glacier energy balance studies.

Potential random errors in the bulk fluxes
We estimated errors of the fluxes derived from the BA method arising from random noise around the measurements both analytically and with Monte Carlo simulations.For error propagation in the functions H and LE = f ( u, (θ, q), z, z 0 , z t,q ) (Eqs.7 and 8), we used a local linearized model, assuming that the errors were independent, normally distributed, and smaller than the partial derivatives (δ are expected measurement errors, and are differences between the measurement height and the surface): Detailed calculation of each term on the right-hand side is found in Appendix C. We assumed the stability corrections were constant under small variations δu, δ θ , δ q and δz, which greatly simplifies the calculations.The effects of this simplification are discussed in Sect.4.2.1.Values of δu, δ θ, δ q and δz were adjusted according to the expected errors for wind speed, air temperature, relative humidity and height, respectively (Table 1).The manufacturer provides a value of δu of 1 % of the reading, but when wind speed is low the error may not tend to zero, so we set it to 0.1 m s −1 .The error in θ was set to 0.35 • C, larger than the expected thermocouple random noise (σ θ = 0.1 • C, Table 1), since θ was affected by surface temperature errors (see Appendix A).The random error in q depends weakly on the random error in air temperature and strongly on the error in relative humidity (see Appendix B).We set the random error in q to 3 %, equal to the error in the relative humidity.The sounding height ranger is quite accurate (±0.01 m), but the surface was irregular at a local scale, so its error was set to ±10 cm.We used δz 0 and δz t values derived from the error analysis of Sicart et al. (2014a): δ ln z 0 = δ ln z t = 1.5.
The Monte Carlo-based error estimate was obtained, for each measurement, by simulating 1000 dispersed measurements with a Gaussian random function whose mean equaled the measurement and whose SD equaled the expected random instrumental error previously mentioned.From these dispersed data, we recalculated 1000 estimates of turbulent fluxes for each run and profile level.Random errors in mean fluxes of a subset were calculated as the interquartile range of the 1000 mean fluxes obtained from the simulated measurements.

Potential systematic error in the bulk fluxes
Systematic errors in the BA method could originate from systematic biases in the measurements.Air temperature can be overestimated due to heating of the shelters by shortwave solar radiation (Huwald et al., 2009).This bias was consid-ered negligible, however, since the air-temperature shelters were continuously artificially aspirated, and residual deviations during clear-sky days were corrected using air temperature derived from the EC systems (Sicart et al., 2014a).Surface temperature measurements, derived from measurements of the longwave emission of the snow or ice surface with pyrgeometers, can be overestimated.This is due to the penetration of shortwave radiation through instrument filters and heat transfers due to heating by shortwave radiation of instrument domes (Yamanouchi and Kawaguchi, 1984).This effect was evident during the day since estimates of surface temperature reached 1.5 • K above the melting point.We tested the impact of this bias on flux estimates by using two methods to correct it.The first method set the surface temperature to 0 • C when the raw estimate based on outgoing longwave radiation exceeded 0 • C (Sicart et al., 2008;Wagnon et al., 2003); we called this method "T s -blocking".It gave an upper boundary for surface temperature.The second method was inspired by Obleitner and De Wolde (1999); we corrected the outgoing longwave emission by a percentage of SW inc .We chose to use a percentage of 0.6 %, half the value of Obleitner and De Wolde (1999), since it yielded a surface temperature near 0 • C when field observations reported snow or ice melt.We called this second method "T s -corrected".
The error related to the ill-defined zero reference level was studied by Sicart et al. (2014a).They showed that this effect had a weak influence on final turbulent heat fluxes since the ablation zone of Zongo Glacier remained quite smooth during the campaign.
We also studied systematic errors related to unmet theoretical requirements that could affect the BA method.First, advection or divergence of fluxes below the sensors could be high when the surface layer was not well developed, especially when a katabatic wind-speed maximum occurred at a low height under low-wind-speed conditions (Litt et al., 2015).This probably led to underestimating surface fluxes when they were assessed too high above the ground, which also affected the EC method.Estimating fluxes with the BA method, applied with increasing measurement heights on the profile mast, helped us to document this issue.Second, low-frequency perturbations affected the surface-layer flow (Sect.2.3).The additional TKE induced intermittently by low-frequency perturbations probably did not scale with mean gradients in the surface layer, potentially causing underestimates of flux (Litt et al., 2015).Comparisons with the EC-based fluxes were used to document this bias.Finally, use of inadequate stability functions could also induce biases.The most commonly used functions lead to underestimating fluxes since they do not account for intermittency under strong stable conditions (e.g.Andreas, 2002;Mahrt, 2007).These cases are related to small fluxes under weak turbulent mixing, of secondary interest for surface energy balance studies.Furthermore, this bias is expected to be small on Zongo Glacier since the mean value of the bulk Richardson number (Ri b ) (Stull, 1988) between the first and fourth (d) Cospectra between the w q time series of both EC systems, using the same notation as for w T .
profile levels was 0.07 (Sicart et al., 2014a).For this low stability, the formulas do not differ much from each other, and intermittency related to weak wind and strongly stable stratification is not a concern.

Random error calculations
The random error in H ec and LE ec increased with wind speed for the Mann and Lenschow (ML) method and did not change significantly with wind speed for the Hollinger and Richardson (HR) method (Fig. 3a and b).Both methods provided similar estimates for low wind speeds below 2-3 m s −1 , which represented the majority of runs (u<2 m s −1 58 % of time; u<3 m s −1 88 % of time).For higher wind speeds, the random error predicted with the ML method increased and was larger than that of the HR method, which remained constant.Large random errors were probably due to poor sampling of large-scale outer-layer structures, which were frequently observed in the downslope and upslope sub-sets when wind speed was high (Sect.2.3).We chose to use the highest error derived from the ML method in the rest of this study.In relative terms, for hourly fluxes, the random errors in H ec (LE ec ) derived from the ML method were between 22 and 60 % (between 36 and 98 %).The same orders of magnitudes but slightly lower values were derived by Vickers et al. (2010) with their method, during the night over a maize field (19 % for H and 23 % for LE).With their filtering method Salesky et al. (2012) found much lower values (10 %) for the random error on H for data collected in California.Interpreting these differences is not straightforward since the context of Vickers et al. (2010) and Salesky et al. (2012) was very different from ours.This comparison suggests that our method tends to maximize the errors.
The HR method was probably not adapted for estimating errors using the EC-mast configuration available, considering the flow characteristics observed.In upslope and downslope subsets, outer-layer eddies that interacted with the surface layer moved through the sensors in 50-100 s at wind speeds of about 3 m s −1 (Litt et al., 2015), which indicated that they were larger than the distance between the two EC systems (roughly 20 m, Fig. 1).Under these conditions, fluxes measured by the two EC systems may not be independent since Table 2. Overview of turbulent fluxes and related random errors calculated with the eddy-covariance (EC) method and bulk-aerodynamic (BA) method using the fifth level of the profile (∼ 2 m height).Time occupied by the subset in the campaign (first column), net turbulent exchange (H +LE) for each wind regime and its ratio (%) to total net turbulent exchange of good-quality runs (GQRs) (second column), share of H and LE in the net exchange and mean flux exchange over the GQRs (third and fourth columns).Relative random errors calculated with the Mann and Lenschow (1994) (for EC fluxes) and Monte Carlo (for BA fluxes) methods in individual H and LE fluxes and net turbulent fluxes (fifth to seventh columns).The two last lines show results obtained for all GQRs., where the subscripts "sub" and "all" refer to averages over the subset considered and to averages over all available runs during the campaign, respectively.The percentages after the ± sign indicate the SD of, or uncertainty in, this ratio.b The share of turbulent energy due to individual H and LE fluxes was calculated as they originated from the same eddies.This is confirmed by analysing the spectral dependence of the covariance of fluxes between the two EC systems (Fig. 3c and d).For low wind speeds in the pure-katabatic subset, the covariance of fluxes from both EC systems was near zero at all frequencies.For high wind speeds in the upslope and downslope subsets, significant covariance was found at low frequency (n<10 −2 ), whereas near-zero covariance was found at high frequency (except for the small sensible heat fluxes in the upslope subset).Fluxes derived from each EC system were independent at short timescales but were not so at the timescale of large eddies.This may explain why HR error followed ML error at low wind speeds but was smaller at high wind speeds.The ML method estimated the errors induced by poor sampling of the large outer-layer structures by using an adapted τ f , whereas the HR method could not estimate it correctly.
On average, the relative random error was 34 % (Table 2) of net EC fluxes in the pure-katabatic subset.Relative random error was the largest (60 %) for the downslope subset because net fluxes were weak for this regime, and high wind speed was associated with large random errors (Fig. 3a and b) due to the presence of outer-layer structures.Relative random error was the lowest (8 %) for the upslope subset because net fluxes were the largest for this regime.Over the campaign, relative random errors derived from the ML method cancelled out to a mean value of 12 %.

Estimates of systematic errors
Relative difference between covariances corrected with the Horst et al. (2015) method and uncorrected covariances for one of the EC systems are shown in Fig. 4, for selected cases whose characteristics were representative of the three windregime subsets.The difference is shown for different attack angles.For most attack angles, corrected fluxes are between 3 and 6 % larger than the uncorrected fluxes, in agreement with findings from Horst et al. (2015).Corrections larger than 6 % are found only for the pure-katabatic case, for both H and LE between −10 • and 0 • attack angles, and for LE in downslope cases above 25 • attack angles.On average over all attack angles, the correction remains small; it is 5.6, 3.7 and 3.6 % for H in the pure-katabatic, upslope and downslope cases, respectively, and 5.7, 3.9 and 2.6 % for LE, respectively.Consequently, we assume that roughly both sensible heat and latent heat fluxes were underestimated by about 6 %, as an upper boundary.
High-frequency losses calculated from the mean cospectra of w with θ and with q were slightly more pronounced for latent than for sensible heat fluxes, probably because latent heat calculations were affected by sensor separation between the C-SAT3 and the LICOR7500 (about 0.30 m).We also checked the ratio between the flux at the shortest timescales of the MRD cospectra (i.e.estimated from two instantaneous measurements at 20 Hz) and the flux at the observed gap scale.If this ratio is low, high-frequency losses must also be low (Reba et al., 2009).MRD cospectral ratios and losses calculated with Fourier cospectra remained small Relative difference between uncorrected covariances (w x , x = θ, q) and covariances corrected for transducer shadowing (w x c ) for different attack angles and selected events inside each of the main wind-regime subsets.Results for (a) w θ and (b) w q during "purekatabatic" (black), "downslope" (blue) and "upslope" (red) events are shown.for the three subsets (∼ 3-4 %, Table 3), showing that highfrequency losses were probably not too large.
The MRD cospectra of w with θ and with q for the three subsets are shown in Fig. 5.In the pure-katabatic subset, no clear separation appeared between short-timescale and longtimescale contributions to fluxes.Starting from between 10 0 and 10 1 s, towards longer timescales, contributions to fluxes evolved erratically around zero (Fig. 5a and b).At the longest timescales, their median contribution fell to zero, and the interquartile range of the dispersion was comparable to the magnitude of the peak at the short timescale (around 10 0 s).The 1 h sampling time was probably too short to capture significant flux at the longer timescales, but on average over several runs this led to a random error rather than to a systematic bias.In the downslope subset, cospectra fell asymptotically to zero towards the longest timescales, as a median over all the runs, and so did dispersion (Fig. 5c and d).The fluxes seemed correctly sampled with a 1 h averaging time.In the upslope wind regime, w and q cospectra were similar to those of the downslope subset; thus, latent heat flux seemed correctly sampled, but sensible heat flux remained weak, and the cospectra of w with θ were small and erratically dispersed around zero (Fig. 5e and f).
These results suggest that systematic errors in EC fluxes led to underestimating the fluxes.The main source of underestimation was probably related to potential underestimates of w in non-orthogonal sonic-anemometers (∼ 6 %, Fig. 4).Putting aside flux divergence between the surface and the sensors and including the high-frequency spectral losses (∼ 4 %, Table 3), underestimates from the EC method could be as large as 10 %.

Random error calculations
At a height of 2 m, the fifth level on the profile mast, the largest contribution to random errors came from roughness length uncertainties since they were poorly known (Fig. 6a).This was also observed for the other measurement heights (not shown).The second-largest random error arose from temperature uncertainties.The errors in u and q were of secondary importance for the fluxes.The random error resulting from height uncertainty was the lowest (green curve, barely visible in Fig. 6a).Analytical and Monte Carlo-based error-calculation methods produced similar results (Fig. 6b).The Monte Carlo method predicted slightly larger uncertainties, probably because it accounted for variations in stability functions due to measurement errors, which were ignored in the analytical calculations.In relative terms (not shown), the difference was <5 %.We used Monte Carlo analysis in the rest of this study since it provided an upper boundary to the error.Random errors decreased with increasing measurement height (Fig. 7), probably because random measurement errors (δu, δθ, δq, δz or δz 0,t,q ) had similar magnitudes, whereas differences ( ) between air variables and surface variables increased.
At the fifth level (2 m height), total relative random error, calculated as the interquartile range of all H + LE estimates M. Litt et al.: Turbulent fluxes over the Figure 5. Multi-resolution decomposition (MRD) cospectra (black lines) of w with θ (upper panels, a, c and e) and of w with q (lower panels b, d and f) from individual runs of "pure-katabatic" (left, a and b), "downslope" (centre, c and d) and "upslope" (right, e and f) wind-regime subsets.Quantiles 25, 50 and 75 of the cospectra at each dyadic timescale calculated over a subset (red curves) and the zero line (blue line) are plotted.Figure 6.Results of random-error calculations using the bulk-aerodynamic method for the fifth level (∼ 2.08 m) of the profile mast.(a) Change in individual random-error terms during the campaign estimated with the analytical method: z 0 error (black), z t,q error (orange), θ error (red), u and q error (blue), and z error (green).(b) Change in the random error in net turbulent flux during the campaign, calculated with the Monte Carlo (black) and analytical (red) methods.obtained from Monte Carlo runs divided by their median, was sometimes as large as 30-40 % ( 20 W m −2 in absolute value, Fig. 6) of net fluxes at hourly time steps.Random errors in repeated measurements of the same variable tend to cancel out; thus, cumulative random errors for the entire campaign were reduced to ∼ 6 % (Table 2).Since net turbulent fluxes were small for the downslope subset and random errors were high, no clear cancellation of the random errors was observed for this subset, and they remained relatively high on average (39 %, Fig. 7 and Table 2).Similarly, the median relative random error was 35 % for the pure-katabatic subset.It was only 4 % for the upslope subset since net turbulent fluxes were high for this subset.

Estimates of systematic errors
In the pure-katabatic subset, the magnitude of BA-based H , LE and their sum decreased as measurement height increased (Fig. 7).At a 2.0 m height (fifth level), sensible (latent) heat flux was 90 % (62 %) of that obtained at 0.70 m (second level).Higher up, at 5 m (eighth level), this flux was , of raw estimates of surface temperature (black), surface temperature derived from the "T sblocking" method (red), surface temperature derived from the "T scorrected" method (blue) and incoming shortwave radiation (green).
only 17 % (16 %), probably due to the shallow depth of the surface layer when a wind-speed maximum was observed at low height, around 2.0 m (Sicart et al., 2014a;Litt et al., 2015).Vertical divergence was also observed in the downslope and upslope subsets but was significant only at heights greater than 2.0 m.We found similar latent and sensible heat fluxes at 2.0 and 0.7 m, but those at 5.0 m were only 70-80 % of those at 0.7 m.When no wind-speed maximum was observed, the surface layer was probably slightly thicker than that in the pure-katabatic regime, due to increased mixing induced by higher wind speeds (Sicart et al., 2014a;Litt et al., 2015), but its extent above the ground probably remained low.Nevertheless, our observations show that measurements at a height of 2 m provide reliable estimates of surface fluxes in downslope and upslope wind regimes.For all wind regimes, the relative thinness of the surface layer probably also led the EC method to underestimate surface fluxes.
During the day and when the surface was not melting, surface temperature from the T s -blocking method was often slightly warmer than that from the T s -corrected method (Fig. 8).In general, the T s -blocking method predicted greater losses of turbulent energy (Fig. 9b).Averaged over the campaign, the magnitude of net turbulent fluxes estimated with the T s -corrected surface temperature was 56 % ( 4 W m −2 ) of that estimated with the T s -blocking surface temperature.This percentage was 64 % ( 10 W m −2 ) in the upslope subset.Because net turbulent fluxes were low in the purekatabatic and downslope regimes, the relative difference was large, but less than 2 W m −2 in absolute value.
A warmer surface reduces the temperature gradient and thus the sensible heat flux, and it also increases q s and thus q, leading to higher sublimation rates.When shortwave radiation influenced longwave radiation measurements, but actual surface temperature was below 0 • C, the T s -blocking method probably overestimated surface temperature because it applied no corrections in such a case (Fig. 8).When uncorrected radiative temperature was above 0 • C, but the surface was not melting, the T s -blocking method set the surface temperature to 0 • C (Fig. 8), also leading to overestimation of the temperature.In contrast, the T s -corrected method is physically based, and the correction was tuned to obtain a surface temperature of 0 • C when we observed melt during field trips.In this sense, the T s -corrected method seems more accurate than the T s -blocking method.Still, the former may still be inaccurate and lead to uncertainties in surface temperature.This shows that determining surface temperature is a critical point when estimating turbulent fluxes with the BA method since it has a strong influence on flux estimates.Since the T s -blocking method provided an upper boundary to surface temperature, it also provided an upper boundary to the magnitude of the net turbulent fluxes derived from the BA method.
At a 2 m height, the BA predicted considerably smaller flux magnitudes than the EC method for all subsets (Fig. 7, grey circles), which suggests that a strong bias affected one or both methods.Correcting for surface-temperature biases in the BA method, likely leading to lower net fluxes in magnitude, and correcting the EC method for its underestimates would lead to an even larger discrepancy between the fluxes from both methods.Random errors in the BA-based and EC-based fluxes were too small to explain this difference (Fig. 7), except for downslope net fluxes.Horizontal variability in fluxes between the masts was probably negligible since the surface remained homogeneous and the fluxes provided by both EC systems were similar (mean differences <1 W m −2 ).The most reliable explanation is that the interaction of outer-layer coherent eddies with the surface-layer flow or oscillations in the katabatic flow induced turbulent mixing that did not scale with local mean gradients, and thus that the BA method did not account for a portion of the fluxes (Litt et al., 2015).Similar studies in nonstationary turbulent flows showed that stability functions may be overestimated, leading to systematic underestimates of fluxes (Cheng et al., 2005;Mahrt, 2007).Alternative formulations for Eqs. ( 7) and ( 8) that account for this influence cannot be derived eas-ily because the assumptions of similarity theory no longer apply to the flow in these situations.Underestimates of net turbulent fluxes were significant only for the upslope subset ( H + LE ba,5 was only 16 W m −2 when H + LE ec was ∼ 32 Wm −2 , Fig. 7).

Net turbulent fluxes and wind regimes
During the 2007 campaign on Zongo Glacier, sublimation was high (−19 to −34 W m −2 , Fig. 9a and Table 2) because the air is very dry at high elevation due to its cold temperature and low density.Sensible heat flux was generally opposite in sign to latent heat flux but, on average, lower in magnitude (13-24 W m −2 Table 2).Net turbulent fluxes resulted in a loss of energy for the glacier (Fig. 9b).
During the night, for downslope and pure-katabatic subsets, temperature inversion in the first few metres above the surface favoured downward sensible heat fluxes.When a katabatic wind-speed maximum was observed at low height, sensible and latent heat fluxes were weak, the closest to zero of all subsets (∼ 10 W m −2 at the fifth level, Figs.7 and 9a), mainly because wind speed was low.In downslope flows, to strong winds, both sensible and latent heat fluxes had large magnitudes (30-55 W m −2 at a 2 m height, Figs.7 and 9a).Net turbulent flux was generally near zero in these two cases because the sensible heat flux was opposed in sign and nearly equal to flux (Table 2).Systematic underestimation of fluxes by the BA method, due to the influence of low-frequency perturbations, was similar in magnitude for H and LE, and partly cancelled out in net turbulent flux.The BA method predicted near-zero net fluxes, while the EC method predicted small losses in the downslope subset and small gains in the pure-katabatic subset.The effect of surface-temperature biases on BA fluxes probably remained negligible because these regimes mainly occurred during the night, when solar radiation was zero.Underestimates of flux by the EC method due to potential underestimates of vertical wind speed probably cancelled out in net fluxes since we considered that they were of the same magnitude for H and LE.
These two types of downslope flows were dominant (68 % of the time) but of secondary importance in net turbulent flux exchange during the campaign: ratios of ∼ 15 % between the mean net turbulent flux over these subsets and the mean over the campaign are found with the EC method (Table 2).Since net turbulent fluxes were small, one may argue that the significant divergence of fluxes (pure-katabatic subset) with height and the large relative random errors (34 % in the purekatabatic and 60 % in the downslope subset) were probably not a concern.However, the high uncertainty in flux leads to non-negligible uncertainties in the contribution of these regimes to net turbulent exchange over the campaign, which remains ill-defined, especially for the downslope regime (Table 2).
For the daytime upslope subset, sensible heat flux was low (a few watts per square metre) because stratification in the first few metres above the surface was near-neutral, but latent heat losses remained high (−25 W m −2 , Figs. 7 and 9a; Table 2) since humidity gradients remained large; thus, net turbulent flux was largely negative (Fig. 9b).Although this regime was observed only 32 % of the time, it contributed significantly to net turbulent exchange during the campaign: ratios of 86 and 100 % between mean net turbulent flux over the subset and that over the campaign are found for the BA and EC methods, respectively (Table 2).Since the net flux magnitude was high, relative random error derived from the ML method remained moderate (8 %), and relative random error in net BA fluxes was low (4 %).Errors in BA-method fluxes due to surface-temperature uncertainties were probably high if the surface was not melting during the daytime regime (Fig. 9b).Systematic errors due to the influence of low-frequency perturbations did not have the same magnitude in the small H as in the significantly negative LE, which led to large systematic bias in net turbulent fluxes in the upslope subset.Underestimates of fluxes H and LE due to underestimating vertical wind speed probably did not compensate each other in net turbulent fluxes, because H and LE were not of the same magnitude, and this error should be considered with caution.The divergence with height was moderate (8 % loss in net fluxes at 2 m, 20 % at 5 m), which suggests that this error is of secondary importance for flux measurements at a 2 m height in upslope flows.
Overall, error analysis shows that the BA method severely underestimated the magnitude of net turbulent fluxes.On average over the campaign, they were only 60 % of net EC fluxes (Table 2), an underestimation most likely due to the inability of the BA method to account for the flux induced by katabatic oscillations or outside-layer interactions with the surface layer.These influences lead to nonstationarity of the flow and flux divergence above the ground and thus to divergence from the required conditions for similarity to hold, which may explain the observed flux underestimations.

Conclusions
We calculated turbulent sensible (H ) and latent (LE) surface heat fluxes during a micrometeorological field campaign deployed at 5080 m a.s.l. in the ablation zone of the tropical Zongo Glacier during 1 month of the austral winter of 2007.Both eddy-covariance (EC) and bulk-aerodynamic (BA) methods were applied.We calculated the related random errors in each flux and qualitatively estimated the main systematic errors.We studied the importance for total turbulent energy transfer during the campaign of the three dominant wind regimes: weak katabatic flows with a wind-speed maximum at low height (∼ 2 m), strong downslope flows without a wind-speed maximum at low height, and moderate daytime upslope flows.We finally studied the influence of errors in total net turbulent fluxes.
In general, fluxes H and LE were a gain and a loss, respectively, of energy for the glacier.Whereas turbulent fluxes had high magnitudes under strong downslope flows, this regime was of lesser importance for net turbulent flux because H and LE were of the same magnitude and canceled out.The katabatic regime also had small fluxes that cancelled out.The highest flux losses occurred during upslope flows because fluxes H were low during the day due to a small temperature difference between the air and the surface, whereas sublimation (LE<0) remained high due to large humidity gradients (Sicart et al., 2005).
For moderate-to-high-wind-speed conditions (>3 m s −1 ) in upslope and strong downslope flows, large outer-layer eddies interacted with the surface flow, and large random errors in the EC method originated from poor statistical sampling of these structures.The random error was moderate when wind speed was low.We showed that the Hollinger and Richardson (2005) method, when applied with the mast configuration of the 2007 campaign, was not adapted for deriving sampling errors due to the presence of large-scale eddy structures because the masts were too close to each other (∼ 20 m).Instead, we used the Mann and Lenschow (1994) method.
When a wind-speed maximum was observed at low height, nonstationarity affected individual fluxes erratically and led to additional random errors in mean EC fluxes.Mean relative random error was 12 % over the campaign.Biases could lead to underestimating flux magnitude by around using the EC method.The largest bias was potential underestimation of vertical wind speed by the non-orthogonal anemometers, which could lead to underestimating fluxes' magnitudes by about 6 %.High-frequency losses remained small, at most ∼ 4%.
For the BA method, the main source of random errors came from uncertainties in roughness lengths; they were large at the hourly timescale (∼ 40 %) but decreased to 5 % for net turbulent fluxes estimated for the entire campaign.Systematic errors were generally high.Surface-temperature errors induced by solar radiation effects probably led to overestimating energy losses in turbulent fluxes, especially in daytime upslope flows.A larger bias arose, however, from nonstationarity induced by interactions of low-frequency perturbations with the surface layer and dominated overall systematic errors.This led to net underestimation of flux magnitude: using the BA method, net turbulent flux over the entire campaign was only 60 % of that measured with the EC method.
Apart from these effects, vertical divergence of fluxes probably occurred due to the shallow depth of the surface layer, which affected both EC and BA methods.Flux divergence had a limited effect below 2 m (<15 %) in strong downslope flows or upslope flows.It was large when a katabatic maximum was observed ( 67 % at 2 m).Nonetheless, lowering the measurement height to reduce this effect would lead to larger uncertainties in BA-based fluxes since it would result in an increase in measurement random error.The results presented here suggest that a height of 1 m for estimating surface turbulence fluxes with the BA method could be a convenient trade-off between these two constraints.
In the context of energy-balance studies on Zongo Glacier, random errors in the BA and EC methods would be large when studying melt processes at short timescales since they were occasionally large for short periods.But random errors would probably be of secondary importance when calculating melt over monthly timescales since random errors cancel out to moderate values.Nonetheless, we showed that the contribution of strong downslope flows to net turbulent exchange was not well defined because of large random errors affecting small net fluxes.In night-time downslope flows, most systematic errors were low or cancelled out, together with the fluxes.Only vertical flux divergence, affecting both BA and EC methods, did not cancel out.In daytime upslope flows, systematic biases in largely negative LE fluxes and small H fluxes did not cancel out.Flux divergence remained moderate but cannot be neglected because the net fluxes were large.Conversely, because of the large net fluxes, relative random errors were lower in these cases.Since this regime exhibited the largest net turbulent exchanges, these issues need further investigation.
The measurements and calculation methods for estimating surface temperature need improvement.Vertical divergence of flux must be studied in more detail; contributions from modelling would be useful, but models are still difficult to apply to katabatic flows disturbed by outer-layer interactions (e.g.Denby and Greuell, 2000).Similar studies are necessary for other glaciers, influenced by different climates, where the contribution of turbulent sensible and latent heat fluxes to surface energy balance differs from that on tropical glaciers, e.g. at high latitudes and on low-altitude glaciers, where sensible heat fluxes are large but sublimation is low and deposition can be high.The role of the errors studied here might differ considerably at sites with distinctly different meteorological characteristics, and the influence of turbulent flux uncertainties on surface energy balance might also differ.An estimation of the turbulent fluxes could be obtained from the estimation of melt rates and of all other energy balance terms.Such an independent estimation may help in understanding the origin of biases evidenced herein.Since ice or snow temperature measurements are challenging due to solar radiation contamination heating of the sensors below the surface (Helgason and Pomeroy, 2012), the potentially large conduction flux below the surface remained unknown during the 2007 campaign on Zongo Glacier.This prevented us from conducting an energy balance closure check.Energy balance closure studies should be conducted on temperate glaciers when the surface is constantly melting, the temperature profile below the surface is isothermal at 0 • C and conduction fluxes are zero.Such studies may improve the estimation of measurement errors in turbulent heat fluxes.z 0,t,q (Eqs.15): (δz q ) 2 , (C2) where we assumed errors in the air density (thus on the air pressure measurements) were negligible.Also, the specific heat and latent heat coefficients are supposed to be accurately known.This expression can be rewritten in terms of relative error: where the terms related to each relative error for each type of measurement appear independently.These terms can thus be studied separately before incorporating them back into the expression for total error.

C2 Wind-speed error
The relative error in wind speed propagates as follows: The error in height measurements likewise affects latent heat flux (changing roughness lengths and stability functions to those for humidity).The equation shows this error decreases as measurement height increases.

C5 Roughness lengths
Uncertainties in roughness lengths measurements propagate as follows: This formula is identical for thermal and humidity roughness lengths in the case of latent heat flux.

C6 Specific humidity error
Random errors in humidity measurements affect latent heat flux as follows:  The second term on the right-hand side of the equation expresses the error in q due to errors in temperature measurements.Due to the quadratic term in the denominator, this term remains smaller than the relative error in relativehumidity measurements.

Figure 1 .
Figure 1.Overview of Zongo Glacier and the instrumental set-up deployed during the 2007 dry season campaign.(a) Picture of the glacier, taken from the moraine at 5080 m a.s.l.The red circle indicates the location of the measurement site.(b) Picture of the profile mast and the eddy-covariance (EC) systems.The automatic weather station, located behind the photographer, is not visible.

Figure 3 .
Figure 3. Random error in fluxes derived from the eddy-covariance (EC) method and covariance of fluxes between the two EC systems.(a) Change in the random error (W m −2 ) in sensible heat flux H over equally spaced intervals of wind speed (width: 1 m s −1 ).Results are from the Hollinger and Richardson (HR) method (blue line and circles) and the Mann and Lenschow (ML) method (red line and circles).Small red circles indicate estimates of the ML method for individual runs.(b) Results of the the HR and ML methods applied to latent heat LE estimates.(c) Cospectra between the w T time series of both EC systems.Results are for downslope (blue), "upslope" (red) and "purekatabatic" (black) subsets.Circles indicate medians of cospectral values over equally spaced logarithmic intervals of normalized frequency.(d)Cospectra between the w q time series of both EC systems, using the same notation as for w T .
Figure 4. Relative difference between uncorrected covariances (w x , x = θ, q) and covariances corrected for transducer shadowing (w x c ) for different attack angles and selected events inside each of the main wind-regime subsets.Results for (a) w θ and (b) w q during "purekatabatic" (black), "downslope" (blue) and "upslope" (red) events are shown.

Figure 7 .Figure 8 .
Figure 7. Mean turbulent flux exchange in the "upslope" (red), "pure-katabatic" (black) and "downslope" (blue) wind regimes estimated with the bulk-aerodynamic method and different profile mast levels, plotted against the height of measurement.Results show dispersion from Monte Carlo simulations.The middle line in box plots indicates the median, the boundaries of the boxes the quantiles 25 and 75 and the whiskers the quantiles 5 and 95.Box plots of eddy-covariance (EC) flux influenced by the error estimated using the Mann and Lenschow method are highlighted in grey circles.Results are for (a) latent heat, (b) sensible heat and (c) net turbulent flux.

Figure 9 .
Figure 9. Turbulent fluxes derived from the for the fifth measurement level on the profile mast (∼ 2.08 m) with different surface-temperature derivation methods over the entire campaign.(a) Time series of fluxes H and LE estimated from surface temperature derived from LW out measurements corrected with the "T s -blocking" method.(b) Net turbulent flux estimated with surface temperature corrected with "T s -blocking" (red) and "T s -corrected" (blue) methods (Sect.3).
relative error in H and LE due to errors in windspeed measurements is simply the relative error in u.C3 Air-temperature errorPropagating the relative error in H due to noise in θ gives -speed error, the error due to measurement errors in temperature is simply the relative error in θ.The error due to noise around temperature indirectly affects LE via the error in measuring humidity; this is discussed later.C4 Sensor height errorRandom error in the fluxes due to random error in measuring the height can be derived as follows:

Table 1 .
Characteristics of the sensors from the automatic weather station, eddy-covariance systems and profile mast.Shown are random errors in measurements provided by the manufacturer and those used in this study.

Table 3 .
Evaluation of high-frequency losses in H ec and LE ec fluxes.Spectral adjustment: percentage of flux losses estimated by adjusting an inertial subrange on the high-frequency end of the mean cospectra of w with θ and of w with q, for the three windregime subsets.Multi-resolution decomposition (MRD) ratio: median ratio calculated over each subset of fluxes calculated at the shortest timescale of the MRD cospectra and the flux calculated at the gap scale of the MRD cospectra.