Random uncertainties of flux measurements by the eddy covariance technique

Large variability is inherent to turbulent flux observations. We review different methods used to estimate the flux random errors. Flux errors are calculated using measured turbulent and simulated artificial records. We recommend two flux errors with clear physical meaning: the flux error of the covariance, defining the error of the measured flux as 1 standard deviation of the random uncertainty of turbulent flux observed over an averaging period of typically 30 min to 1 h duration; and the error of the flux due to the instrumental noise. We suggest that the numerical approximation by Finkelstein and Sims (2001) is a robust and accurate method for calculation of the first error estimate. The method appeared insensitive to the integration period and the value 200 s sufficient to obtain the estimate without significant bias for variety of sites and wide range of observation conditions. The filtering method proposed by Salesky et al. (2012) is an alternative to the method by Finkelstein and Sims (2001) producing consistent, but somewhat lower, estimates. The method proposed by Wienhold et al. (1995) provides a good approximation to the total flux random uncertainty provided that independent cross-covariance values far from the maximum are used in estimation as suggested in this study. For the error due to instrumental noise the method by Lenschow et al. (2000) is useful in evaluation of the respective uncertainty. The method was found to be reliable for signal-to-noise ratio, defined by the ratio of the standard deviation of the signal to that of the noise in this study, less than three. Finally, the random uncertainty of the error estimates was determined to be in the order of 10 to 30 % for the total flux error, depending on the conditions and method of estimation.


Introduction
The eddy covariance (EC) method is the most direct and defensible way to measure vertical turbulent fluxes of momentum, energy and gases between the atmosphere and biosphere.Considering an optimal measurement setup and a standardised scheme for post-field processing of the measured EC raw data, we can assume that the systematic error is minimised, and then the random error of the fluxes is typically dominating the EC flux measurement uncertainty at short timescales.The accuracy of flux random error estimates becomes important for interpretation of measurements especially when detecting small fluxes in terms of turbulent exchange or signal-to-noise ratio (SNR) of the instrumentation.Moreover, it is desirable to estimate the total random uncertainty for each averaging period as well as to separate it into the main components, e.g.one-point sampling error and instrumental noise (Businger et al., 1986).For the uncertainty due to instrumental noise, the method proposed by Lenschow et al. (2000) has been recently applied to EC measurements not only for energy and CO 2 (Mauder et al., 2013;Mammarella et al., 2015) but also for CH 4 (Peltola et al., 2014) and N 2 O fluxes (Rannik et al., 2015).Few authors (Detto et al., 2011;Schmidt et al., 2012;Sturm et al., 2012;Peltola et al., 2013;Deventer et al., 2015) have used the method proposed by Billesbach (2011) as a means of estimating the random instrumental noise.This approach, also called "random shuffle method", consists of randomly shuffling one of the data records in time and then estimating the error as covariance between the two decorrelated time series.
Recently Langford et al. (2015) analysed in detail the uncertainties related to flux detection from the EC data with low SNR.The authors evaluated the impact of the time-lag Ü. Rannik et al.: Random uncertainties of flux measurements determination and called for caution since under low SNR condition the traditional methods of maximising the crosscovariance function can lead to a systematic bias in determined fluxes.The study also reviewed the approaches for estimation of flux random errors.For quantifying the flux uncertainty Langford et al. (2015) suggest using the method by Wienhold et al. (1995) and, following Spirig et al. (2005), suggest multiplying the flux error standard deviation by a factor of three to obtain the limit of detection at 99 % confidence level.The method by Lenschow et al. (2000) to calculate the effect of instrumental noise on the flux error was also validated for data with low SNR by Langford et al. (2015).They compared the method with estimates derived from the rootmean-square (RMS) deviation of covariance of white noise and vertical velocity records and found that the error was not sensitive to the type of distribution of the noise and the RMS approach was consistent with the method by Lenschow et al. (2000).
In the current study we review available methods for the random error estimation of turbulent fluxes, which are widely used by the flux community.We perform calculation and analysis of flux errors by considering different error formulations described in Sect. 2.
We use the measured natural turbulent records for (i) quantitative comparison of the error estimates by Finkelstein and Sims (2001), Salesky et al. (2012), Wienhold et al. (1995), Lenschow et al. (2000) and Billesbach (2011) and (ii) evaluation of sensitivity of error estimates on numerical approximations and calculation details.Based on the analysis we provide recommendations regarding the choice of the flux random error estimates, together with calculation guidelines for numerical evaluation.
In addition, we generate artificial records with predefined statistical properties characteristic to atmospheric turbulence to (iii) evaluate the flux error estimates with high accuracy.Numerically evaluated error estimates are compared with the analytical predictions to validate the theoretical expressions for different error estimates.From simulated time series the calculated error estimates allow us also to (iv) evaluate the uncertainty of the flux random errors.

Theory
Turbulent fluxes averaged over a limited time period have random errors because of the stochastic nature of turbulence (Lenschow et al., 1994;Rannik et al., 2006) as well as due to noise present in measured signals (Lenschow and Kristensen, 1985).

Random error of the flux
The random error of the flux defined by F = w s = (w − w )(s − s ) , where the angle brackets denote ensemble averaging, w the vertical wind speed and s the scalar, can be evaluated as the standard deviation of the covari-ance, hereafter in the paper denoted by δ F , being the measure of 1 standard deviation of the random uncertainty of turbulent flux observed over an averaging period T .Theoretically, there are several ways to approximate the same error estimate.
For stationary time series, in the limit T → ∞, the flux random error can be expressed by using the instantaneous flux ϕ = w s = (w − w)(s − s) statistics according to Wyngaard (1973) and Lenschow et al. (1994): where the overbar denotes time averaging.
The integral timescale (ITS) of ϕ, τ ϕ , is defined according to where is the autocovariance function of ϕ, t the time delay and σ 2 ϕ is the variance of ϕ.Equation ( 2) can be used directly to estimate the timescale τ ϕ by integration of the auto-covariance function of ϕ, calculated from the high-frequency data records.Rannik et al. (2009) estimated the timescale τ ϕ and compared with simple parameterisations used in practical applications.The timescale τ ϕ was converted to a corresponding normalised frequency n ϕ by using n ϕ = f ϕ (z−d) U and τ ϕ = 1 2πf ϕ , where z is the observation height, d the displacement height and U denotes the average wind speed.For unstable conditions the value 0.24 (for temperature fluxes) to 0.27 (aerosol particle number concentration fluxes) was obtained for n ϕ .It was established that normalised frequency was not a function of stability under unstable conditions.Under stable stratification the frequency n ϕ was determined to increase with stability, which was parameterised by The flux uncertainty estimate according to the spectral or the Fourier method is defined as (Lenschow and Kristensen, 1985;Rannik and Vesala, 1999) where the spectrum of time series x = w, s, S x , can be represented as the squared magnitude of the Fourier transform (FT) of x, S x = |FT(x)| 2 , with normalisation over frequencies f is assumed as The cross-spectrum can be represented as S ws = FT(w) FT(s) * , where * stands for complex conjugate.Perhaps the most frequently used method to estimate the flux error is the method equivalent to the spectral method proposed by Lenschow and Kristensen (1985) in the time domain as R s for scalar s) represents the auto-covariance and R ws (t ) = (w(t) − w)(s(t + t ) − s) the cross-covariance functions.Finkelstein and Sims (2001) The Eq. ( 5) above is numerically approximated by Finkelstein and Sims (2001) as

Method by
with suitably chosen value of m, where n equals the amount of data points within the averaging period and w w (p) = the authors showed that where the variance of the filtered flux is defined as σ 2 ϕ = ( ϕ − ϕ) 2 using the identity ϕ = ϕ, and C ϕ is a constant introduced for convenience as For practical application the authors suggested choosing a series of time filter widths t distributed evenly at logarithmic scale between t, min and t, max and to estimate the coefficient C ϕ by least squares according to Eq. ( 8), i.e. the re- t .The values t, min = 10 Z/U and t, max = T /10 were recommended (Z = z − d denotes the measurement height relative to displacement height d).
At the limit t = T Eq. ( 8) becomes equivalent to Eq. ( 1) and the flux error by the method of Salesky et al. (2012) is given by where we denote by Ĉϕ the least square estimate of C ϕ .For calculation of the error estimate δF, S we used the Matlab code available online by Salesky et al. (2012).
The value Ĉϕ enables to obtain also an estimate of the ITS τ ϕ from Eq. ( 9).The same applies also to the method by Finkelstein and Sims (2001).If calculating the flux variance σ 2 ϕ = (ϕ − ϕ) 2 in addition to the error estimate δF, FS , the ITS can be obtained from the definition of the error in Eq. ( 1).Thus both methods have the advantage of not requiring the calculation of the ITS for flux random error estimation, which can be uncertain, and can be used for inferring the ITS.
All the error estimates presented in Sect.2.1 are different methods for evaluation of the same flux error provided that the averaging period T is much larger than the ITSτ ϕ .

Flux random error due to instrumental noise
Random uncertainty of the observed covariance due to presence of noise in instrumental signal, assuming the white noise with variance independent of frequency, gives essentially the lowest limit of the flux that the system is able to measure.Such uncertainty estimate can be expressed in its simplest form as (e.g.Mauder et al., 2013) where σ 2 w and σ 2 c denote the variances of the vertical wind speed and scalar concentration and σ 2 n_w, f and σ 2 n_c, f represent the instrumental noise variance of the vertical wind speed and scalar concentration measurements as observed at frequency f , respectively.Following the atmospheric similarity relationships (under near-neutral conditions) σ w = 1.25u * and σ c = 1.25c * , where u * denotes the friction velocity and c * the flux concentration defining the flux via F = u * c * , the relative flux error due to instrumental noise can be expressed as where signal-to-noise ratios were defined by SNR x = σ x σ n_x, f .
The noise level of the sonic anemometers is typically a few hundredths of m s −1 (see e.g.Rannik et. al. (2015) and the results of this study, Sect.4.3.1 below).Typical values of σ w vary between 0.1 and 1 m s −1 .This yields for the sonic anemometers the SNR values from 1 (under very low turbulence conditions) up to a few tens.The relative flux error resulting from the sonic anemometer noise is thus estimated to be up to 3 % for 30 min averaging period.This is small enough to be negligible from the practical point of view.
In most practical cases the instrumental noise is of importance for atmospheric compounds with low concentrations and fluxes.Therefore, we further assume negligible contribution of anemometer noise and express the flux error due to gas analyser noise only by using σ 2 n, f instead of σ 2 n_c, f for simplicity.For the white type of noise, which is typical for measurement instrumentation, the signal noise at frequency f can be expressed through its value at frequency 1 Hz by assuming averaging of the signal over periods f −1 and 1 s, respectively.This enables us to rewrite Eq. ( 13) as δ F, N = σ w σ n, f= 1 Hz √ T . The method to estimate noise contribution in second-to fourth-order moments in atmospheric measurements was derived rigorously by Lenschow et al. (2000) and applied to EC fluxes by Mauder et al. (2013).Lenschow et al. (2000) derived the method to estimate the instrumental random noise variance σ 2 n, f from the auto-covariance function of the measured turbulent record close to zero-shift, enabling us to determine the respective error for each half-hour flux averaging period under field conditions.In this study, the auto-covariance is linearly extrapolated to lag zero using the auto-covariance values at lags 1. ..5 (at 10 Hz frequency sampling rate) and the difference between this extrapolation and the observed auto-covariance value at lag zero (i.e. the variance of the time series) is the variance related to instrumental noise.The lag interval from 0.1 to 0.5 s was chosen as a compromise between accuracy and precision of the variance estimate.This method relies on the property of the noise that it is not correlated with the true signal variation.In the following, the noise variance estimate obtained according to the Lenschow et al. (2000) approach is denoted by σ 2 n, f and the respective flux error according to Eq. ( 13) by δF, N .
2.3 Other flux random error estimates 2.3.1 Method by Wienhold et al. (1995) Wienhold et al. (1995) use a method to calculate "the error in the flux determination, the flux detection limit", calculating the standard deviation of the covariance function R ws be-tween the intervals −50 ≤ t N ≤ −40 and 40 ≤ t N ≤ 50 as δF, W = 1 2p which is expected to underestimate the total flux error presented by Eq. ( 5).However, it produces larger error values compared to the flux error due to instrumental noise as defined by Eq. ( 13) and compared to the error estimate by Billesbach (2011, Eq. (18) below).We also demonstrate below (Sect.4.1) that the subsequent values of R ws are not statistically independent and the numerical estimate Eq. ( 15) converges to the theoretical expression Eq. ( 16) only if independent values of R ws are selected for calculation of the standard deviation as the error estimate δF, W . Billesbach (2011) Billesbach (2011) proposed a method to calculate the flux error estimate, which according to the authors was "designed to only be sensitive to random instrument noise".The error is calculated according to

Method by
where j ∈ [1. ..n] but the values are in the random order.We used 20 repetitions when calculating the uncertainty estimates with this method in order to obtain robust estimates (Billesbach, 2011).Random shuffling of the time series s makes it uncorrelated in time and decorrelates s (assumed to consist of the sum of turbulent signal and instrumental noise) from w, resulting in two independent variables.The variance of the product of two independent variables ϕ = w s is equal to the product of the variances σ 2 ϕ = σ 2 w σ 2 s .Therefore, the error of the average of ϕ over period T sampled at frequency where the scalar time series variance is the sum of the variances of turbulent scalar concentration and that of the noise, i.e. σ 2 s = σ 2 c + σ 2 n, f .Thus, the error according to the method by Billesbach (2011) is given by Eq. ( 18), which becomes equivalent to Eq. ( 13) when replacing σ s with σ n, f .Therefore the method interprets the variance of turbulent variation as a part of the random noise and produces an error estimate that overestimates the flux error due to instrumental noise only.Also, since turbulence spectrum does not follow the property described by Eq. ( 14), the error estimate according to Eq. ( 17) becomes dependent on the choice of the frequency f and cannot therefore be considered as a robust method to estimate the flux error.Lenschow and Kristensen (1985) have shown that the auto-covariance function for the Poisson type of noise has the form R s lowing the formal presentation in Eq. ( 5) and considering that w and s are independent due to shuffling, making the term corresponding to the product of the cross-covariances vanish, the error estimate for the method by Billesbach (2011) becomes which after integration becomes equivalent to Eq. ( 18).

The random error of the ensemble average flux
If an average over fluxes F i (i = 1. ..N ) is calculated, each of these representing a flux value observed over averaging period T and being characterised by an error δ F, i , then the error of the average flux 3 Materials and methods

Sites and measurements
Measurements from three different and contrasting sites with different surface properties and observation heights of about 23 m (forest site in Hyytiälä, SMEAR II), 2.7 m (Siikaneva fen site) and 1.5 m (Lake Kuivajärvi) above surface were used to evaluate the flux error estimates for June 2012 (July 2012 at Kuivajärvi) for measured temperature, carbon dioxide (CO 2 ), water vapour (H 2 O) and methane (CH 4 , Siikaneva only) fluxes.

SMEAR II (forest site)
The first set of measurements was done at the SMEAR II station (Station for Measuring Forest Ecosystem-Atmosphere Relationships), Hyytiälä, Southern Finland (61  nik (1998), and a more detailed description of the station and the measurements is in Hari and Kulmala (2005).

Kuivajärvi (lake site)
The second dataset is taken from Lake Kuivajärvi (61 • 50 N, 24 • 17 E), located close to the Hyytiälä Forestry Field Station and SMEAR II Station.Lake Kuivajärvi is a small humic boreal lake extending about 2.6 km in northwest to southeast direction, and it is a few hundred metres wide (surface area is 0.63 km 2 ).The measurement platform, firmly anchored from all the four corners, is located approximately 1.8 and 0.8 km from the northern end and southern end, respectively.Turbulent fluxes of momentum, heat, CO 2 and H 2 O are measured by an EC system (located on the above mentioned platform), which includes an ultrasonic anemometer (Metek USA-1, GmbH, Elmshorn, Germany) to measure the three wind velocity components and sonic temperature and the enclosed path infrared gas analyser LI-7200 (LI-COR Inc., Lincoln, NE, USA) that measures CO 2 and H 2 O concentrations.The data are sampled at 10 Hz, and the gas inlet is at 1.7 m above the water surface close to the sonic anemometer.More details about the site and measurements can be found in Mammarella et al. (2015).

Siikaneva (fen site)
The third dataset was collected at Siikaneva fen site (61 • 49.961 N, 24 • 11.567 E).The EC data used in this study were measured with a 3-D sonic anemometer (Metek USA-1, GmbH, Elmshorn, Germany) and one closed-path analyser (LI-7000, LI-COR Inc., Lincoln, NE, USA) for CO 2 and H 2 O.The sonic anemometer and the gas inlet was situated at 2.75 m above peat surface and the air was drawn to the analyser through a 16.8 m long heated sampling line.CH 4 mole fraction was measured with a closed-path gas analyser (FMA, Los Gatos research, USA).Further details about the site and measurements can be found in Peltola et al. (2013).

Flux processing
Turbulent fluxes and other statistics reported in the study were calculated over 30 min averaging period by block averaging approach (i.e.no detrending was applied if not mentioned otherwise) using the EddyUH software (Mammarella et al., 2016).
Prior to flux calculation, raw data despiking, conversion of CO 2 and H 2 O from wet to dry mole fraction and twodimensional coordinate rotation of wind vector was performed (Kaimal and Finnigan, 1994).The fluxes were corrected for frequency response underestimation at low and high frequencies by using the co-spectral transfer functions calculated according to Rannik et al. (1999) and Mammarella et al. (2009) respectively, together with in situ parameterisation of the co-spectral model derived from ensemble mean of measured temperature co-spectra (Mammarella et al., 2009(Mammarella et al., , 2015;;Peltola et al., 2013).
The measured data (wind speed and concentration records) were quality screened for spikes (all 30 min periods with a single data point exceeding physically meaningful value excluded) and, according to Vickers and Mahrt (1997), by applying the following statistics and selection thresholds: data with concentration skewness outside (−2, 2), kurtosis outside (1, 8) or Haar mean and Haar variance exceeding 2 were rejected.In addition, flux data were rejected if the second coordinate rotation angle was outside the range (−15 • , 15 • ) and at Kuivajärvi data were rejected if the wind was not blowing along the lake (directions 345-135 • and 170-290 • ).Additional quality screening was performed for flux stationarity by using the threshold value 1 (Foken and Wichura, 1996).

Superimposing Gaussian noise to the measured records
SNR is defined in the current study as SNR = . Records with low noise level (sonic temperature at the three sites) were used to evaluate the performance of the Lenschow et al. (2000) and Billesbach (2011) methods (Sect. 4.3.1 and 4.3.2, respectively), superimposing the measured signal with Gaussian noise.Natural variability of records in combination with different noise levels σ n, f = 10 Hz (0.025, 0.152 and 0.30 K) led to a range of SNRs from about 0.3 to 20, enabling us to determine the range and threshold of SNR where the Lenschow et al. (2000) method can be used.In addition, the temperature signal was high-pass filtered with a simple first-order low-pass filter in order to simulate scalar measurements (CO 2 , H 2 O, CH 4 , etc.) with a closed-path analyser performing as the low-pass filter to measured signal.Lowpass filtering was executed using 0.1, 0.3 and 0.6 s time constants.However, since the results for different time constants did not differ qualitatively, we present only the results for the time constant 0.3 s (Sect.4.3.1).

Simulation of artificial records
We generated artificial records with pre-defined statistical properties characteristic to atmospheric turbulence.Gaussian probability density functions were assumed for vertical wind speed and concentration time series.The atmospheric surface layer (ASL) similarity relationships were assumed for the variances of the records and the timescales of the autocorrelated processes were defined via normalised frequencies (Appendix A).
The analysis was carried out as following.First, we calculated the flux errors according to analytical expressions Eqs.(A4), (A5) and ( 18) for the three flux errors δ F , δ F, W and δ F, B , respectively.
Second, we calculated from the repeated simulated artificial records time series (N = 10 000) the fluxes and evaluated the error estimates according to where denotes ensemble averaging over a large number of records with duration T = 30 min.In order to estimate δ F, W , we retained auto-correlation of the series w and c but assumed no cross-correlation, i.e. with α wc = 0.In order to obtain the estimate for δ F, B we assumed uncorrelated time series w and c by taking α w = 0, α c = 0 and α wc = 0 in Eqs.(A1) and (A2).Third, we evaluated the error estimates δF, FS , δF, S , δF, W and δF, B according to Eqs. ( 6), ( 10), ( 15) and ( 17), respectively.These error values allow us to evaluate also the uncertainty of the random flux errors via calculation of the variance of the error estimates: 4 Results

Evaluation of error estimates based on simulated time series
The random errors of the covariance time series as presented in Appendix A can be derived analytically as the total error (Eq.A4) and as the error for the Wienhold et al. (1995) method (Eq.A5).Note that the flux error estimates by Finkelstein and Sims (2001) and Salesky et al. (2012) represent theoretically the same error given in general case by Eq. ( 5) and specifically for simulated time processes by Eq. (A4).Assuming the ASL similarity forms, the relative flux errors defined by δ F |F | −1 vary with wind speed and stability (Fig. 1).
The error δ F, W ignores the covariance part of the error expression and is therefore slightly smaller.Nevertheless, the method by Wienhold et al. (1995) provides good approximation of the total flux error δ F .The error according to the Billesbach ( 2011) method (δ F, B ) is much lower (estimated Table 1.Flux error estimates δ F , δ F, W and δ F, B obtained from theory and simulations.Height to wind speed ratios ZU −1 = 10 s (Z denotes the height above the displacement height i.e.Z = z − d).Number of simulated realisations N = 10 000.τ w and τ c denote the timescales of simulated series w and c in Eq. (A1).The values in the parenthesis indicate the uncertainty of the error estimates (obtained from Eq. 22) relative to the average values in percent i.e. 100 %×σ δF δF −1 . Anal.
Eq. ( 18) Simul.c σ F Eq. ( 21) Eq. ( 17) a Assuming auto-correlated but independent variables w(t) and c(t), with α wc = 0. b The method was modified such that in Eq. ( 15) the covariances R ws (t )were calculated with 10 s intervals within the lag ranges from −300 to −100 and +100 to +300 s, making essentially the estimates independent.c Assuming independent random variables w(t) and c(t), with α w = 0, α c = 0 and α wc = 0.
-1 -0.5 0 0.5 1 0  according to Eq. 18) and the relative error does not show dependence on wind speed.
According to ASL similarity functions the relative flux error is largest at near-neutral stability and decreases both for unstable and stable conditions.In stable cases the errors are smaller because the turbulent spectrum is shifted towards higher frequencies, resulting in more efficient averaging (over the same period T ) and reduced relative random uncertainty.This is similar to the effect of wind speed where higher wind speed implies higher-frequency turbulence and lower relative random uncertainty.In unstable cases the normalised frequencies n w and n c are independent of stability and the stability dependence of the relative error is caused by the functions in Eq. (A6).
We further analysed in detail the flux errors for the conditions characterised by ZU −1 = 10 s and three stability cases ZL −1 = −1, 0, 1. Table 1 indicates that the simulated error estimates according to Eq. ( 21) are close to the analytical expectations.Also the numerical methods by Finkelstein and Sims (2001), Salesky et al. (2012) and Billesbach (2011) produce similar values.However, the error estimated according to the Wienhold et al. (1995) method, δF, W in Eq. ( 15), deviates from what we expected theoretically (Eqs.16 and A5).The underestimation is particularly evident under neutral conditions when the relative error is largest.We modified the method to essentially decorrelate the individual R ws values by using values after every 10 s time shift and obtained much better correspondence to the expectations according to theory and simulations.Therefore we believe the error estimate δF, W can underestimate the flux error not only because it omits the www.atmos-meas-tech.net/9/5163/2016/Atmos.Meas.Tech., 9, 5163-5181, 2016 covariance term but also because the subsequent values of R ws are not independent, leading to underestimation of the total variability.The variability of the error estimates (calculated according to Eq. 22) is around 10 to 20 % for δF, FS and is slightly larger for δF, S and δF, W .The uncertainty of the random error estimate is larger for unstable and neutral conditions and smaller for stable.

Comparison of flux errors
The flux uncertainty estimates increase approximately linearly with the flux magnitude at all sites (Fig. 2 and Table 2).As expected from theory, the error estimate δF, FS yielded largest values for the uncertainty.The same magnitude would be expected from the estimates δF, S .However, the Salesky et al. (2012) method gave approximately ten to twenty percentages smaller uncertainty estimates than the Finkelstein and Sims (2001) method (Table 2 and Fig. 3).The possible reasons for this will be discussed in Sect. 5.The error δF, W gave few tens of percentages smaller uncertainty estimates compared to δF, FS (CO 2 : forest 29 %, fen 35 %, lake: 45 %; H 2 O: forest 33 %, fen 23 %, lake 15 %; CH 4 : fen 18 %).This is related to the difference between Eqs. ( 5) and ( 16): the estimate δF, FS included a cross-covariance term, whereas δF, W estimated the flux uncertainty related only to the auto-covariance term.Also, the method by Wienhold et al. (1995) as defined by Eq. ( 15) likely underestimates the error as discussed in Sect.4.1.The error δF, B obtained according to Lenschow et al. (2000) gave systematically higher uncertainty estimates than δF, N , as predicted in Sect.2.3.As δF, N is expected to estimate the flux uncertainties due to instrumental noise, δF, B is clearly a different error estimate.
Based on the linear regression statistics presented in Table 2, a few findings can be emphasised.The intercept values are small (compared to the flux error magnitudes) and imply that flux uncertainties tend to vanish with no turbulent exchange.Generally the relative flux error is larger over the forest (slope 0.14 and 0.18 for CO 2 and LE) compared to fen site (respective slopes 0.08 and 0.09).Surprisingly, the relative flux error is largest for CO 2 over the lake.This could be due to advective conditions for CO 2 .During the calm conditions, in particular at nights, CO 2 is expected to drain downhill towards the lake and accumulate, causing the concentration to increase and induce variation which is not related to local exchange over the lake.Additional variance in concentration record would impact also the flux error estimate.
Table 2. Linear fitting parameters between absolute value of the flux and flux uncertainty estimates (uncertainty = |flux| × S + I ) obtained from measurements at three sites.Robust fitting method was used to minimise the effect of outliers.

CO 2
Latent heat flux CH 4 The error due to instrumental noise ( δF, N ) is weakly correlated with the flux value, as expected from theory.The method δF, B gives correlated estimates with fluxes and this is expected due to correlation between the concentration variance and the flux (turbulent exchange naturally gives rise to concentration fluctuations).
For qualitative comparison with the behaviour of our theoretical model based on ASL similarity theory we constructed a plot similar to the one presented in Sect.4.1 (Fig. 1); see Fig. 4. The figure illustrates that the observed behaviour holds: the relative flux errors increase with increasing (z − d)U −1 (i.e. with lower wind speeds) and the stability dependence looks similar.Our theoretical model predicted the highest relative errors for near-neutral conditions.This holds for CO 2 , but in the case of H 2 O the peak is shifted towards the stable stratification side.This apparent dissimilarity of scalars could be the result of different source-sink behaviour.
The method δF, B produces wind speed invariant relative error estimates and again similar behaviour with stability as presented in Fig. 1.It is noteworthy that in Fig. 4 the estimates δF, FS and δF, S compare well whereas the estimate δF, FS tended to produce higher values, as represented by regression statistics in Table 2 and Fig. 3.This implies that the median values compare well whereas the larger error estimates that affect the regression slopes more tend to deviate.
In most cases the flux random uncertainty is dominated by the stochastic nature of turbulence and the instrumental noise is a minor part of the total uncertainty.The relative contribution of instrumental noise to the total uncertainty can be assessed by comparing the flux error estimates δF, FS and δF, N .2012) method ( δF, S , Eq. ( 10) in this study) plotted against the error estimated according to the Finkelstein and Sims ( 2001) method ( δF, FS , Eq. ( 6) in this study) for forest site.
Figure 5 shows that for CO 2 flux measurements at the fen site the instrumental noise causes around 5. ..10 % of the flux random uncertainty indicating that instrumental noise level is low enough for flux measurements at the given site.For sites with very low fluxes the situation might be opposite and the instrumental noise becomes limiting in detection of surface exchange.The calculated σ 2 n, f values were grouped according to SNR and ITS τ ϕ .Note that the ITS for ϕ is smaller roughly by a factor of 3 compared to the timescale of concentration c (cf. Eqs. 3 and A8).The results are shown in Fig. 6.At all three sites, the method by Lenschow et al. (2000) overestimated the noise variance if SNR > 3, the ITS was small and temperature signal was not low-pass filtered (cf.Fig. 6a-c).In these cases the signal was high compared to noise (SNR > 3), large part of the signal was at high frequencies (implying small ITS) and the high frequencies were not attenuated.When the ITS was larger, the noise was estimated more accurately, especially at the lake and fen sites.In contrast, if SNR < 3, the accuracy of the noise estimation did not significantly depend on the ITS, and the relative error of the noise estimation was in general within 10 %.
If the temperature signal was low-pass filtered before superimposing the signal with Gaussian noise, the accuracy of the noise estimation was improved (cf.Fig. 6d-f).When the SNR was below 5, the noise variance was estimated successfully (relative error within ±30 %), regardless of the ITS.Increasing the ITS improved the results, especially at the fen site (cf.Fig. 6e).In addition, stronger high-frequency signal damping generally increased the accuracy of the noise estimate (not shown).This result suggests that the Lenschow et al. (2000) method is more suitable for signals which are high-pass filtered (measurements with closed-path gas analysers) than for less attenuated signals (measurements with open-path gas analysers and sonic anemometers).
On the whole, the accuracy of the Lenschow et al. (2000) method depends on how strong is the signal relative to noise at high frequencies, since the noise is estimated using small time shifts close to the auto-covariance peak.The signal-tonoise ratio at high frequencies decreases if (i) the total SNR   2011) method ( δF, B , Eq. ( 17) in this study) plotted against flux uncertainty caused by instrumental noise (a) and a combination of noise and signal variances (b).If the method would estimate purely the instrumental noise, the points in (a) would follow 1 : 1 line.Instead, as predicted by Eq. ( 18), the points follow 1 : 1 line in (b).Sonic temperature data from the forest site were used with different level of noise added as shown in the legend.
measurement height) and the high-frequency variation in the signal is dampened.However, for measurements close to the ground with an open-path analyser the method does not perform equally well under all conditions.The reliability of the Lenschow et al. (2000) method in estimating the instrumental noise from the signal is determined by the SNR, as illustrated in the Fig. 6.Therefore, a priori knowledge on the instrumental precision characteristics is needed when analysing the outcome of the method.We analysed in more detail the performance of the Lenschow et al. (2000) method in Fig. 7.In case of low SNR and high ITS the method estimates the true noise with relatively small bias (Fig. 7a, d).For the same SNR but low ITS value the variance is strongly overestimated (Fig. 7b, e).We argued earlier that low-pass filtered signals enable to obtain better noise variance estimates (Fig. 6).However, in case of high SNR (Fig. 7c, f) the method leads to significant underestimation of the true variance.Thus filtered signals (instruments with not perfect frequency response) are not always preferred in terms of the method's ability to determine the signal noise.
In Table 3 we report the estimated signal noise statistics for the instruments used in the current study by defining reliable values according to criterion SNR < 3. The fraction of reliable estimates is low for the instruments with high precision characteristics.For example, the method does not typically work for estimation of the precision of horizontal wind speed measurements of the sonic anemometers.However, the method performed to produce the noise characteristics for the vertical wind speed components at all sites.This is at least partly due to the fact that the vertical wind speed variance Estimated from time series with 2 h (7200 s) duration to avoid large uncertainties of the auto-covariance function at long lags.Normalisation is done with the average between 400 and 600 s.Lines show medians and the areas interquartile ranges around the medians.The time constant used in the high-pass filtered case is given in the plot.Sonic temperature data from the forest site were used.
is smaller than that of the horizontal component, resulting in lower SNR for the vertical component.The estimated signal noise characteristics are in good correspondence with instrument specifications (where available) except for the CH 4 analyser, which had much better precision value than reported by the manufacturer.
Table 3. Amount and magnitude of reliable noise estimates obtained by using Lenschow et al. (2000) method.The reported values are estimates of 1 SD at 10 Hz sampling rate (i.e.σn, f = 10 Hz ).Noise estimate was considered reliable if SNR < 3, where SNR was calculated using the modal value of σn, f = 10 Hz obtained from the measurements.Also the performance specifications reported by instrument manufacturers are given for comparison.a Peak-to-peak value.For the comparison with 1 SD precision characteristic should be divided by a factor of 3. b Converted from a value reported for 1 Hz data using Eq. ( 14)

Random error according to Billesbach (2011)
Billesbach ( 2011) introduced the so-called "random shuffle" method to estimate the instrumental noise from EC measurements.However, as argued in Sect.2.3 this method does not estimate the flux error due to the instrumental noise since it mixes turbulent variation with noise and thus the error corresponding to the sum of the variances of the turbulent signal and the noise is deduced by the method.This is exemplified by Fig. 8a and b: the error estimated with the "random shuffle" method ( δF, B ) is equal to the calculation using the combined variances of signal and noise (Fig. 8b) and not to the calculation using the variance of the noise only (Fig 8a).
Thus the method cannot be expected to be suitable for estimation of the flux error due to (instrumental) noise in the signal.

Integration time in the Finkelstein and Sims (2001) method
Under different wind speed and stability conditions the ITS of turbulence varies and therefore it becomes relevant what would be the appropriate integration time for the method by Finkelstein and Sims (2001).The integration time of the flux error δF, FS is studied by varying m in Eq. ( 6) up to 1500 s.Further, to see the influence of different high-pass filtering techniques used in EC flux calculation, we applied the method to time series of sonic temperature from the forest site with following detrending options: mean removal (no detrending applied), linear detrending and auto-regressive highpass filtering with time constant 200 s.In general, the flux error estimate increases with integration time up to about 300 s (Fig. 9).At short integration times the high-pass filtered time series converge faster to the limit at longer integration times.This indicates contribution of low-frequency part of the spectra to the error estimates.Possibly part of this low-frequency variation is contributed by non-stationarity of the series.At larger integration times than about 300 s the error estimates essentially do not change.We choose further to normalise the error estimates with average value over the interval from 400 to 600 s.Plots for different sites (varying the observation level and surface type) and wind speed and stability influences as reflected by ITS classes indicate that integration time 200 s serves as an optimal choice for all conditions (Fig. 10).This would guarantee less than 10 % systematic underestimation of the flux error even in case of 25 % largest ITS values (ITS > 75th percentile) for and lake sites.The figure also illustrates that the cross-covariance term in Eq. ( 5) contributes 10 to 30 % of the error estimate, suggesting that the method by Wienhold et al. (1995), which ignores this term, underestimates the error by the same fraction.Ccov term in Eq. ( 6) Acov term in Eq. ( 6) Figure 10.Normalised uncertainty estimate according to Finkelstein and Sims (2001), δF, FS , as a function of integration time used in Eq. ( 6).Contribution of the cross-covariance (blue area) and auto-covariance (green area) terms in the Eq. ( 6) are shown separately.The values were normalised with average δF, FS at integration times between 400 and 600 s and grouped in two ITS classes before plotting.No detrending of time series was used.

Discussion and conclusions
Commonly applied random error estimates of turbulent fluxes were tested and compared in this study.The methods proposed by Finkelstein and Sims (2001) and Salesky et al. (2012), the error estimates δF, FS and δF, S according to Eqs. ( 6) and ( 10) respectively, approximate the random flux error defined as 1 standard deviation of the random uncertainty of turbulent flux observed over an averaging period T .The method by Salesky et al. (2012) is a novel approach that introduces a filtering method in error evaluation.However, the method stems from the same theory and is therefore expected to result in the same error estimates as Finkelstein and Sims (2001).We performed evaluation of performance of both methods with stationary artificial records as well as field measurements.Under stationary conditions (represented by simulated records) the method by Salesky et al. (2012) slightly underestimated the expected theoretical values (Table 1).The median flux error statistics obtained from measured data compared well (Fig. 4), whereas the obtained regression slope statistics suggested that the method by Finkelstein and Sims (2001) produced systematically higher error values than the Salesky et al. (2012) method (Table 2).It is also important that the methods did not produce distinctly different estimates which could be observed as outliers in Fig. 3.We interpret this that both of the methods perform well but the filtering embedded in the method by Salesky et al. (2012) effectively reduces the flux uncertainty introduced by non-stationarity in natural records.The methods by Finkelstein and Sims (2001) and Salesky et al. (2012) do not require estimation of the ITS and have an advantage compared to the methods which do so be-cause of uncertainty of the estimation.It was demonstrated by Salesky et al. (2012) that, under non-stationary conditions, their method performed better than the other methods involving direct evaluation of the ITS, namely the methods by Lumley and Panofsky (1964) and Lenschow et al. (1994).The methods proposed by Lenschow and Kristensen (1985) and in a discrete form by Finkelstein and Sims (2001) require just integration of the products of the covariance functions over sufficiently long time exceeding the ITS.While avoiding direct calculation of the ITS, the impact of nonstationarity is indirectly embedded in the method by Finkelstein and Sims (2001).As discussed above, we believe that the impact of non-stationarity is reduced due to the filtering in Salensky et al. (2012) method.The theory of turbulent flux and flux error calculation relies on the assumption of the stationarity.Violation of this assumption introduces additional uncertainty in time-average statistics including the flux random errors.Therefore superior performance of any of these two methods discussed above is not evident, and we suggest that both methods are reliable alternatives for flux error evaluation.Wienhold et al. (1995) defined an error estimate ( δF, W , Eq. ( 15) in this study), calculating the standard deviation of cross-covariance function over the lag interval far from the maximum.They called the error estimate as the detection limit of the flux.It was shown in the current study that the error estimate δF, W is in a good correspondence with δF, FS even though it does not rigorously define the same flux error.The method δF, W underestimates the flux random uncertainty by a few tens of percent due to the fact that it ignores the covariance part of the estimate in Eq. ( 5).We also demonstrated in this study that the error estimate δF, W as formulated by et al. (1995) underestimates the true flux uncertainty due to the fact that the cross-covariance estimates at neighbouring lags are not independent.To overcome this deficiency we suggest calculating the flux error variance from the cross-correlation values over longer lag interval but separated in time.For example, in the numerical exercise we chose the cross-covariance values with 10 s intervals within the lag ranges from −300 to −100 and +100 to +300 s.The modified approach reproduced the flux error values close to theoretical expectations whereas the original method underestimated the theoretical value up to 26 % (from three studied cases).
An alternative to one-point statistical estimation of the flux random errors as described in this study (Sect.2.1) is the two-tower approach, where the flux random error is evaluated by using the difference of the fluxes measured at two EC towers (e.g.Hollinger et al., 2004).The method assumes statistically similar observation conditions with independent realisations of turbulence at the two towers.Since the conditions are difficult to realise because of spatial correlation in measurements (e.g.Rannik et al., 2006), we suggest that the one-point statistical approach provides rigorous but more convenient method to estimate the flux random errors.Nevertheless, the two-tower approach was shown to give close results to the method by Finkelstein and Sims (2001) when similar weather conditions at the two sites were included in the analysis (Post et al., 2015).
The error estimates with very clear physical meaning are the total error resulting from stochastic nature of turbulence due to limited sampling in time and/or space, the methods by Finkelstein and Sims (2001) and Salesky et al. (2012), and to a good approximation the method by Wienhold et al. (1995), as well as the random error due to instrumental noise only.To estimate the latter from the field measurements (not from laboratory experiment) Lenschow et al. (2000) suggested calculating the signal noise variance from the difference between the signal auto-covariance at zero lag and the extrapolated value of the auto-covariance function to zero lag.The noise variance enables to calculate the flux error according to Eq. ( 11), which gives essentially the flux uncertainty under hypothetical conditions of no turbulent exchange (and thus variability) of scalar concentration (or vertical wind speed if the error due to anemometers noise is considered).Billesbach (2011) proposed the flux error estimate based on the product of vertical wind speed and concentration fluctuations, randomly re-distributing one of the series (denoted by δF, B , Eq. ( 17) in this study).The method was called the "shuffling method" and the authors proposed that the method was designed to only be sensitive to random instrument noise.We point out in this study that the method effectively adds the variance of turbulent scalar variation to noise variance and therefore the method is not equivalent to (overestimating) the method proposed by Lenschow et al. (2000) and also not to the Finkelstein and Sims (2001) method by strongly underestimating the total flux uncertainty.
Different flux error estimates have been assigned the meaning of the flux detection limit.For example, Wienhold et al. (1995) called their method as "detection limit of the flux".Billesbach (2011) suggested that the method they introduced "was sensitive to random instrument noise".The method by Lenschow et al. (2000) estimates the flux value that the system is able to detect within an averaging period T under hypothetical conditions of no turbulent variation of concentration.This error estimate serves as the theoretical lowest detection limit of the EC system.However, under natural turbulent exchange conditions the flux random uncertainty is contributed in addition to signal noise also by the stochastic nature of turbulence and the total flux error is larger, also meaning that the detection limit is larger than compared to the error introduced by the instrumental noise.Respectively, Langford et al. (2015) have defined "limit of detection" as 3 × δF, W (the uncertainty according to Wienhold et al. (1995) method) to give the flux measurement precision within 99 % confidence interval.By default most of the publications refer to 1 standard deviation of flux random variability (which corresponds to 68 % confidence intervals assuming normal distribution) when talking about the flux precision or random errors.If different confidence level is aimed, as by Langford et al. (2015), this should be explicitly stated.
The flux detection limit has been used also in conjunction with other flux measurement techniques.For example, in the case of chamber measurements the flux detection limit has been used to denote the flux error arising from all possible error sources.The traditional way to perform chamber measurements is to determine the gas concentration at several time moments during the chamber operation.In such data collection the sources of uncertainty are the imprecision related to gas sampling (either manual or automatic) as well as instrumental uncertainty (e.g.Venterea et al., 2009), leading to a measurement precision which is called a detection limit of chamber-based flux measurement system.It has to be noted that the flux detection limit of the chamber systems depends on several factors such as the type of the chamber and respective sampling method, the precision of the instrument, chamber dimensions and operation time.Therefore the flux detection limit of the chamber-based systems (which accounts for all possible sources of uncertainty) is comparable to the total stochastic error of the EC fluxes.
We also studied the performance of the Lenschow et al. (2000) and Finkelstein and Sims (2001) flux error estimation methods over different ecosystems and observation conditions.The performance of the Lenschow et al. (2000) method is affected by the SNR and the ITS of turbulence.We established that the method provides reliable estimates for SNR < 3 (in the statistical sense, single biased values can occur).However, no criterion based on the ITS could be provided as the results deviated among sites.
Application of the EC method requires stationarity of time series within averaging period (e.g.Foken and Wichura, Ü. Rannik et al.: Random uncertainties of flux measurements 1996).results in higher random uncertainty of the flux value and therefore the stationarity requirement has to be fulfilled if each 30 min or 1 h average value is expected to be statistically significant.We tested sensitivity of the flux errors derived by the Finkelstein and Sims (2001) method on integration time and high-pass filtering of the fluxes performed by mean removal, linear detrending and auto-regressive filtering.It was observed that the flux error increased with integration time up to about 300 s, revealing the influence of the low-frequency (possibly non-stationary) signal variance on the flux estimates.The high-pass filtered time series were less affected.For consistency, the flux errors should be calculated based on the same time series (in terms of filtering) as used for flux calculation.Apart from the impact of the low-frequency contribution to flux errors (and fluxes), which we believe is related to non-stationarity of the conditions, we observed that, in order to obtain δF, FS with good accuracy, integration of Eq. ( 6) over 200 s is sufficient for wide range of sites as well as observation conditions.Finkelstein and Sims (2001) originally performed summation over 20 s and suggested that the results changed less than 1 to 2 % for summation over 10 to 40 s for the dataset they used.Our results suggested that longer summation period is needed for robust determination of the error in case of tower-based measurements over variety of surfaces and wide range of observation conditions.
The EC fluxes are uncertain due to stochastic nature of turbulence by about 10 to 20 % under typical observation conditions.By using simulation of time series with statistical properties similar to natural records we deduce that the flux error estimates in turn are uncertain by about 10 to 30 %.

Data availability
Data are freely available upon request from the authors.

Figure 2 .
Figure 2. Flux error estimates δF, FS , δF, S , δF, W , δF, B and δF, N vs. the flux magnitude.Individual outliers were left outside the panels in order to show the majority of the data better.

Figure 3 .
Figure 3. Error estimate calculated based on the Salesky et al. (2012) method ( δF, S , Eq. (10) in this study) plotted against the error estimated according to the Finkelstein and Sims (2001) method ( δF, FS , Eq. (6) in this study) for forest site.

Figure 4 .Figure 5 .
Figure 4. Dependence of relative errors on (z − d)U −1 and (z − d)L −1 at the SMEAR II forest site.Data were binned before plotting (symbols connected by lines: medians; areas: interquartile range).For (a) and (c) only periods with |(z − d)L −1 | < 0.1 were used, and for (b) and (d) periods when 2 m s −1 < U < 4 m s −1 were selected.

Figure 6 .
Figure 6.Relative difference between the estimated ( σ 2 n, f ) and the actual (σ 2 n, f ) signal noise variance as a function of SNR and ITS calculated as z−d 2πn m U , where n m was found according to Eq. (3).The data were bin averaged before plotting.The values for τ in the subplots show the time constant used in low-pass (auto-regressive) filtering of the original temperature time series (i.e.before superimposing the noise).The dashed lines highlight the zero line (no systematic error in the noise estimate), whereas the dotted lines show the ±30 % thresholds.The top row shows data with no low-pass filtering.

Figure 7 .Figure 8 .
Figure 7. Three example cases for application of theLenschow et al. (2000) method with presented spectra (upper panels) and autocorrelation functions (lower panels).Case 1 (a, d): no filtering, low SNR, high ITS.Case 2 (b, e): no filtering, low SNR, low ITS.Case 3 (c, e): filtering used, high SNR, moderate ITS.Red lines correspond to the original signal; black lines correspond the original signal (highpassed filtered in case 3) superimposed with the noise; blue lines are the linear fits to the auto-correlation functions for lags from 0.1 to 0.5 s; green line in spectral plots indicates the +1 slope on log-log representation.

Figure 9 .
Figure9.Normalised errors δF, FS up to integration time 1500 s.Estimated from time series with 2 h (7200 s) duration to avoid large uncertainties of the auto-covariance function at long lags.Normalisation is done with the average between 400 and 600 s.Lines show medians and the areas interquartile ranges around the medians.The time constant used in the high-pass filtered case is given in the plot.Sonic temperature data from the forest site were used.
• 51 N, 24 • 17 E; 181 m a.s.l.).The station is surrounded by extended areas of coniferous forests and the tower of the EC measurements is located in a 50-year-old (in 2012) Scots pine (Pinus sylvestris L.) forest with dominant tree height of 17 m.The EC measurements were performed at 23 m height, approximately 6 m above the forest canopy.The wind speed components and sonic temperature were measured by an ultrasonic anemometer (Solent Research 1012R2, Gill Ltd, UK), and fast response CO 2 and H 2 O mole fraction by an infrared gas analyser (LI-6262, LI-COR Inc., Lincoln, NE, USA).A description of the measurement in micrometeorological context at SMEAR II station can be found in Ran- They are RMS values of 10 Hz data if not otherwise noted.