Articles | Volume 14, issue 7
Atmos. Meas. Tech., 14, 5089–5106, 2021
https://doi.org/10.5194/amt-14-5089-2021
Atmos. Meas. Tech., 14, 5089–5106, 2021
https://doi.org/10.5194/amt-14-5089-2021

Research article 28 Jul 2021

Research article | 28 Jul 2021

The high-frequency response correction of eddy covariance fluxes – Part 2: An experimental approach for analysing noisy measurements of small fluxes

The high-frequency response correction of eddy covariance fluxes – Part 2: An experimental approach for analysing noisy measurements of small fluxes
Toprak Aslan1, Olli Peltola2, Andreas Ibrom3, Eiko Nemitz4, Üllar Rannik1, and Ivan Mammarella1 Toprak Aslan et al.
  • 1Institute for Atmospheric and Earth System Research (INAR)/Physics, Faculty of Science, University of Helsinki, P.O. Box 68, 00014 Helsinki, Finland
  • 2Climate Research Programme, Finnish Meteorological Institute, P.O. Box 503, 00101 Helsinki, Finland
  • 3Dept. of Environmental Engineering, Technical University of Denmark (DTU), Lyngby, Denmark
  • 4Edinburgh Research Station, UK Centre for Ecology and Hydrology (UKCEH), Bush Estate, Penicuik EH26 0QB, UK

Correspondence: Toprak Aslan (toprak.aslan@helsinki.fi)

Abstract

Fluxes measured with the eddy covariance (EC) technique are subject to flux losses at high frequencies (low-pass filtering). If not properly corrected for, these result in systematically biased ecosystem–atmosphere gas exchange estimates. This loss is corrected using the system's transfer function which can be estimated with either theoretical or experimental approaches. In the experimental approach, commonly used for closed-path EC systems, the low-pass filter transfer function (H) can be derived from the comparison of either (i) the measured power spectra of sonic temperature and the target gas mixing ratio or (ii) the cospectra of both entities with vertical wind speed. In this study, we compare the power spectral approach (PSA) and cospectral approach (CSA) in the calculation of H for a range of attenuation levels and signal-to-noise ratios (SNRs). For a systematic analysis, we artificially generate a representative dataset from sonic temperature (T) by attenuating it with a first order filter and contaminating it with white noise, resulting in various combinations of time constants and SNRs. For PSA, we use two methods to account for the noise in the spectra: the first is the one introduced by Ibrom et al. (2007a) (PSAI07), in which the noise and H are fitted in different frequency ranges, and the noise is removed before estimating H. The second is a novel approach that uses the full power spectrum to fit both H and noise simultaneously (PSAA21). For CSA, we use a method utilizing the square root of the H with shifted vertical wind velocity time series via cross-covariance maximization (CSAH,sync). PSAI07 tends to overestimate the time constant when low-pass filtering is low, whilst the new PSAA21 and CSAH,sync successfully estimate the expected time constant regardless of the degree of attenuation and SNR. We further examine the effect of the time constant obtained with the different implementations of PSA and CSA on cumulative fluxes using estimated time constants in frequency response correction. For our example time series, the fluxes corrected using time constants derived by PSAI07 show a bias between 0.1 % and 1.4 %. PSAA21 showed almost no bias, while CSAH,sync showed bias of ±0.4 %. The accuracies of both PSA and CSA methods were not significantly affected by SNR level, instilling confidence in EC flux measurements and data processing in set-ups with low SNR. Overall we show that, when using power spectra for the empirical estimation of parameters of H for closed-path EC systems the new PSAA21 outperforms PSAI07, while when using cospectra the CSAH,sync approach provides accurate results. These findings are independent of the SNR value and attenuation level.

1 Introduction

Vertical turbulent fluxes of momentum, energy, and gases between the atmosphere and the biosphere measured by the eddy covariance (EC) technique are subject to both low- and high-frequency losses (Foken2008; Aubinet et al.2012). The physical limitations in instrument response times, spatial separation of instruments, line averaging, and air transport through the sampling tubes cause high-frequency losses (Aubinet et al.1999).

The EC sampling system acts as a low-pass filter on the flux, and the signal loss must be compensated for with the frequency response correction (FRC) during post-processing. The first step in the FRC is the description of the effect of the low-pass filtering of the measurement system, and for this, the transfer function approach has been widely used since it was first proposed by Moore (1986). The joint transfer function (H) that describes the low-pass filtering of the whole EC system can be determined theoretically or experimentally (Foken2008; Aubinet et al.2012). The theoretical approach involves various specific transfer functions that are estimated to represent different causes of flux loss. Conversely, in the experimental approach H is estimated from in situ measurements. Due to its simplicity, many studies have implemented the theoretical approach, which typically works well with fluxes measured by open-path EC systems, as well as momentum fluxes and sensible heat fluxes measured by sonic anemometers (Aubinet et al.2012). However, for complex EC systems the necessary information to calculate H is not available and needs to be estimated empirically. In addition, the time response of the system can vary with relative humidity (Ibrom et al.2007a), tube ageing (Mammarella et al.2009), and variations in the flow regime in the tube. Thus, the theoretical approach is not preferred for gas fluxes measured with closed-path EC systems, for which the experimental approach is therefore recommended (Aubinet et al.2012; Sabbatini et al.2018; Nemitz et al.2018).

In the experimental approach for closed-path systems, H is usually estimated from either the measured power spectra (i.e. PSA) or cospectra (i.e. CSA) of sonic temperature and the mixing ratio of the target gas (χ). Different studies use either PSA (Ibrom et al.2007a; Nordbo et al.2011; Fratini et al.2012; Sabbatini et al.2018) or CSA (Aubinet et al.1999; Humphreys et al.2005; Mammarella et al.2009; Peltola et al.2013). Also, some software packages used for EC flux calculation are based on PSA (e.g. EddyPro, see LI-COR Biosciences2020), while others are based on CSA (e.g. EddyUH, see Mammarella et al.2016).

Interestingly, there has not been much debate to date whether to use power spectra or cospectra to determine the time constant of the H (or response time), which characterizes the EC system's high-frequency response. Only recently, Wintjen et al. (2020) investigated the optimal method for high-frequency response correction, for fluxes of nitrogen compounds, recommending CSA. Ibrom et al. (2007a) argued for using PSA as the vertical wind speed (w) does not contain any relevant information for the spectral attenuation of the gas collection and data acquisition system, allowing us to describe sensor-related attenuation independently. This should in principle provide a better estimation of the time constant of the gas analysis system because the spectral data are not mingled with other components, such as, for example, sensor separation. In this approach, the effect of sensor separation is then treated explicitly with an additional correction step (Horst and Lenschow2009). In most cases the instrumental noise becomes visible in the high-frequency range of χ power spectra, and this has to be dealt with before the time constant of the gas sampling system can be estimated. The noise removal procedure is not well established, and this represents a major uncertainty in the PSA approach. This paper explores this uncertainty and proposes a novel, more robust approach to account for the noise in PSA.

On the other hand, this noise is often assumed not to correlate with the fluctuations in w, and therefore it does not contribute to the cospectra between w and χ. If this holds, this makes the CSA attractive for the estimation of the time constant because then the noise would effectively disappear from the measured signal. Yet, the use of the CSA relies on the correct determination of the time lag between w and χ, which may be difficult in the case of small fluxes due to noisier cross-correlation function, making the search for the absolute maximum harder (Langford et al.2015). Due to the above-mentioned reasons, empirically determined H can be a source of uncertainty for the FRC (Lee et al.2004). Additionally, the CSA approach inadvertently accounts for the phase shift caused by low-pass filtering, a topic discussed in our companion paper (Peltola et al.2021).

EC measurements conducted under low-flux conditions result in relatively high signal noise, i.e. low signal-to-noise ratio (SNR) (Smeets et al.2009). EC fluxes with low SNR are normally found in many ecosystems especially for methane (CH4) and nitrous oxide (N2O), as well as for other non-H2O and non-CO2 gas species and aerosol particles, in which any implementation of the FRC becomes uncertain (Rannik et al.2015; Nemitz et al.2018; Oosterwijk et al.2018). Low SNR can also be observed for CO2 in specific ecosystems (e.g. lakes) and seasons (e.g. senescence, winter dormancy) for CH4 over well-drained soils or peatlands during winter and for N2O in long periods outside the high emission periods (e.g. fertilizer applications, freeze–thaw cycles, or rain events), all of which, although small, significantly contribute to the long-term flux budgets and, hence, must be corrected to reduce the systematic bias. Thus, the investigation of uncertainties in commonly used FRC methods is of great importance for obtaining unbiased, harmonized, and continuous time series of gas fluxes measured by EC technique.

To our knowledge, the uncertainty in fluxes caused by the use of the PSA and the CSA has not been investigated systematically so far, motivating this study, which hypothesizes that the success of the PSA and CSA usage in FRC depends on the attenuation condition and the level of SNR. Consequently, we expect to see substantially different time constants, correction factors, and eventually different overall magnitudes of correction estimates with respect to the attenuation and SNR conditions. To test this hypothesis, we need a scalar dataset which represents different attenuation levels and noise conditions. Assuming spectral similarity between scalars, we apply different levels of attenuation and noise to sonic temperature time series (T) in order to generate a proxy representing an attenuated gas concentration dataset (e.g. CH4, N2O) with known characteristics. We use a first-order low-pass filter which solely depends on a single time constant (τLPF) to attenuate the signal. Then, we systematically contaminate the signal with white noise. We then assess which analysis approach most closely retrieves the true time constant used to degrade the flux in the first place. First, we try to retrieve the system time constants using the PSA and CSA and then compare those with original values (i.e. τLPF). Second, in order to demonstrate how variation in time constant estimation further affects the cumulative fluxes, we run a low-pass filter over 1-month-long time series of T and correct the attenuation with the FRC of Fratini et al. (2012) via implementing the time constants calculated in the first step. In Sect. 2, the theory of experimental FRC is summarized. In Sect. 3, materials used in this study and methods are explained. Results and discussion are then interpreted in Sect. 4.

2 Theory

2.1 Background of methods typically used to determine the system time response

In order to calculate the true unattenuated (i.e. frequency-response corrected) covariance (wχcorr) with the transfer function method, the measured covariance (wχmeas) is multiplied by a correction factor (Fcorr):

(1) w χ corr = w χ meas F corr .

One way to calculate Fcorr is to estimate the ratio of the integrated unbiased and biased cospectra as a function of frequency. In order to define the cospectrum of the unattenuated scalar under the assumption that the normalized cospectrum of all scalars has the same form (i.e scalar similarity), either the cospectra model (see Mammarella et al.2009) or the measured cospectra (see Fratini et al.2012) of T are used as a reference. Many studies used the surface layer models described by Kaimal et al. (1972) and based on Kansas experiments (see Moore1986). The attenuated cospectra are obtained by multiplication of the reference cospectra with the transfer function (H), which characterizes the filtering of the EC system.

Another way to calculate Fcorr is to simulate the attenuation with a recursive filter in time rather than in frequency space, in which Fcorr is defined as the ratio of the unattenuated and attenuated covariances (see Goulden et al.1997). In addition, based on the same approach, an experimental method was proposed by Ibrom et al. (2007a) to parameterize the correction factor separately for stable and unstable stratifications using meteorological data. For this method, a further modification was later suggested by Fratini et al. (2012) for the processing of large fluxes.

Regardless of the method chosen, H needs to be obtained either theoretically or empirically before Fcorr can be calculated. In the empirical approach, H can be determined using in situ measurements as a ratio of the normalized power spectra (for PSA) or cospectra (for CSA) of the attenuated scalar to those of an unattenuated scalar, e.g. T. In both approaches, in order to reduce the uncertainties in the low-frequency part of the spectra and to fulfil the assumption of spectral similarity, data must be selected from periods with rigorous stationary turbulent mixing. In addition, the power spectra and cospectra of T and χ should be normalized with their standard deviations so that they can be compared with each other (see their Eq. 2 in Ibrom et al.2007a).

For the PSA, H can be calculated using the power spectra of χ and T (Eq. 2), in which the effect of sensor separation should additionally be treated via the method proposed by Horst and Lenschow (2009). For PSA, H is derived as

(2) H PSA ( f ) = [ S χ ( f ) σ χ 2 ] [ S T ( f ) σ T 2 ] - 1 ,

where Sχ indicates the power spectrum of measured target gas mixing ratio, ST is the power spectrum of T, and variances (σχ2 and σT2) are calculated across the frequency range over which no attenuation occurs. Instrumental noise often becomes dominant in the high-frequency range of the power spectrum and also contributes to σχ (see blue line in Fig. 1).

Thus, prior to the calculation of Eq. (2), the noise contribution to the power spectra of Sχ should be removed (Ibrom et al.2007a). Finally, the frequency dependence of HPSA can be described through a sigmoidal curve which is characterized by the time constant (τ) of the measurement system (for more details see Peltola et al.2021, and references therein):

(3) H emp ( f ) = 1 1 + ( 2 π f τ ) 2 .

For PSA, Ibrom et al. (2007a) updated Eq. (2) by introducing a normalization factor which is used to secure the spectral similarity especially for small fluxes. As described in their Eq. (6), the time constant and the normalization factor (Fn) are obtained by fitting the following equation to the dampened and noise-free χ data:

(4) S χ ( f ) σ χ 2 = S T ( f ) σ T 2 F n 1 1 + ( 2 π f τ ) 2 .

In our study, we follow the same procedure (hereafter PSAI07) for the time constant calculation for the PSA, which is summarized in Sect. 3.2. In addition to PSAI07, we used a new comprehensive method for PSA, which is summarized in Sect. 2.2.

Alternatively, for the CSA, H is calculated as

(5) H CSA ( f ) = [ Co w χ ( f ) w χ ] [ Co w T ( f ) w T ] - 1 ,

where f is the natural frequency, Cowχ indicates the cospectrum of measured w and target gas mixing ratio χ, CowT is the cospectrum of measured kinematic heat flux, and wχ and wT are covariances calculated across a frequency range where the cospectra are not attenuated.

For CSA, τ is obtained via fitting Hemp to Eq. (5); however, there is an ongoing debate on whether the correct transfer function for the cospectra would be Hemp instead of Hemp (Moore1986; Eugster and Senn1995; Horst1997, 2000; Fratini et al.2012; Hunt et al.2016), which is related to the low-pass filtering time lag (i.e. phase shift) as discussed thoroughly by Peltola et al. (2021). They showed that the phase shift effect can be approximated well when square root is used. We therefore opt to apply Hemp to CSA.

2.2 Estimating the time constant from a noise-contaminated power spectrum

In the PSAI07 application to noisy data, the time constant is typically obtained via two separate fitting procedure steps following the approach of Ibrom et al. (2007a) as illustrated in Fig. 1 for a noisy spectrum. First, in order to remove the noise, the noise part of the attenuated noisy power spectra (solid blue line) is fitted with a line of an unconstrained slope (dashed blue line) within the marked-frequency domain in the high-frequency end, and then the line is extended towards lower frequencies and subtracted from attenuated noisy power spectra, yielding noise-free power spectra (black line) represented with the term Sχ(f)σχ2 in Eq. (4). The unattenuated and noise-free spectra of T (i.e. term ST(f)σT2 in Eq. 4) is represented with the red line. Second, the time constant is obtained via fitting Eq. (4) within the marked-frequency domain.

Regarding the noise removal procedure, this approach can be problematic as we will show in Sect. 4.1 below because if the resulting line (dashed blue line in Fig. 1) has a slope less than unity, its extrapolation to lower frequencies erroneously removes the true signal. In addition, the frequency domains used for fittings should be determined visually, requiring expertise in micrometeorology and signal processing, which limits its effective application in the research community.

Here we introduce a new alternative approach (hereafter PSAA21) that overcomes above-mentioned shortcomings as we will show in Sect. 4.1 below. It performs the estimation of the time constant and accounts for instrumental noise without removal in a single non-linear comprehensive fitting step. In this approach, the time constant is yielded via fitting Eq. (7) to directly attenuated noisy power spectra (solid blue line). The fitting is performed across the entire frequency domain (Fig. 1), indicating that the visual inspection is not needed. The brief derivation of Eq. (7) is shown below.

Equation (4) can be extended to the following equation, which includes also the noise component of the scalar:

(6) f S χ ( f ) σ χ 2 = f S T ( f ) σ T 2 F n 1 1 + ( 2 π f τ ) 2 + f S χ , n ( f ) σ χ 2 ,

where Sχ,n(f) is the power spectrum of the noise in χ. Here it is assumed that the noise and signal in measured χ time series are uncorrelated and hence two independent and additive components of the time series. In the case of white noise, Eq. (6) can be simplified to

(7) f S χ ( f ) σ χ 2 = f S T ( f ) σ T 2 F n 1 1 + ( 2 π f τ ) 2 + f b ,

where b is the ratio between noise variance and variance used to normalize the χ power spectrum, which is shown linearly in log–log scale. All terms in the equation have been multiplied by f because this is the standard normalization used to depict the spectral density functions (Fig. 1). The detailed derivation of Eq. (7) can be found in Appendix A.

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f01

Figure 1A diagram illustrating fitting procedures for PSA methods. Shown are the spectra of unattenuated and noise-free temperature (red line) and spectra of low-pass filtered and noisy scalar before (solid blue line) and after (black line) noise removal. For PSAI07, the noise is detected via fitting a line (dashed blue) to the high-frequency end of noisy scalar over the frequency range highlighted. Then, it is extended towards lower frequencies and subtracted from the noisy spectrum, yielding noise-free spectra. Later, the time constant is calculated via fitting Eq. (4) to noise-free spectra over the frequency range highlighted. For PSAA21, the time constant is obtained from one comprehensive fitting of Eq. (7) to noisy spectra over the whole frequency range highlighted.

Download

3 Materials and methods

3.1 Sites description and measurements

Two datasets measured with sonic anemometer in an EC set-up from the Siikaneva fen site were mainly used in this study. The site is located in southern Finland (6149.9610 N, 2411.5670 E; 160 m a.s.l.). The data were measured with 10 Hz sampling frequency using a 3-D sonic anemometer (model USA-1; METEK GmbH, Elmshorn, Germany). Further details about the site and measurements can be found in Peltola et al. (2013).

The first dataset (D1) was used for the time constant calculation (see Sect. 3.2). It contained 70 half-hourly EC data records measured in fully turbulent daytime conditions in the period from May to September 2013 with an average sensible heat flux of 114.3 W m−2, friction velocity of 0.3 m s−1, and wind speed of 2.1 m s−1.

The second dataset (D2) was used for the cumulative flux calculation based on T (see Sect. 3.3). This dataset was measured between 1 and 30 May 2013 and consisted of 1440 half-hourly periods for which fluxes were calculated.

Lastly, to demonstrate the performances of different methods in real-world data, we used CO2 dataset (D3) measured with 10 Hz sampling frequency using an infrared gas analyser (LI-7000, LI-COR, Lincoln, NE, USA). The measurements were simultaneously done with D1 at the same set-up at 2.75 m above the peat surface, where the centre of the sonic anemometer was displaced 25 cm vertically above the intake of gas analyser. The air was drawn to the analyser through a 16.8 m long heated inlet sampling line.

3.2 Data processing for time constant estimation

The data processing flow for all variants of PSA and CSA is summarized in Fig. 2. In order to generate the artificial dataset, which represents various known levels of SNR and attenuation, we first applied commonly used EC data processing procedures, i.e. de-spiking, 2-D coordinate rotation of the wind velocity vector, and linear de-trending to D1. We then degraded each half-hourly T time series with a first order low-pass filter in the spectral domain (see Sect. 3.2.1) and contaminated it with a prescribed amount of white noise in the time domain (see Sect. 3.2.2). For PSA, we calculated power spectra of T, following Sabbatini et al. (2018), normalized it by the total variance calculated within the frequency range of 0.0012 and 0.05 Hz, and averaged it into a logarithmically equally distanced frequency base. Then we took an ensemble average of the 70 power spectra. For PSAI07, we first removed the noise from the power spectra (see Sect. 2.2), then retrieved the time constant (i.e. τPSAI07) via fitting Eq. (4) within a frequency range, the lower limit of which is 0.01 Hz, while the optimal higher limit is defined via visual inspection. For PSAA21, we obtained the time constant (i.e. τPSAA21) via fitting Eq. (4) using the entire frequency domain.

In practice there is no time lag between T and w as both variables are measured with the same instrument. However, low-pass filtering introduces a time lag in the scalar of interest in addition to any physical time lag that may be caused by transport through sampling lines and/or sensor separation (Massman2000; Ibrom et al.2007b; Peltola et al.2021). In our case, since the scalar of interest is represented by the attenuated T, the time lag occurring between T and w was adjusted via the maximization of the cross-covariance for the time constant estimation in CSA. Next, we calculated cospectra of T and w, normalized them with the total covariance calculated within the frequency range of 0.0012 and 0.05 Hz, and averaged the data into exponentially spaced frequency bins. Then we calculated ensemble averages. Lastly, we derived the time constants for the original dataset (i.e. τCSAH,sync) via fitting the square root of Eq. (3). The fit was estimated for the frequency range from 0.01 to 2 Hz.

In summary, for three methods (two PSAs and one CSA), we assessed 45 different conditions each, combining five different attenuation levels with nine different SNRs. We repeated the same procedure 100 times to account for the uncertainty associated with the white noise generation and thus obtained 100 different values for the time constant for all attenuation and SNR levels for both PSA and CSA. The relevant results are shown in Sect. 4.2.

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f02

Figure 2Flow chart of the data processing for time constant calculation using the dataset D1: τLPF is the time constant of the first order filter, τPSAI07 and τPSAA21 represent the estimated time constants with power spectra (PSA) and τCSAH,sync with cospectra (CSA), and H is the spectral transfer function.

Download

3.2.1 Low-pass filtering

The dynamic performance of any EC measurement system can be approximated with a linear first-order non-homogeneous ordinary differential equation (Massman and Lee2002):

(8) τ LPF d χ O d t + χ O ( t ) = χ I ( t ) ,

where χO is the output of a scalar sensor, χI(t) is the true scalar concentration (i.e. input), and τLPF is the characteristic time constant of the sensor response (Massman and Lee2002). The spectral response (hLPF(ω)) of such a system can be obtained by the Fourier transform of the ratio of the output signal to the input signal, i.e. χO/χI (Horst1997):

(9) h LPF ( ω ) = 1 1 - j ω τ LPF ,

where ω=2πf and j=-1.

The desired (i.e. low-pass filtered) output data in the frequency domain can be estimated as

(10) Z O = h LPF ( ω ) Z I ,

where ZI is the Fourier transform of χI. From this, the low-pass filtered time series χO can be acquired by applying the inverse Fourier transform to ZO and taking the real part. In practice, the complex conjugate of hLPF is used to derive the correct temporal lag (scalar lag with respect to w).

We followed this procedure to filter T with five different τLPF values, i.e. 0.1, 0.2, 0.3, 0.4, and 0.5 s, corresponding to fc values of 1.60, 0.80, 0.53, 0.40, and 0.32 Hz, respectively.

3.2.2 Noise superimposition

In this study we use Gaussian white noise, which has equally distributed spectral densities across all frequencies (Stull2012), to contaminate the filtered signal in time-space. In order to generate time series with varying levels of SNR, we first generated white noise with unit standard deviation and multiplied the white noise with the standard deviation of the original T time series with different ratios (e.g. from 0.1 to 0.9). This represents the amount of noise compared to the amount of signal (e.g. from 10 % to 90 %). We then added this noise to the filtered signal. As a result, we obtained SNR values of 10.0, 5.0, 3.3, 2.5, 2.0, 1.6, 1.4, 1.2, and 1.1, which are equal to the ratio of the standard deviation of T and the standard deviation of white noise.

3.3 Data processing when estimating long-term budgets

Data processing steps when estimating the long-term budgets using the D2 dataset are summarized in Fig. 3.

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f03

Figure 3Flow chart of the data processing for the cumulative flux calculation using T of dataset D2. FPSAI07, FPSAA21, FCSAH,sync, and FREF are the cumulative fluxes corrected by following Fratini et al. (2012), in which Fcorr is calculated by implementing the time constants (for PSA: τPSAI07 and τPSAA21; for CSA: τCSAH,sync; for reference: τLPF) calculated in previous section.

Download

We applied regular EC data processing, which included de-spiking, coordinate rotation, and de-trending. Next, the T time series were deteriorated with values of τLPF that varied between 0.1 and 0.5 s to mimic a realistic range of scalar attenuations (mimicking the conditions of, for example, CH4 or N2O). Later, the low-pass-induced time lag was accounted for via maximization of the cross-covariance, which was followed by the calculation of the covariances. Then we applied the frequency response correction using Eq. (1), in which Fcorr was estimated using the method proposed by Fratini et al. (2012)1:

(11) F corr = f = f min f max CO ( f ) d f f = f min f max CO ( f ) H emp ( f ) d f ,

where CO equals the current T cospectrum when the absolute sensible heat fluxes exceeded 15 W m−2. For small fluxes we used a site-specific cospectral model (see Appendix C) for CO instead of parameterization of Fcorr proposed by Ibrom et al. (2007a). To make sure that the analysis was not affected by low data quality, we removed the fluxes with low friction velocity (u*<0.2 m s−1), unrealistic sensible heat fluxes and non-stationary conditions (Foken and Wichura1996), eliminating 554 half-hourly data points out of 1440 (ca. 38 %).

The time constant (τ) in Hemp is estimated by either PSA (i.e. τPSAI07 and τPSAA21) or CSA (i.e. τCSAH,sync), and this yielded cumulative fluxes FPSAI07, FPSAA21, and FCSAH,sync, respectively. For simplicity, in the FRC only median values of time constants of PSA and CSA ensembles estimated were used for each combination of SNR and attenuation. The reference fluxes (FREF) were also estimated with Fcorr using τLPF. We present the differences between cumulative fluxes as relative differences with respect to the reference (FREF) in per cent. In particular, as an example, the relative bias for PSAI07 is calculated as 100(FPSAI07-FREF)/FREF while for CSAH,sync as 100(FCSAH,sync-FREF)/FREF for several attenuation and SNR conditions. The relevant results are shown in Sect. 4.3.

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f04

Figure 4Effect of several low-pass filtering (τLPF=0.1, 0.3, and 0.5 s) and SNR (5, 3.3, and 2.5) on spectra of sonic temperature (T) of 70 half-hourly data records illustrated in logarithmic scale, where f is natural frequency, ST is spectral density, and σT2 is the variance of T. Shown are the spectra of raw measured sonic temperature (red lines) and of the artificially deteriorated (i.e. low-pass filtered and noise superimposed) sonic temperature before (dark blue) and after (black) noise removal following Ibrom et al. (2007a) through the subtraction of the linear fit (dashed blue lines) to the high-frequency end of the deteriorated spectra. The vertical lines mark the frequency range used for fitting for noise removal. The lower thresholds of the frequency range are 3, 2.3, and 2 Hz for the attenuation levels of 0.1, 0.3, and 0.5 s, respectively.

Download

4 Results and discussion

4.1 Time constant estimation with PSA

In order to illustrate the important steps in the data analysis in PSA (e.g. low-pass filtering, noise superimposition, and noise removal only for PSAI07), here we show how the shape of power spectra changes on a logarithmic scale for a few SNR values (5, 3.3, and 2.5) and low-pass filtering conditions (Fig. 4). The low-pass filter time constants (τLPF) of 0.1, 0.3, and 0.5 s result in fc values beyond which the signals become attenuated of ca. 1.6, 0.5, and 0.3 Hz. As the time constant increases, fc decreases, causing stronger attenuation. Referring to the upper panel of Fig. 4, for constant attenuation (τLPF=0.1 s), if the SNR decreases, noise becomes more visible and the line fit to the noise becomes more consistent. Similarly, referring to the left panel of Fig. 4, for constant SNR (e.g. for SNR=5), from top to bottom, as attenuation increases, the slope of the fit increases, and therefore the goodness of fit increases (Table 1). As discussed in Sect. 2.2, white noise causes a slope of 1, and thus any discrepancies from 1 result in overestimated noise removal. Based on Fig. 4 alone, it is evident that removing the noise from power spectra can be done with higher accuracy when the high-frequency attenuation increases and/or SNR decreases.

Table 1Results of the noise removal procedure applied in the power spectra approach following Ibrom et al. (2007a) (PSAI07): the values indicate the slopes of the fitted line to the high-frequency end of the spectrum in Fig. 4, with the coefficient of determination (R2) shown in the parenthesis as a function of τLPF and SNR. Note that for accurate noise removal the slopes should equal 1.

Download Print Version | Download XLSX

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f05

Figure 5Time constants calculated using the power spectrum approach, comparing (a) PSAI07 and (b) PSAA21 in several low-pass filtering conditions (τLPF= 0.1–0.5 s) as a function of signal-to-noise ratio (SNR) over the range 10–1.1, which corresponds to the amount of noise (e.g. 10 %–90 %). The solid black curve represents the median, while the shaded area represents interquartile ranges. The dashed line corresponds to the expected value (τLPF), which was used for the artificial low-pass filtering.

Download

The results of the time constant estimation are shown in Fig. 5a and 5b for PSAI07 and PSAA21, respectively. Results are presented as medians with 25th and 75th percentile ranges for the repeat simulations. It can be seen that for PSAI07 the interquartile range (IQR) expands as the amount of noise increases, while for PSAA21 it is small and almost constant with only a slight increase.

For PSAI07, the optimal frequency domains used for noise fitting were detected as 3, 2.6, and 2.3 to 5 Hz for the different attenuation levels of 0.1, 0.2, and 0.3 s, respectively. For attenuation levels of 0.4 and 0.5 s, we used the frequency range of 2 to 5 Hz.

PSAI07 overestimates the time constants with improving accuracy from low-attenuation to high-attenuation conditions regardless of the SNR values. The overestimation is likely due to the noise removal procedure, which further attenuates the high-frequency end of the spectra by removing the part of the signal together with noise (i.e. the noise is fitted with a slope <1). In this approach, the accuracy of the noise removal procedure can be improved by visual inspection and adjustment of the fitting range to provide slopes close to 1 and thus better fitting parameters. However, especially for low-attenuation conditions (e.g. 0.1, 0.2 s), the frequency range is not sufficient to detect the noise statistically, meaning that the linear fitting method is not ideal for differentiating in the spectral power the noise contribution from the real variations due to turbulence. In addition to the shortcoming of the linear fitting, the visual inspection of detecting the frequency ranges for both noise removal and the H fitting constitutes another uncertainty source due to its subjectivity. It requires expertise on the topic and is hard to automate when using software used for flux calculations (e.g. EddyPro). Moreover, in some cases, the optimal visual inspection might not be sufficient as the low-attenuation results in our study suggest. In our case we could not further improve the accuracy even though the exact attenuation and SNR level were known, which is not the case with real-world data.

PSAA21 successfully estimates the time constants regardless of the attenuation and SNR level. This is due to using the whole frequency range for fitting without separating the superimposed attenuated signal and noise. The most important advantage of the PSAA21 is that it does not require visual inspection.

We assumed that the noise contaminating the signal is white noise, which may not always be the case in real-world data. Thus, the accuracy of the PSAA21 depends on knowledge of the nature of the noise, which should be determined in advance and implemented in Eq. (6) by adjusting the last term that characterizes the noise. This said, in many EC studies, the type of noise is attributed to white noise (e.g. Launiainen et al.2005; Peltola et al.2014; Rannik et al.2015; Gerdel et al.2017; Wintjen et al.2020)2 but brown noise in other studies (Wintjen et al.2020). A simple approach to characterize the type of noise is either by examining the high-frequency end of spectra, which is similar to our study, or through the Allan variance, e.g. Werle et al. (1993). Alternatively, we conducted a brief investigation into the type of noise by comparing the power spectra of measurements with very low SNR and those of white and blue noise (see Appendix B). It should be noted, however, that there are situations when noise is complex and difficult to predict a priori.

4.2 Time constant estimation with CSA

Figure 6 illustrates the time-lag-corrected cospectra of three low-pass filtered cases (i.e. τ=0.1, 0.3, and 0.5 s) with the most noisy conditions (i.e. lowest SNR of 1.1) in addition to the unattenuated and noise-free cospectra (i.e. original). It shows that the white noise contamination did not cause linear increase in the high-frequency end of the cospectra, enabling the time-constant calculation without additional procedure related to noise removal.

CSAH,sync successfully estimates the time constants (Fig. 7). The variation around the expected values can be attributed to the shortcomings in the maximization of the cross-covariance used for the time-lag correction, the precision of which is limited to the sampling interval (e.g. 0.1 s). Hence, the importance of the precision of time-lag detection can be noted as a shortcoming of CSA. The use of square root might be the another reason for the variation as it is an approximation, not an exact representation of the phase shift effect (Peltola et al.2021).

Lastly, as expected, the level of noise does not affect the accuracy of the CSA method at any SNR level as the random noise does not correlate with w, which is a pivotal advantage of the CSA in general. It does, however, add random noise to the result that is greater than in the improved PSAA21 approach.

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f06

Figure 6Normalized ensemble cospectra for the original and various attenuated time series (i.e. τ=0.1, 0.3, and 0.5 s). White noise was added at a signal-to-noise ratio (SNR) of 1.1. The red line represents the unattenuated and noise-free cospectrum, while the dashed, dotted-dashed, and dotted lines show the noise-contaminated and low-pass filtered cospectra for the three attenuation levels.

Download

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f07

Figure 7Time constants calculated using the cospectral approach, i.e. CSAH,sync, for several low-pass filtering conditions (τLPF= 0.1–0.5 s) as a function of signal-to-noise ratio (SNR) over the range 10–1.1, which corresponds to the amount of noise (e.g. 10 %–90 %). The solid black line represents the medians, while the shaded areas represent interquartile ranges. The dashed line corresponds to the expected value (τLPF), which was used for the artificial low-pass filtering.

Download

4.3 Effect of time constant variations on cumulative fluxes

Figures 8 and 9 illustrate how variation in the estimation method of the time constant affected the cumulative fluxes in comparison to reference fluxes (see Sect. 3.3) calculated using dataset D2. Figure 10 shows the correction factors (Fcorr) for the case with τ=0.3 s and SNR of 2. The values vary between 1 and 1.3, indicating the frequency response correction as high as 30 %.

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f08

Figure 8Relative biases of the cumulative fluxes derived with different PSA methods, i.e. PSAI07 (blue) and PSAA21 (black), compared with the reference flux as a function of SNR for different attenuation timescales (0.1–0.5 s), calculated as, for example, 100(FPSAI07-FREF)/FREF.

Download

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f09

Figure 9Relative biases of the cumulative fluxes derived with the CSAH,sync compared with the reference flux as a function of SNR for different attenuation timescales (0.1–0.5 s), calculated as, for example, 100(FCSAH,sync-FREF)/FREF.

Download

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f10

Figure 10Correction factors (Fcorr) of half-hourly fluxes calculated with different approaches, i.e. PSAA21 (black cross), PSAI07 (blue point), CSAH,sync (green point), and the reference (red point) for the case with τ=0.3 s and SNR of 2.

Download

PSAI07-based fluxes showed a bias varying between 0.1 % and 1.5 %, whereas PSAA21 showed almost no bias, reflecting the more accurate time constant estimation previously shown. Fluxes based on CSAH,sync were close to the expected value with the bias of ±0.4 % with negligible response to the SNR level. These findings are consistent with the observed biases in time constant estimation, meaning that where the time constant and the low-pass filtering were overestimated (e.g. with the PSAI07 especially for τ=0.1 s and CSAH,sync for τ=0.1 and 0.5 s), the spectral correction factor and thus the fluxes were overestimated, too. In summary, the findings indicate that using the PSAA21 for the time constant estimation provides quite accurate FRC. On the other hand, for the CSA method, CSAH,sync reasonably approximates the time constants, hence FRC. These results are in agreement with the analysis presented in Peltola et al. (2021), showing that CSAH,sync approximates the effect of phase shift on the estimation of the time constant and flux correction factor. See Peltola et al. (2021) for more details.

4.4 Review of typical signal-to-noise ratios and response times encountered during closed-path flux measurements

The range of attenuation and SNR conditions reported in the literature is rather wide and varies depending on ecosystem type, the scalar of interest, data processing, a configuration of instruments, and set-up of the EC system. Ibrom et al. (2007a) examined fluxes of water vapour and CO2 measured over a temperate forest, both of which were disturbed by noise in the high-frequency range of power spectra (see their Fig. 2). They identified an fc of CO2 as low as 0.325 Hz (i.e. ca. τ=0.5 s), indicating strong attenuation, while H2O showed even stronger attenuation which increased with relative humidity to up to 0.010 Hz (i.e. ca. τ=16 s) due to the sorption/desorption effects on the sampling line internal walls. Langford et al. (2015) reviewed SNR values of various gases published in many studies. In their comprehensive analysis (their Fig. 4), the majority of the half-hourly datasets of isoprene and acetone, measured with a quadrupole-based proton transfer reaction-mass spectrometer (PTR-MS, Ionicon Analytik GmbH, Austria) above broad-leaf woodland, roughly showed SNRs of 0.5 and 0.3, respectively, and benzene measured (with the same instrument) in the urban environment showed an SNR of 0.3. Additionally, N2O measured with an older tunable diode laser (Aerodyne Research Inc., Billerica, MA, USA) over managed grassland showed an SNR of 1. Rannik et al. (2016) examined the random uncertainties in fluxes measured over forest, lake, and peatland ecosystems in the boreal region. According to their unpublished calculations, the SNRs of CO2 and H2O measured with an infrared gas analyser (LI-6262, LI-COR, Lincoln, NE, USA) were 2.6 and 10.3 for the forest site, respectively, while the SNRs of CO2, H2O, and CH4 measured with two closed-path analysers (i.e. LI-7000, LI-COR, Lincoln, NE, USA, for CO2, H2O, and FMA, and Los Gatos Research, Los Gatos, CA, USA for CH4) were 2.7, 24.3, and 5.4, respectively. At the same forest site, Kohonen et al. (2020) investigated fluxes of carbonyl sulfide (COS) with a newer generation Aerodyne quantum cascade laser spectrometer (QCLS) (Aerodyne Research Inc., Billerica, MS, USA) that also measures mole fractions of CO2 and H2O. They reported high attenuation for their measurement set-up (time constant of 0.68 s for their EC system) and high noise disturbance in their power spectra (their Fig. 7) with an SNR of about 1. During the measurement period, COS showed small fluxes near the detection limit, causing uncertainty in the calculation of H. Thus, they calculated the time constant using the cospectra of CO2, assuming that both fluxes are affected by the same attenuation.

In addition, N2O fluxes measured over urban areas (Järvi et al.2014) and CO2 fluxes over lake ecosystems (Mammarella et al.2015) can also be examples of low SNR conditions.

Given the wide range of SNR and attenuation conditions summarized above, we analysed only a limited range of SNR and attenuation. Also, the impact of the system response time depends on the position of the spectral peak frequency which changes not only with wind speed but also measurement height and surface roughness. Nevertheless, the analysis of the different approaches showed a systematic behaviour with respect to SNR level and attenuation conditions. This provides the opportunity to extend the results of this study beyond the examined values and guides the selection of the right method to find the relevant H.

The key constraint of the study was that we artificially simulated the various attenuation and SNR conditions. Thus, demonstrating the performance of the methods with real-world data is of great importance. Accordingly, here we provide an example via processing CO2 data D3 with the SNR of 2.3. The power spectra and cospectra of CO2 and T are shown in Fig. 11. The estimated time constants were 0.095, 0.075, and 0.173 s for PSAI07, PSAA21, and CSAH,sync, respectively. The details and statistics of fitting procedures are summarized in Table 2. The PSAA21 fit successfully approximates the noisy CO2 power spectra (R2=0.993), enabling us to assume that the noise is white. However, the slope of the linear noise fit of PSAI07 differs from 1 (0.82). We believe that this is due to the low attenuation so that the noise cannot become visible similarly as in the case of the artificial dataset with low attenuation shown in Fig. 4, resulting in overestimation as it erroneously removes part of the signal together with noise. The differences in time constants obtained via PSA and CSA approaches might be attributed to sensor separation, uncertainties in time-lag correction, and the use of square root in CSA.

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f11

Figure 11(a) Power spectra of (1) sonic temperature (red), (2) CO2 from the Siikaneva site measured simultaneously with D1 (black), and (3) noise-removed CO2 (dashed blue); shaded areas show the frequency range used for noise fitting for PSAI07, and time constant fitting for PSAI07 and PSAA21. (b) Cospectra of CO2 (black) and T (red) with vertical wind speed (w).

Download

Table 2Time constants obtained via different approaches (PSAI07, PSAA21, and CSAH,sync) for CO2 data from the Siikaneva site and relevant statistics.

Download Print Version | Download XLSX

5 Conclusions

Here we investigated the limitations of two commonly used approaches to empirically estimate the eddy-covariance (EC) transfer function needed for the frequency response correction of measured fluxes by analysing a temperature flux time series which was synthetically degenerated mimicking slow-response set-ups and noisy sensors. The first approach (PSA) is based on the ratio of measured power spectra, while the second (CSA) is based on the ratio of measured cospectra. For PSA, we examined two alternative approaches of accounting for the white noise contribution to the power spectra: i.e. PSAI07 and PSAA21. The latter is newly introduced here and does not require noise removal prior to the fitting of the response function. For CSA, we examined one approach utilizing the square root of the transfer function (H) with shifted w time series via maximization of the cross-covariance, i.e. CSAH,sync. We generated an artificial dataset using T with differing degrees of low-pass filtering (simulated damping) and additional random noise. The advantage of using artificial datasets is that the real values of frequency attenuation, SNR, and physical time lag are all known, allowing the precise comparison of the estimations and expected values. PSAI07 overestimated the noise contribution and consequently the signal loss and time constant for low-attenuation conditions, but better performance was found as attenuation increases. The new PSAA21 approach successfully estimated the time constants regardless of the attenuation and SNR level, identifying the noise and signal comprehensively and providing no bias and the lowest random uncertainty in the case of noisy data. However, the approach assumed that the signal is contaminated by white noise, but this is not necessarily always the case. Hence, prior to the calculation of the time constant with this method, the nature of noise must be known. CSAH,sync showed a slight deviation from expected values. That can be attributed to lack of precision in time-lag estimation and the use of square root.

We then examined the effect of the different approaches to estimate the time constants on the cumulative fluxes: fluxes corrected using the PSAI07-based time constants showed the bias varying between 0.1 % and 1.5 % in comparison with reference fluxes (see Sect. 3.3), for which PSAA21 showed almost no bias. By contrast, fluxes corrected using CSAH,sync showed the bias of ±0.4 %. This analysis, however, does not account for the shortcomings of the method, i.e. Fratini et al. (2012), and thus should be read as the sole effect of time constant estimation on the final fluxes. The comparison of commonly used methods are comprehensively done in the companion paper.

The SNR did not affect the accuracy of either PSA or CSA approaches, alleviating concerns on EC flux measurements with low SNR levels.

In summary, for the empirical estimation of parameters of H of closed-path EC systems, our findings showed that PSAA21 is the most accurate, precise, and robust method when power spectra are used. This finding was independent of SNR and degree of attenuation. On the other hand, using the square root of H in CSAH,sync, which is a good approximation of the attenuation of cospectra in the presence of phase shifts (Peltola et al.2021), provided the reasonably correct estimation of τ and Fcorr when cospectra are computed after with time-lag quantification by cross-covariance maximization.

Finally, given the constraints of this study, we encourage additional studies based on the real attenuation and SNR conditions, investigating also other types of noise contamination to provide a step forward in efforts to standardize the EC method, which is of great importance to avoid systematic biases of fluxes and improve comparability between different datasets.

Appendix A: Derivation of the PSAA21

In this section, we derive the new approach for the PSA (i.e. PSAA21) in detail. We assume that the measured χ contains two independent and additive components: the attenuated turbulent signal and the noise. We further assume that scalar similarity holds (i.e. power spectra of all scalars follow a similar shape). After these assumptions the power spectrum (Sχ) of χ can be written as

(A1) S χ ( f ) = σ χ 2 σ T 2 H S T ( f ) + S χ , n ( f ) ,

where σχ2 and σT2 are the variances of χ and T related to turbulent signal (no attenuation or noise), ST(f) is the power spectrum of T, H describes the attenuation of χ due to imperfect instrumentation, and Sχ,n is the noise in the χ measurements. Note that here the noise in the T measurements was neglected as being small except in very low turbulence conditions (when eddy-covariance measurements become problematic anyway and fluxes are removed by the u* filter). Assuming that the instrument measuring χ can be approximated by a first-order linear sensor and reordering terms, this yields

(A2) S χ ( f ) σ χ 2 = S T ( f ) σ T 2 F n 1 1 + ( 2 π f τ ) 2 + S χ , n ( f ) σ χ 2 ,

where the proportionality constant Fn was introduced due to the effect that noise adds to σχ2 and thus biases the normalization (see Ibrom et al.2007a). The spectral density is best shown in log–log space as a wide range of frequencies, and spectral densities are well displayed (Stull2012), with natural frequency (f) on the x axis and spectral density multiplied by f on the y axis. After multiplication by f, Eq. (A2) becomes

(A3) f S χ ( f ) σ χ 2 = f S T ( f ) σ T 2 F n 1 1 + ( 2 π f τ ) 2 + f S χ , n ( f ) σ χ 2 .

As described in Sect. 2.2, fSχ,n(f)σχ2 can be approximated by a model that is linear in frequency in logarithmic space as fSχ,n(f)σχ2=alog(f)+log(b). In the case of white noise, where a=1, this results in a simplified linear model as fSχ,n(f)σχ2=log(f)+log(b) and in non-logarithmic space as elog (f)+log (b), which is fb. Substitution into Eq. (A3) yields

(A4) f S χ ( f ) σ χ 2 = f S T ( f ) σ T 2 F n 1 1 + ( 2 π f τ ) 2 + f b .

It is worth mentioning that this method can be used to retrieve the variance of noise as well since b equals the ensemble-averaged noise power spectra (Sχ,n) divided by σχ2. However, we do not further examine it as it is not the main interest of this study.

Appendix B: Identification of instrumental noise

Here we perform a short analysis to characterize the type of noise in an example high-frequency CH4 time series by comparing the spectra of a turbulent scalar with those of artificially generated white and blue noise.

Methane fluxes from an upland forest site can be considered to be greatly affected by instrumental noise because fluxes are small and near the detection limit. Therefore, we used a forest methane mixing ratio dataset to identify the type of noise by comparing its spectra with the spectra of white and blue noise which was artificially generated and with the same standard deviation of the methane dataset.

The data were collected at the SMEAR II station (Station for Measuring Forest Ecosystem–Atmosphere Relations), Hyytiälä, southern Finland (6151 N, 2417 E; 181 m a.s.l.), which is a class 1 Integrated Carbon Observation System (ICOS) ecosystem station. The station is surrounded by extended areas of coniferous forests, and the EC tower is located in a 55-year-old (in 2017) Scots pine (Pinus sylvestris L.) forest with a dominant tree height of 19 m. The measurements were performed with 10 Hz sampling frequency at a height of 23 m, i.e. approximately 4 m above the forest canopy. A fast response laser absorption spectrometer (G2311-f, Picarro, Santa Clara, CA, USA) was used to measure the CH4 mixing ratio.

Power spectra of methane were calculated using measurements from 10 July at 12:00–14:30 UTC+2. The high-frequency range of the normalized frequency-weighted CH4 power spectrum follows the f+1 power law scaling, which is consistent with white noise and inconsistent with, for example, blue noise (Fig. B1).

https://amt.copernicus.org/articles/14/5089/2021/amt-14-5089-2021-f12

Figure B1Normalized spectra of methane concentration (red), white noise (black), and blue noise (blue) (fSs/σs2), all of which have the same standard deviation. The solid straight line is f+1 line.

Download

Appendix C: The site-specific cospectral model used in flux calculations

The cospectral model used in flux calculations for small fluxes (i.e. absolute value of sensible heat flux smaller than 15 W m−2) was calculated following Horst1997:

(C1) f Co ( f ) w χ = 2 π n / n m 1 + ( n / n m ) 2 ,

where n is the normalized frequency, and nm is the cospectral peak frequency, which is derived from in situ measurements in Siikaneva as

(C2) n m = 0.09 , z - d L 0 , 0.09 ( 1 + 4.5 ( z - d L ) 0.78 ) , z - d L > 0 ,

where z is the measurement height, d is the displacement height, and L is the Obukhov length, a measure of atmospheric stability.

Code and data availability

Data to reproduce Figs. 4, 5 6, 7, 8, 9, and 11 and code to demonstrate PSAA21 are available in Aslan et al. (2021, https://doi.org/10.5281/zenodo.4753827).

Author contributions

IM and OP designed the study, and TA processed and analysed the data. UR investigated the distortion caused by the quadrature spectrum in CSA. UR, TA, and AI investigated the limitations of the noise removal procedure of PSAI07. AI developed the PSAA21 following an idea that emerged during a stimulating discussion between all authors and investigated the square-root conundrum in CSA. TA wrote the manuscript with contributions from all co-authors.

Competing interests

The authors declare that they have no conflict of interest.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

This study was supported by ICOS and through the European Commission. Toprak Aslan is grateful to the Finnish National Agency For Education and the Vilho, Yrjö, and Kalle Väisälä Foundation for their kind support for funding. Olli Peltola is supported by the postdoctoral researcher project (decision 315424) funded by the Academy of Finland, and Eiko Nemitz acknowledges support by the Natural Environment Research Council award number NE/R016429/1 as part of the UK-SCAPE programme delivering National Capability.

Financial support

This research has been supported by the European Commission, H2020 Research Infrastructures (RINGO (grant no. 730944)), the Finnish National Board of Education (grant no. KM-18-10792), the Academy of Finland (decision no. 315424), the Natural Environment Research Council (award no. NE/R016429/1), and the Väisälän Rahasto.

Open-access funding was provided by the Helsinki University Library.

Review statement

This paper was edited by Glenn Wolfe and reviewed by George Burba and Marc Aubinet.

References

Aslan, T., Peltola, O., Ibrom, A., Nemitz, E., Rannik, Ü., and Mammarella, I.: The high frequency response correction of eddy covariance fluxes. Part 2: An experimental approach for analysing noisy measurements of small fluxes, Zenodo [data set], https://doi.org/10.5281/zenodo.4753827, 2021. a

Aubinet, M., Grelle, A., Ibrom, A., Rannik, Ü., Moncrieff, J., Foken, T., Kowalski, A. S., Martin, P. H., Berbigier, P., Bernhofer, Ch., Clement, R., Elbers, J., Granier, A., Grünwald, T., Morgenstern, K., Pilegaard, K., Rebmann, C., Snijders, W., Valentini, R., and Vesala, T.: Estimates of the annual net carbon and water exchange of forests: the EUROFLUX methodology, Adv. Ecol. Res., 30, 113–175, 1999. a, b

Aubinet, M., Vesala, T., and Papale, D.: Eddy covariance: a practical guide to measurement and data analysis, Springer Science & Business Media, Dordrecht, The Netherlands, ISBN 978-94-007-2350-4, 2012. a, b, c, d

Biosciences, L.-C.: EddyPro software instruction manual, LI-COR Inc., Lincoln, Nebraska, USA, 2020. a

Eugster, W. and Senn, W.: A cospectral correction model for measurement of turbulent NO2 flux, Bound.-Lay. Meteorol., 74, 321–340, 1995. a

Foken, T.: Micrometeorology, vol. 2, Springer, Berlin, Germany, https://doi.org/10.1007/978-3-642-25440-6, 2008. a, b

Foken, T. and Wichura, B.: Tools for quality assessment of surface-based flux measurements, Agr. Forest Meteorol., 78, 83–105, 1996. a

Fratini, G., Ibrom, A., Arriga, N., Burba, G., and Papale, D.: Relative humidity effects on water vapour fluxes measured with closed-path eddy-covariance systems with short sampling lines, Agr. Forest Meteorol., 165, 53–63, 2012. a, b, c, d, e, f, g, h

Gerdel, K., Spielmann, F. M., Hammerle, A., and Wohlfahrt, G.: Eddy covariance carbonyl sulfide flux measurements with a quantum cascade laser absorption spectrometer, Atmos. Meas. Tech., 10, 3525–3537, https://doi.org/10.5194/amt-10-3525-2017, 2017. a

Goulden, M. L., Daube, B. C., Fan, S.-M., Sutton, D. J., Bazzaz, A., Munger, J. W., and Wofsy, S. C.: Physiological responses of a black spruce forest to weather, J. Geophys. Res.-Atmos., 102, 28987–28996, 1997. a

Horst, T.: A simple formula for attenuation of eddy fluxes measured with first-order-response scalar sensors, Bound.-Lay. Meteorol., 82, 219–233, 1997. a, b, c

Horst, T.: On frequency response corrections for eddy covariance flux measurements, Bound.-Lay. Meteorol., 94, 517–520, 2000. a

Horst, T. and Lenschow, D.: Attenuation of scalar fluxes measured with spatially-displaced sensors, Bound.-Lay. Meteorol., 130, 275–300, 2009. a, b

Humphreys, E. R., Andrew Black, T., Morgenstern, K., Li, Z., and Nesic, Z.: Net ecosystem production of a Douglas-fir stand for 3 years following clearcut harvesting, Glob. Change Biol., 11, 450–464, 2005. a

Hunt, J. E., Laubach, J., Barthel, M., Fraser, A., and Phillips, R. L.: Carbon budgets for an irrigated intensively grazed dairy pasture and an unirrigated winter-grazed pasture, Biogeosciences, 13, 2927–2944, https://doi.org/10.5194/bg-13-2927-2016, 2016. a

Ibrom, A., Dellwik, E., Flyvbjerg, H., Jensen, N. O., and Pilegaard, K.: Strong low-pass filtering effects on water vapour flux measurements with closed-path eddy correlation systems, Agr. Forest Meteorol., 147, 140–156, 2007a. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Ibrom, A., Dellwik, E., Larsen, S. E., and Pilegaard, K.: On the use of the Webb–Pearman–Leuning theory for closed-path eddy correlation measurements, Tellus B, 59, 937–946, 2007b. a

Järvi, L., Nordbo, A., Rannik, Ü., Haapanala, S., Riikonen, A., Mammarella, I., Pihlatie, M., and Vesala, T.: Urban nitrous-oxide fluxes measured using the eddy-covariance technique in Helsinki, Finland, Boreal Environ. Res., 19, 108–121, 2014. a

Kaimal, J. C., Wyngaard, J., Izumi, Y., and Coté, O.: Spectral characteristics of surface-layer turbulence, Q. J. Roy. Meteor. Soc., 98, 563–589, 1972. a

Kohonen, K.-M., Kolari, P., Kooijmans, L. M. J., Chen, H., Seibt, U., Sun, W., and Mammarella, I.: Towards standardized processing of eddy covariance flux measurements of carbonyl sulfide, Atmos. Meas. Tech., 13, 3957–3975, https://doi.org/10.5194/amt-13-3957-2020, 2020. a

Langford, B., Acton, W., Ammann, C., Valach, A., and Nemitz, E.: Eddy-covariance data with low signal-to-noise ratio: time-lag determination, uncertainties and limit of detection, Atmos. Meas. Tech., 8, 4197–4213, https://doi.org/10.5194/amt-8-4197-2015, 2015. a, b

Launiainen, S., Rinne, J., Pumpanen, J., Kulmala, L., Kolari, P., Keronen, P., Siivola, E., Pohja, T., Hari, P., and Vesala, T.: Eddy covariance measurements of CO, Boreal Environ. Res., 10, 569–588, 2005. a

Lee, X., Massman, W., and Law, B.: Handbook of micrometeorology: a guide for surface flux measurement and analysis, vol. 29, Springer Science & Business Media, Dordrecht, The Netherlands, ISBN 1-4020-2264-6, 2004. a

Mammarella, I., Launiainen, S., Gronholm, T., Keronen, P., Pumpanen, J., Rannik, Ü., and Vesala, T.: Relative humidity effect on the high-frequency attenuation of water vapor flux measured by a closed-path eddy covariance system, J. Atmos. Ocean. Tech., 26, 1856–1866, 2009. a, b, c

Mammarella, I., Nordbo, A., Rannik, Ü., Haapanala, S., Levula, J., Laakso, H., Ojala, A., Peltola, O., Heiskanen, J., Pumpanen, J., and Vesala, T.: Carbon dioxide and energy fluxes over a small boreal lake in Southern Finland, J. Geophys. Res.-Biogeo., 120, 1296–1314, 2015. a

Mammarella, I., Peltola, O., Nordbo, A., Järvi, L., and Rannik, Ü.: Quantifying the uncertainty of eddy covariance fluxes due to the use of different software packages and combinations of processing steps in two contrasting ecosystems, Atmos. Meas. Tech., 9, 4915–4933, https://doi.org/10.5194/amt-9-4915-2016, 2016. a

Massman, W. and Lee, X.: Eddy covariance flux corrections and uncertainties in long-term studies of carbon and energy exchanges, Agr. Forest Meteorol., 113, 121–144, 2002. a, b

Massman, W. J.: A simple method for estimating frequency response corrections for eddy covariance systems, Agr. Forest Meteorol., 104, 185–198, 2000. a

Moore, C.: Frequency response corrections for eddy correlation systems, Bound.-Lay. Meteorol., 37, 17–35, 1986. a, b, c

Nemitz, E., Mammarella, I., Ibrom, A., Aurela, M., Burba, G. G., Dengel, S., Gielen, B., Grelle, A., Heinesch, B., Herbst, M., Hörtnagl, L., Klemedtsson, L., Lindroth, A., Lohila, A., McDermitt, D. K., Meier, P., Merbold, L., Nelson, D., Nicolini, G., Nilsson, M. B., Peltola, O., Rinne, J., and Zahniser, M.: Standardisation of eddy-covariance flux measurements of methane and nitrous oxide, Int. Agrophys., 32, 517–549, 2018. a, b

Nordbo, A., Launiainen, S., Mammarella, I., Leppäranta, M., Huotari, J., Ojala, A., and Vesala, T.: Long-term energy flux measurements and energy balance over a small boreal lake using eddy covariance technique, J. Geophys. Res., 116, D02119, https://doi.org/10.1029/2010JD014542, 2011.  a

Oosterwijk, A., Henzing, B., and Järvi, L.: On the application of spectral corrections to particle flux measurements, Environ. Sci.: Nano, 5, 2315–2324, 2018. a

Peltola, O., Mammarella, I., Haapanala, S., Burba, G., and Vesala, T.: Field intercomparison of four methane gas analyzers suitable for eddy covariance flux measurements, Biogeosciences, 10, 3749–3765, https://doi.org/10.5194/bg-10-3749-2013, 2013. a, b

Peltola, O., Hensen, A., Helfter, C., Belelli Marchesini, L., Bosveld, F. C., van den Bulk, W. C. M., Elbers, J. A., Haapanala, S., Holst, J., Laurila, T., Lindroth, A., Nemitz, E., Röckmann, T., Vermeulen, A. T., and Mammarella, I.: Evaluating the performance of commonly used gas analysers for methane eddy covariance flux measurements: the InGOS inter-comparison field experiment, Biogeosciences, 11, 3163–3186, https://doi.org/10.5194/bg-11-3163-2014, 2014. a

Peltola, O., Aslan, T., Ibrom, A., Nemitz, E., Rannik, Ü., and Mammarella, I.: The high-frequency response correction of eddy covariance fluxes – Part 1: An experimental approach and its interdependence with the time-lag estimation, Atmos. Meas. Tech., 14, 5071–5088, https://doi.org/10.5194/amt-14-5071-2021, 2021. a, b, c, d, e, f, g, h, i

Rannik, Ü., Haapanala, S., Shurpali, N. J., Mammarella, I., Lind, S., Hyvönen, N., Peltola, O., Zahniser, M., Martikainen, P. J., and Vesala, T.: Intercomparison of fast response commercial gas analysers for nitrous oxide flux measurements under field conditions, Biogeosciences, 12, 415–432, https://doi.org/10.5194/bg-12-415-2015, 2015. a, b

Rannik, Ü., Peltola, O., and Mammarella, I.: Random uncertainties of flux measurements by the eddy covariance technique, Atmos. Meas. Tech., 9, 5163–5181, https://doi.org/10.5194/amt-9-5163-2016, 2016. a

Sabbatini, S., Mammarella, I., Arriga, N., Fratini, G., Graf, A., Hörtnagl, L., Ibrom, A., Longdoz, B., Mauder, M., Merbold, L., Metzger, S., Montagnani, L., Pitacco, A., Rebmann, C., Sedlák, P., Šigut, L., Vitale, D., and Papale, D.: Eddy covariance raw data processing for CO2 and energy fluxes calculation at ICOS ecosystem stations, Int. Agrophys., 32, 495–515, 2018. a, b, c

Smeets, C. J. P. P., Holzinger, R., Vigano, I., Goldstein, A. H., and Röckmann, T.: Eddy covariance methane measurements at a Ponderosa pine plantation in California, Atmos. Chem. Phys., 9, 8365–8375, https://doi.org/10.5194/acp-9-8365-2009, 2009. a

Stull, R. B.: An introduction to boundary layer meteorology, vol. 13, Springer Science & Business Media, Dordrecht, The Netherlands, ISBN 978-94-009-3027-8, 2012. a, b

Werle, P., Mücke, R., and Slemr, F.: The limits of signal averaging in atmospheric trace-gas monitoring by tunable diode-laser absorption spectroscopy (TDLAS), Appl. Phys. B, 57, 131–139, 1993. a

Wintjen, P., Ammann, C., Schrader, F., and Brümmer, C.: Correcting high-frequency losses of reactive nitrogen flux measurements, Atmos. Meas. Tech., 13, 2923–2948, https://doi.org/10.5194/amt-13-2923-2020, 2020. a, b, c

1

Here we preferred using the square root to describe the true transfer function as it is a good approximation when maximization of the cross-covariance is used for the time-lag correction as shown by Peltola et al. (2021).

2

This also includes Ibrom et al. (2007a), who falsely interpreted the (white) noise as blue. This conclusion, however, was misguided by neglecting the fact that the power spectra were multiplied by f.

Short summary
Vertical turbulent fluxes of gases measured by the eddy covariance (EC) technique are subject to high-frequency losses. There are different methods used to describe this low-pass filtering effect and to correct the measured fluxes. In this study, we analysed the systematic uncertainty related to this correction for various attenuation and signal-to-noise ratios. A new and robust transfer function method is finally proposed.