Research article 28 Jul 2021
Research article  28 Jul 2021
The highfrequency response correction of eddy covariance fluxes – Part 1: An experimental approach and its interdependence with the timelag estimation
 ^{1}Climate Research Programme, Finnish Meteorological Institute, P.O. Box 503, 00101 Helsinki, Finland
 ^{2}Institute for Atmospheric and Earth System Research (INAR)/Physics, Faculty of Science, University of Helsinki, P.O. Box 68, 00014 Helsinki, Finland
 ^{3}Dept. of Environmental Engineering, Technical University of Denmark (DTU), Lyngby, Denmark
 ^{4}Edinburgh Research Station, UK Centre for Ecology and Hydrology (UKCEH), Bush Estate, Penicuik EH26 0QB, UK
 ^{1}Climate Research Programme, Finnish Meteorological Institute, P.O. Box 503, 00101 Helsinki, Finland
 ^{2}Institute for Atmospheric and Earth System Research (INAR)/Physics, Faculty of Science, University of Helsinki, P.O. Box 68, 00014 Helsinki, Finland
 ^{3}Dept. of Environmental Engineering, Technical University of Denmark (DTU), Lyngby, Denmark
 ^{4}Edinburgh Research Station, UK Centre for Ecology and Hydrology (UKCEH), Bush Estate, Penicuik EH26 0QB, UK
Correspondence: Olli Peltola (olli.peltola@fmi.fi)
Hide author detailsCorrespondence: Olli Peltola (olli.peltola@fmi.fi)
The eddy covariance (EC) technique has emerged as the prevailing method to observe the ecosystem–atmosphere exchange of gases, heat and momentum. EC measurements require rigorous data processing to derive the fluxes that can be used to analyse exchange processes at the ecosystem–atmosphere interface. Here we show that two common postprocessing steps (timelag estimation via crosscovariance maximisation and correction for limited frequency response of the EC measurement system) are interrelated, and this should be accounted for when processing EC gas flux data. These findings are applicable to EC systems employing closed or enclosedpath gas analysers which can be approximated to be linear firstorder sensors. These EC measurement systems act as lowpass filters on the time series of the scalar χ (e.g. CO_{2}, H_{2}O), and this induces a time lag (t_{lpf}) between vertical wind speed (w) and scalar χ time series which is additional to the travel time of the gas signal in the sampling line (tube, filters). Timelag estimation via crosscovariance maximisation inadvertently accounts also for t_{lpf} and hence overestimates the travel time in the sampling line. This results in a phase shift between the time series of w and χ, which distorts the measured cospectra between w and χ and hence has an effect on the correction for the dampening of the EC flux signal at high frequencies. This distortion can be described with a transfer function related to the phase shift (H_{p}) which is typically neglected when processing EC flux data. Based on analyses using EC data from two contrasting measurement sites, we show that the lowpassfilteringinduced time lag increases approximately linearly with the time constant of the lowpass filter, and hence the importance of H_{p} in describing the highfrequency flux loss increases as well. Incomplete description of these processes in EC data processing algorithms results in flux biases of up to 10 %, with the largest biases observed for short towers due to the prevalence of smallscale turbulence. Based on these findings, it is suggested that spectral correction methods implemented in EC data processing algorithms are revised to account for the influence of lowpassfilteringinduced time lag.
 Article
(2715 KB)  Companion paper
 BibTeX
 EndNote
The eddy covariance (EC) is the standard micrometeorological technique for measuring vertical turbulent fluxes of momentum, heat and gases in the atmospheric surface layer (Aubinet et al., 2012). Under certain conditions (flat terrain, homogeneous and stationarity turbulent flows, source/sink homogeneity, absence of chemical sources/sinks), the measured EC fluxes represent direct and continuous estimates of the surface exchange of energy and matter between the surface and the atmosphere and are then widely used to estimate energy and water balances, as well as the carbon budget of different type of ecosystems (Baldocchi, 2003; Vesala et al., 2008; Mammarella et al., 2015).
Gas flux measurements are made using a threedimensional sonic anemometer and a gas analyser, which are able to provide fastresponse measurements of turbulent fluctuations of vertical wind velocity and gas concentration (Aubinet et al., 2012). All EC systems, including openpath, enclosedpath or closedpath gas analysers, act as lowpass filters and thus cause a systematic bias to flux estimates. Flux loss at high frequency is due to the incapability of the measurement system to detect smallscale variation. The main reasons for cospectral attenuation are inadequate frequency response of the sensor, sensor separation and line averaging, as well as, in closedpath systems, the sampling of air through tubes and filters.
The frequency response correction is usually performed based on a priori knowledge of the system transfer function and the unattenuated cospectrum:
Here CF is the estimated spectral correction factor, Co^{w,χ} the normalised unattenuated cospectrum between scalar (χ) and vertical wind speed (w), f the frequency, H_{lpf} the total transfer function describing the lowpass filtering, and f_{1} and f_{2} the integration limits that are determined by the length of the averaging period and Nyquist frequency, respectively. The correction, performed by multiplying the measured covariance by the factor CF, always increases the absolute flux magnitude. Note that the highpass transfer function (Rannik and Vesala, 1999) is not included in Eq. (1) as here we focus only on the lowpass filtering effect of the EC system. The lowpass filtering transfer function H_{lpf} can be derived either theoretically or empirically (Foken et al., 2012a). The correction is in general different for momentum flux, sensible and latent heat fluxes, and the various gas fluxes and is specific to each EC system. In the theoretical approach, the atmospheric surface layer cospectral models (Moncrieff et al., 1997) are used, and H_{lpf} is calculated as a convolution of specific transfer functions representing different causes of flux loss, the equations for which can be found in Moncrieff et al. (1997) and Moore (1986). This approach works well for correcting the fluxes of momentum and sensible heat, as well as gas fluxes measured by openpath systems. Further, Horst (1997) and Massman (2000) have proposed alternative theoretical approaches, providing an analytical estimation of correction factors.
Alternatively, the empirical approach can be used, in which the model cospectra and H_{lpf} are estimated using in situ measurements. This is the preferable approach for closedpath systems, in which additional effects related to the inlet dust filters and rain caps (Aubinet et al., 2016; Metzger et al., 2016), potential variations in the sampling line flow rate, and absorption and/or desorption processes inside the sampling tube (Ibrom et al., 2007a; Nordbo et al., 2014; Runkle et al., 2012) may contribute to the EC system cutoff frequency, introducing transfer functions which are difficult to estimate a priori.
For this approach, different methods have been proposed for retrieving H_{lpf} from the measured power spectra or cospectra of the sonic temperature T and the target gas dry mole fraction (Fratini et al., 2012; Ibrom et al., 2007a; Mammarella et al., 2009; Nordbo et al., 2014). While the use of power spectra is the method recommended by the ICOS methodology (Sabbatini et al., 2018; Nemitz et al., 2018) and is implemented in the EddyPro software package, other studies and software packages (EddyUH) have supported the use of measured cospectra (Aubinet et al., 2000; Mammarella et al., 2009, 2016) for the empirical retrieval of H_{lpf}. The companion paper (Aslan et al., 2021) has investigated methodological issues and flux uncertainty related to these two methods under different attenuation conditions and signaltonoise ratio scenarios.
In EC systems, particularly those using closedpath gas analysers, the measurements of vertical wind velocity and gas concentration are not colocated, and a time delay between the two signals exists (Aubinet et al., 2012). The standard procedure is to determine this time delay for each averaging period (typically 30 min) by maximising the related crosscovariance function (in absolute terms) and taking the corresponding lag as the true signal delay. However, this time shift between the two signals depends not only on the system setup (e.g. air sample travel time from the inlet to the infrared gas analyser, IRGA, sampling cell, separation distance between the inlet and the centre of the sonic anemometer path), but it can additionally reflect the lowpass filtering of the measured signal (Massman, 2000; Ibrom et al., 2007a, b).
In this study, we investigate how the lowpassfilteringinduced phase shift affects the estimation of the highfrequency flux loss, and we show the implications that occur when H_{lpf} time constants are empirically derived from measured power spectra or cospectra and when related CFs values are applied to covariances estimated from the crosscovariance maximum. Finally, we present a method to take the phase shift effect into account when processing EC data. Towards these aims, we use EC data collected above a forest canopy (Hyytiälä) and a wetland (Siikaneva) in Finland covering a large range of attenuation conditions and integral turbulent timescale characteristics.
2.1 Transfer functions
The measured crossspectrum (Cr_{m}) between the fluctuations of vertical wind speed and the attenuated and lagged scalar can be described by (e.g. Massman, 2000)
where Z_{w} and Z_{χ} are Fourier transforms of the time series of vertical wind speed (w^{′}) and scalar (χ^{′}) fluctuations which are perfectly in phase and not attenuated, h_{χ} is the Fourier transform of the filter function that describes the response of the instrument measuring the scalar time series, $j=\sqrt{\mathrm{1}}$, and ϕ_{phys}=ωt_{phys}, where t_{phys} is a constant time shift (s) between w^{′} and χ^{′} and ω=2πf is the angular frequency, with f representing the frequency (Hz). The superscript ^{*} denotes complex conjugation. Equation (2) describes the crossspectra measured with a typical EC system employing a closedpath gas analyser in which the scalar fluctuations are attenuated (h_{χ}) and delayed in time (t_{phys}) with respect to w^{′} due to the gas sampling system (tubes, filters) and horizontal sensor separation. The signal travel time in the gas sampling line (t_{phys}) could be approximated from the volume of the sampling line and flow rate. Note that we assume that the sensor measuring w^{′} is perfect without any distortions in the measured w^{′} time series. This assumption does not limit the generality of the findings below. We also neglect any influence of sensor separation (Horst and Lenschow, 2009) or highpass filtering on flux attenuation.
Following, for example, Horst (1997) and Massman (2000), the scalar sensor is approximated to be a linear firstorder sensor for which h_{χ} can be written as
where τ is the time constant (s) of the filter, also called response time. Note that ${Z}_{\mathrm{w}}{Z}_{\mathit{\chi}}^{*}$ is the true crossspectrum (Cr) between unattenuated w^{′} and χ^{′} time series which are in phase, and it has a real (cospectrum Co) and an imaginary (quadrature spectrum Q) part (Cr = ${Z}_{\mathrm{w}}{Z}_{\mathit{\chi}}^{*}$ = Co+jQ). Recall also that the integral of Co equals the unattenuated covariance $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{\mathit{\chi}}^{\prime}}$, i.e. the vertical turbulent flux of scalar χ. Following previous studies (Horst, 1997; Massman, 2000), the quadrature spectrum (Q) can be assumed to be zero for stationary turbulent flow and hence Cr = ${Z}_{\mathrm{w}}{Z}_{\mathit{\chi}}^{*}$ ≈ Co. Now, using Euler's formula and some derivation (see Appendix A), we find
where the first term on the right equals the measured cospectrum (Co_{m}) and the second term the measured quadrature spectrum (Q_{m}). Note that whilst Q was assumed to be negligible, Q_{m} is in fact not equivalent to zero even when ϕ_{phys}=0, i.e. when the travel time of scalar χ^{′} signal in the sampling system is zero. This means that w^{′} and the attenuated ${\widehat{\mathit{\chi}}}^{\prime}$ (the ^{∧} denotes attenuation throughout the text), depending on frequency, can be out of phase due to lowpass filtering of Z_{χ} with h_{χ}. This has already previously been noted (e.g. Massman, 2000; Massman and Lee, 2002; Wintjen et al., 2020). Now, simply based on the real part of Cr_{m}, we can define a transfer function for the firstorder system (H) as
and a transfer function (H_{p}) for a generic phase shift ϕ as (Massman, 2000)
These two transfer functions together describe how the measured cospectrum (Co_{m}) deviates from the cospectrum calculated from ideal measurements free from any attenuation and time shifts (Co).
Similarly, forming the crossspectrum of the attenuated scalar with itself yields
where ${Z}_{\mathit{\chi}}{Z}_{\mathit{\chi}}^{*}$ is the unattenuated power spectrum. From this it follows that the same transfer function (i.e. H) applies to both power spectrum and cospectrum in the case that there is no phase shift between w^{′} and χ^{′} (i.e. H_{p} equals 1) and the quadrature spectrum Q is zero.
2.2 Timelag determination via crosscovariance maximisation
Ideally the time lag between w^{′} and ${\widehat{\mathit{\chi}}}^{\prime}$ would be estimated from the known dimensions of the gas sampling system and flow rate, in addition to other components causing time delay between signals. However, it is in practise difficult to keep track of all the components causing the time delay, and hence a practical (yet sometimes flawed; e.g. Langford et al., 2015) solution has been to estimate the time lag between w^{′} and ${\widehat{\mathit{\chi}}}^{\prime}$ via crosscovariance maximisation. This can be considered to be equivalent to finding a time shift t (and hence ϕ=ωt) that maximises the integral of Co_{m} (i.e. maximises the covariance between w^{′} and ${\widehat{\mathit{\chi}}}^{\prime}$).
The original aim of this crosscovariance maximisation is to account for the time lag between the time series w^{′} and ${\widehat{\mathit{\chi}}}^{\prime}$ that is induced by sampling lines (air filters, tubing) and horizontal sensor separation, in other words to find such a time lag that removes the ${e}^{j{\mathit{\varphi}}_{\text{phys}}}$ term in Eq. (2). After successfully accounting for the travel time (t_{phys}), the attenuation of the measured cospectrum Co_{m} (and hence flux) would be described only by H instead of H and H_{p} (see Eq. 4) since the timelag estimation sets H_{p}=1. However, the filter described by Eq. (3) results in a phase shift and hence produces an additional time lag between w^{′} and ${\widehat{\mathit{\chi}}}^{\prime}$, i.e. t_{lpf}. This can be seen from Eq. (4): Q_{m} differs from zero when ϕ=0 (see also Massman, 2000; Massman and Lee, 2002; Wintjen et al., 2020). Crosscovariance maximisation inadvertently derives the sum of t_{phys} and t_{lpf}, which implies that the time lag determined by crosscovariance maximisation is not the desired transport time lag. Using the sum of t_{phys} and t_{lpf} instead of t_{phys} induces a nonnegligible negative ϕ in Eq. (4). Hence when shifting the time series according to the crosscovariance maximisation, the transfer function related to the phase shift (H_{p}) can no longer be neglected and should be calculated using $\mathit{\varphi}=\phantom{\rule{0.125em}{0ex}}\mathit{\omega}{t}_{\text{lpf}}$, where t_{lpf} is the bias in the time lag when estimated via crosscovariance maximisation, in other words the lowpassfilteringinduced time lag (note the minus sign since H_{p} accounts for the overestimated time lag).
As the maximal possible cospectral energy content per unit frequency can be described with the amplitude spectrum (${A}_{\mathrm{m}}=\sqrt{{\text{Co}}_{\mathrm{m}}^{\mathrm{2}}+{Q}_{\mathrm{m}}^{\mathrm{2}}}$), it can be assumed that after the crosscovariance maximisation (time series shifted by t_{phys}+t_{lpf}), Co_{m} can be approximated by A_{m}. It is straightforward to derive A_{m} from Eq. (4):
Hence A_{m} is attenuated with $\sqrt{H}$ instead of H which describes the attenuation of Co_{m} in the case that there is no phase shift between w^{′} and ${\widehat{\mathit{\chi}}}^{\prime}$ as shown above. This suggests that after crosscovariance maximisation $\sqrt{H}$ approximates HH_{p} and the correct transfer function for cospectrum calculated after crosscovariance maximisation. However, $\sqrt{H}$ is not equivalent to HH_{p}. This relates to the fact that Q_{m} cannot be nullified with a constant time shift since lowpass filtering induces a frequencydependent time shift ($\frac{\mathrm{arctan}(\mathit{\omega}\mathit{\tau})}{\mathit{\omega}}$; see Massman, 2000). If Q_{m} could be nullified, then Co_{m}=A_{m} (see Eq. 10) and the attenuation of cospectra could be described accurately with $\sqrt{H}$. Note also that Eq. (10) applies universally, independent of ϕ.
The dependence between lowpassfilteringinduced time lag (t_{lpf}) and filter response time (τ) was estimated using the approximation $H{H}_{\mathrm{p}}\approx \sqrt{H}$ or ${H}_{\mathrm{p}}\approx \frac{\mathrm{1}}{\sqrt{H}}$. For this analysis the phase transfer function H_{p} was approximated by using the series expansion of the terms in H_{p} (Eq. 6). This approximation resulted in ${H}_{\mathrm{p}}\approx \mathrm{1}\frac{\mathrm{1}}{\mathrm{2}}(\mathit{\omega}{t}_{\text{lpf}}{)}^{\mathrm{2}}+\mathit{\omega}\mathit{\tau}(\mathit{\omega}{t}_{\text{lpf}})$. Similarly, $\frac{\mathrm{1}}{\sqrt{H}}$ can be approximated by $\mathrm{1}+\frac{\mathrm{1}}{\mathrm{2}}(\mathit{\omega}\mathit{\tau}{)}^{\mathrm{2}}$. Equating these two approximations up to the second order yields a quadratic equation in which the frequency dependence cancels out: ${t}_{\text{lpf}}^{\mathrm{2}}\mathrm{2}\mathit{\tau}{t}_{\text{lpf}}+{\mathit{\tau}}^{\mathrm{2}}=\mathrm{0}$. The solution of this simple equation gives t_{lpf}=τ. Thus, up to the secondorder approximation, which holds roughly up to the frequency value ω=1, the lowpassfilteringinduced time lag equals the time constant of the filter. Note, however, that at this frequency range (ω between 0 and 1) the filtering effect is very weak. The dependence between t_{lpf} and τ can be determined also by numerically solving ${H}_{\mathrm{p}}\approx \frac{\mathrm{1}}{\sqrt{H}}$ (which equals $\mathrm{cos}(\mathrm{2}\mathit{\pi}f{t}_{\text{lpf}})\mathit{\omega}\mathit{\tau}\mathrm{sin}(\mathrm{2}\mathit{\pi}f{t}_{\text{lpf}})\approx \sqrt{\mathrm{1}+{\left(\mathrm{2}\mathit{\pi}f\mathit{\tau}\right)}^{\mathrm{2}}}$) and using an assumption that t_{lpf} is proportional to τ, i.e. t_{lpf}=Cτ, where C is a proportionality constant. These analyses are continued in Sect. 4.1. Note that the dependency between t_{lpf} and τ derived here and in Sect. 4.1 was not used in the estimation of correction factors (see Sect. 3.2.1) but merely to analyse the factors controlling the dependency.
3.1 Measurement sites
Measurements from the Siikaneva fen and the Hyytiälä forest site (SMEAR II) were used. Both stations are part of the Integrated Carbon Observation System (ICOS) measurement station network.
The SMEAR II station is situated in southern Finland (61^{∘}51^{′} N, 24^{∘}17^{′} E; 181 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, above sea level). The station is surrounded by extended areas of coniferous forests, and the EC tower is located in a 57yearold (in 2019) Scots pine (Pinus sylvestris L.) forest with a dominant tree height of ca. 19 m. The EC measurements were performed with 10 Hz sampling frequency at 27 m height, whereas the zero plane displacement height (d) was 14 m. The wind speed components and sonic temperature were measured by a threedimensional ultrasonic anemometer (Solent Research HS1199, Gill Ltd, Lymington, UK), while carbon dioxide and water vapour mixing ratios were measured by an infrared gas analyser (LI7200, LICOR Biosciences, Lincoln, NE, USA). A 0.77 m long tube (4 mm inner diameter) was used to sample air for the gas analyser with a nominal flow rate of 10 L min^{−1}. The centre of the sonic anemometer was displaced 20 cm horizontally and 1 cm vertically from the intake of the gas analyser. The measurement setup was designed to closely follow the ICOS EC measurement protocols (Franz et al., 2018; Rebmann et al., 2018). The measurements were conducted between May and August 2019.
The Siikaneva fen site is located in southern Finland (61^{∘}49.9610^{′} N, 24^{∘}11.5670^{′} E; 160 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), consisting mainly of sedges (Eriophorum vaginatum, Carex rostrata, C. limosa) and Sphagnum species, namely S. balticum, S. majus and S. papillosum, with a height of ca. 10–30 cm. Further details about the site can be found in Riutta et al. (2007), Peltola et al. (2013) and Rinne et al. (2018). The EC measurements were conducted between May and August 2013. The data used in this study were measured with 10 Hz sampling frequency using a threedimensional sonic anemometer (Metek USA1, GmbH, Elmshorn, Germany) and a closedpath analyser (LI7000, LICOR Biosciences, Lincoln, NE, USA). The sonic anemometer and the gas inlet were situated 2.75 m above the peat surface, and the air was drawn to the analyser through a 16.8 m long heated inlet sampling line. The centre of the sonic anemometer was displaced 25 cm vertically above the intake of the gas analyser.
A short dataset (hereafter ${D}_{{\mathrm{S}}_{\mathrm{1}}}$) measured between 09:30 and 11:30 (UTC + 2) on 16 June 2013 at the Siikaneva site was used in the method validation (see Sect. 4.2), whilst longterm datasets from both Siikaneva (hereafter ${D}_{{\mathrm{S}}_{\mathrm{2}}}$) and Hyytiälä (hereafter D_{H}) were used to evaluate the accuracy of different spectral correction methods and their effect on CO_{2} and H_{2}O fluxes (see Sect. 4.3).
3.2 Data processing
Highfrequency EC data from the two sites were processed in order to evaluate (1) the accuracy of different spectral correction methods and (2) their effect on gas (CO_{2} and H_{2}O) flux estimates at these sites. The processing steps followed commonly accepted routines: (1) data were despiked using a method based on the running median filter (Brock, 1986; Starkenburg et al., 2016), (2) coordinates were rotated using doublerotation which aligned one of the horizontal wind components with the mean wind, setting the mean vertical and cross wind components to zero, (3) gas (CO_{2} and H_{2}O) mole fractions were converted pointbypoint to be relative to dry air if not already done internally by the gas analyser (LI7200), (4) turbulent fluctuations were extracted from the measurements using blockaveraging, and (5) time lags between gas or filtered T (see Sect. 3.2.2) and vertical wind time series were accounted for using crosscovariance maximisation. Power spectra and cospectra were calculated after these processing steps. Spectral corrections and the related correction factors (CFs) were calculated using four methods described below; see Sect. 3.2.1. After processing, the flux time series were quality filtered by removing periods with unrealistic sensible heat fluxes or friction velocities and highly nonstationary fluxes (Foken and Wichura, 1996). Sitespecific friction velocity thresholds (0.17 and 0.3 m s^{−1} for Siikaneva and Hyytiälä, respectively) were used to discard gas flux data during lowturbulence periods. Furthermore, at Hyytiälä wind directions between 150 and 230^{∘} were discarded since here the mast construction disturbed the airflow. At Siikaneva an omnidirectional sonic anemometer was used mounted on the top of a mast, and hence clearly disturbed wind directions were not identified. After quality filtering 4210 and 4722 30 min periods were available for analysis from Hyytiälä and Siikaneva, respectively.
Atmospheric stability was evaluated using the Obukhov length (L):
where u_{*} is the friction velocity, v_{k}=0.4 the von Karman constant, g the acceleration due to gravity, and θ_{v} virtual potential temperature. Overline denotes temporal averaging and primes deviations from the mean. From the Obukhov length, the stability parameter was calculated as $\mathit{\zeta}=\frac{zd}{L}$, where z is the measurement height and d is the zero plane displacement height.
3.2.1 Empirical estimation of the lowpass filtering of the signal
The total transfer function of an EC setup can be estimated empirically from measured w–$\widehat{\mathit{\chi}}$ and w–T cospectra:
or, alternatively, from the power spectra of the measured scalar and T:
where ${\text{Co}}^{w,\widehat{\mathit{\chi}}}$ denotes the cospectral density between w and $\widehat{\mathit{\chi}}$, and $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{\widehat{\mathit{\chi}}}^{\prime}}$ represents the covariance between w and $\widehat{\mathit{\chi}}$. ${S}_{\widehat{\mathit{\chi}},\widehat{\mathit{\chi}}}$ and ${\mathit{\sigma}}_{\widehat{\mathit{\chi}}}^{\mathrm{2}}$ denote the power spectral density and variance of $\widehat{\mathit{\chi}}$, respectively, and the corresponding terms are defined for T. When estimating H_{e,cs} and H_{e,ps}, the covariances and variances were calculated from unattenuated frequencies (between 5 × 10^{−3} and 5 × 10^{−2} Hz) (Foken et al., 2012a; Sabbatini et al., 2018). An additional normalisation factor (F_{n}) was also incorporated in order to take into account the fact that the variances and covariances in Eqs. (12) and (13) may be subject to various highfrequency losses (Ibrom et al., 2007a). Here, w and T were approximated to be free from any lowpass filtering effects, and thus Co^{w,T} and S_{T,T} could be used as a reference. H_{e,cs} and H_{e,ps} were calculated from ensembleaveraged spectra from measurement periods that fulfilled the following criteria: flux stationarity (Foken and Wichura, 1996) was below 0.3, u_{*} > 0.1 m s^{−1}, and $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{T}^{\prime}}$ > 0.02 K m s^{−1}. For CO_{2} it was also required that its turbulent flux was below −0.05 $\mathrm{\mu}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{mol}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, and for H_{2}O it was required that the flux was directed upwards. For H_{2}O the transfer functions were determined in relative humidity (RH) bins, and the τ and t_{lpf} values used in calculating CF were estimated from exponential fits made to the values obtained in the RH bins (similar as for τ in Mammarella et al., 2009). In the case of CO_{2} and H_{2}O fluxes, the power spectra were subjected to noise removal prior to utilising them in Eq. (13) by assuming that the signal was contaminated by white noise and that at the highest frequencies the power spectra contained only noise. See the influence of different noise removal techniques on the estimation of highfrequency attenuation in our companion paper (Aslan et al., 2021).
Often both H_{e,cs} and H_{e,ps} are approximated by H with a value for τ that describes the lowpass filtering of the EC setup according to Eq. (5). Then the measured flux is corrected for lowpass filtering by multiplying it with CF (see Eq. 1). There are different ways to estimate Co^{w,χ} for Eq. (1). Here we limited the analysis only to periods when the absolute sensible heat fluxes were higher than 15 W m^{−2} and estimated Co^{w,χ} in Eq. (1) with the measured Co^{w,T} following Fratini et al. (2012).
Sabbatini et al. (2018)Aubinet et al. (2000)Fratini et al. (2012)We used four methods to estimate EC system response times (and additionally t_{lpf} in Method 4) and then to calculate CF (see Table 1). All the four methods are subject to similar limitations; i.e. they are applicable only when the attenuation is not strong at frequencies close to the cospectral peak frequency. Method 1 follows the ICOS EC data processing protocol (Sabbatini et al., 2018) and is implemented, for example, in EddyPro based on Hunt et al. (2016) (see also https://www.licor.com/env/support/EddyPro/topics/whatsnew.html, last access: 24 June 2021). Method 2 follows Aubinet et al. (2000) and is implemented in EddyUH (Mammarella et al., 2016). Method 3 follows Fratini et al. (2012) and is implemented in older versions of EddyPro, and Method 4 is the only one that explicitly accounts for the temporal lag caused by the lowpass filtering action of the system and is introduced in this study. In Method 4, both parameters τ and t_{lpf}, were estimated by fitting HH_{p} to H_{e,cs}. Throughout the study, crosscovariance maximisation was used to account for the lag between scalar and vertical wind speed as typically done in the global flux measurement network. Note that also other techniques have been used to estimate the time lag (e.g. Hunt et al., 2016; Taipale et al., 2010) and that the selection of the timelag estimation method may have an effect on the correct shape of the cospectral transfer function (see Sect. 2). For instance, the crosscovariance moving average method to estimate the time lag introduced in Taipale et al. (2010) should be related to Method 4 as it is also based on crosscovariance maximisation, whereas for the approach in Hunt et al. (2016) the correct cospectral transfer function is H (with the caveat that the physical time lag is estimated accurately).
3.2.2 Lowpass filtering of temperature data
In order to evaluate the performance of the different spectral correction methods presented in Sect. 3.2.1, highfrequency T time series were attenuated with different values of τ in order to mimic attenuated scalar time series. Similar to the procedure used by Aslan et al. (2021), T time series were converted to frequency domain via a Fourier transform, multiplied by Eq. (3) and converted back to time domain with the inverse Fourier transform. Three different values for τ were used: 0.1, 0.3 and 0.7 s. The correct value for CF was obtained from the ratio between $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{T}^{\prime}}$ and $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{\widehat{T}}^{\prime}}$, where $\widehat{T}$ was the time series attenuated with the methodology described above. Note that here $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{\widehat{T}}^{\prime}}$ was calculated using crosscovariance maximisation in order to mimic regular EC gas flux processing, meaning that the lag caused by lowpass filtering was taken into account. This correct value for CF was then compared against CF estimated with the four methods described in Sect. 3.2.1 in order to evaluate the accuracy of the methods.
4.1 Timelag dependence on attenuation and turbulence timescales
In Sect. 2.2, we explained why the following approximate relationship must hold up to some bound at the highfrequency end: ${H}_{\mathrm{p}}\approx \mathrm{1}/\sqrt{H}$. Using the assumption that t_{lpf}=Cτ, this equation was solved numerically for the constant C as a function of frequency f. Figure 1 (upper panel) illustrates that such a solution is dependent on frequency and varies over a relatively large interval up to the frequency $\frac{\mathrm{2}}{\mathrm{3}}\frac{\mathit{\pi}}{\mathit{\tau}}$ (the relationship holds well up to such a frequency for any τ). This suggests that t_{lpf} and τ are not directly proportional, but their relation depends on frequency. In accordance with the approximations with the Taylor expansion (Sect. 2.2), at low frequencies the constant C equals 1. For illustration, the lower panel of Fig. 1 presents the left and right side terms of the approximate equality (${H}_{\mathrm{p}}\approx \mathrm{1}/\sqrt{H}$) calculated using a value of 0.63 for the coefficient C. In spite of the large variation of the exact coefficient with frequency shown in the upper panel in Fig. 1, the selected constant value produces very good correspondence for the whole frequency range, suggesting that C can be approximated to be constant, although strictly speaking it is not. The value of C=0.63 was obtained by estimating the average time lag in a least square sense over equally spaced frequencies in logarithmic scale over the interval from $\mathit{\tau}=\mathrm{0.01}\frac{\mathit{\pi}}{\mathit{\tau}}$ to $\frac{\mathrm{2}}{\mathrm{3}}\frac{\mathit{\pi}}{\mathit{\tau}}$ for a range of time constants from 0.05 to 1 s. We observed that the optimum coefficient was independent of τ and equal to 0.63 (with twodigit accuracy). Note that these analyses were independent of the actual turbulent signal as they were obtained using the transfer functions only (H and H_{p}).
In order to evaluate the influence of turbulent signal on the dependence between t_{lpf} and τ, the integral in Sect. 2.2, i.e. Eq. (9), was solved numerically with various combinations of time lag, τ and $(zd)/U$ in order to evaluate which timelag value resulted in maximum value for the integral for a given τ and $(zd)/U$ combination. This approach is similar to finding the time lag between signals using crosscovariance maximisation, and it identified the dependence of the lowpassfilteringinduced lag (t_{lpf}) on τ and the turbulence timescale ($(zd)/U$). For this numerical experiment, τ varied between 0.05 and 1 s, $(zd)/U$ between 1 and 10 s, and Co was calculated based on a model (similar to the one in Kristensen et al., 1997, but fitted to Siikaneva observations). Based on this analysis, t_{lpf} depended approximately linearly on τ (Fig. 2) in accordance with the analysis above, yet the dependence had an additional weak $\frac{\mathit{\tau}U}{zd}$ term: ${t}_{\text{lpf}}=a\mathit{\tau}b\frac{\mathit{\tau}U}{zd}+c$, where a=0.73, b=0.15 s and c=0.02 s. However, this numerical experiment was repeated also for very highattenuation levels; the dependence found above failed to fully recover t_{lpf} values at high attenuation, and this bias increased with τ and inverse of turbulence timescale (i.e. $U/(zd)$). For example at τ=2 s and $(zd)/U=\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$ the equation above gave t_{lpf}=1.2 s, whereas based on the numerical integration value 1.0 s was obtained, meaning that the dependence above was not able to accurately describe the lowpassinduced time lag at these highattenuation levels. This is likely due to the fact that also the shape of the cospectrum Co in Eq. (9) has an effect on t_{lpf} dependence on τ. At very highattenuation levels the peak of the cospectrum is also attenuated, and this results in a different t_{lpf} dependence on τ than the equation above which describes the dependence when attenuation takes place at high frequencies (i.e. inertial subrange).
The discrepancy between the two dependencies between t_{lpf} and τ obtained above was likely due to the fact that the latter took into account also the variability in the turbulent flux with frequency, whereas the former was based purely on transfer functions. Nevertheless, these results indicate that, for a given site, t_{lpf} can be approximated to be constant at a specific attenuation (τ) level, hence supporting the estimation of t_{lpf} as a fitting parameter in Method 4 when fitting HH_{p} to H_{e,cs} (see Sect. 3.2.1 and Table 1).
4.2 Assessment of the flux loss correction methods using an attenuated T time series
The four methods to estimate the CF (Table 1) were evaluated using an artificially attenuated turbulent T time series comparing different values for τ. Figure 3 shows example crosscovariance functions from the Siikaneva site with three different values of τ applied to the same 2 h time series.
The peak of the crosscovariance between w and $\widehat{T}$ shifted to longer positive times (i.e. $\widehat{T}$ lags w), and the crosscovariance maximum decreased as τ increased. This is related to the phase shift and flux attenuation caused by the lowpass filter, respectively. The lowpass filtering was also evident in the corresponding power spectra and cospectra. The resulting empirical transfer functions in Fig. 4 show that power spectra were more attenuated than the cospectra when cospectra were calculated from data in which the time lag derived from crosscovariance maximisation was used to shift the time series (as is typically done when processing EC data). Moreover, the shape of the empirical transfer function derived from the cospectra differed from Eq. (5). This difference in shape was related to H_{p} and the phase shift caused by the lowpass filter. Note that the phase shift resulted in negative values for HH_{p} and H_{e,cs} at high frequencies, meaning that the related attenuated cospectrum had both positive and negative values. This implies that attenuated cospectra should not be presented on log–log scales (negative values cannot be presented on log scale). Note also that HH_{p} is not exactly equal to $\sqrt{H}$ at all frequencies. In particular, HH_{p} shows negative values at high frequencies, whereas $\sqrt{H}$ does not. Hence the approximation of HH_{p} with $\sqrt{H}$ will likely fail under very strong attenuation. When the correction factor was calculated with Method 4, which takes into account the lowpassfilteringinduced phase shift, the attenuation factors ($\mathrm{1}/\text{CF}$) agree with the normalised crosscovariance maxima (see Fig. 3). This indicates that the attenuation of crosscovariance maxima was accurately estimated with this method. However, when H_{p} was neglected and the attenuation factors were calculated with H only, then they agreed with crosscovariance values at zero lag like those predicted by theory (Sect. 2). Note that the time lag discussed here only represents the lowpassfilter effect and not other effects that also alter the time lag. This is only possible because we work with degraded temperature time series, in which phase shifts through lowpass filtering are the only cause for time delays.
The value for τ used in attenuating the T time series was further varied between 0.05 and 1 s and used for the same 2 h time series. The bias in CF calculated with Method 1 increased with τ (Fig. 5a). The other methods did not produce as significantly biased estimates for CF when using this example dataset, although methods 2 and 3 somewhat underestimated CF (approximately 4 % underestimation when τ = 0.9 s), and this underestimation increased with τ. These findings are in agreement with Fig. 4: $\sqrt{H}$ was above and H below the empirical transfer function calculated from the cospectrum (H_{e,cs}), meaning that they underestimated and overestimated the attenuation, respectively. These findings are also in accordance with the results of Aslan et al. (2021) with noisy measurements. Note that Method 2 estimated close to correct value for CF (Fig. 5a) despite biased estimates for the response time (Fig. 5b) due to the compensation of two errors (biased τ and incorrect shape for the cospectral transfer function). Method 4 was able to reproduce the correct value for τ from H_{e,cs} (Fig. 5), indicating that the difference between H_{e,cs} and H_{e,ps} was indeed related to H_{p} and the lowpassfilteringrelated phase shift. The time lags estimated by Method 4 (i.e. by fitting HH_{p} to H_{e,cs}) agreed with the ones estimated by maximising the crosscovariance (Fig. 5c), also corroborating this conclusion. It should be noted that the step changes in the data that are very visible in Fig. 5c but to some extent also in Fig. 5a and b are caused by the finite resolution of the underlying data of w and T of 0.1 s.
In order to evaluate the four methods to estimate CF further, turbulent T data were lowpass filtered with three values of τ (0.1, 0.3 and 0.7 s) over several summer months (May–August) at two contrasting EC sites, Hyytiälä and Siikaneva. Then the related flux attenuation was corrected with the four methods, and the results were compared against a reference flux calculated from unattenuated T data. Table 2 summarises the relative differences between the four correction methods and the reference. In general, the values obtained with this extended dataset are in line with Fig. 5a for Siikaneva, although slight differences can be noted. The performance of all four methods decreases with increasing τ at both sites, and Method 1 showed the worst performance, which is in line with the example shown in Fig. 5 and follows the predictions from theory (Sect. 2). Of the analysed cases, the biggest bias (+8.6 %) was found when applying Method 1 at Siikaneva with τ = 0.7 s. In general, the performance of the methods was worse at Siikaneva than at Hyytiälä. This was likely related to the lower measurement height at Siikaneva (z−d = 2.7 m) than at Hyytiälä (z−d = 13 m) and thus the comparably larger contribution of high frequencies to the covariance. At Siikaneva, smaller eddies dominated the turbulent transport, and hence the effect of inaccuracies in describing the highfrequency attenuation of the signal was amplified compared with Hyytiälä.
The importance of turbulence scale on the performance of the four methods was evaluated by stratifying the data according to stability and plotting them against $\mathit{\tau}/\text{ITS}$ (Fig. 6), where ITS is the integral timescale of ${w}^{\prime}{T}^{\prime}$ calculated from the autocovariance of ${w}^{\prime}{T}^{\prime}$ time series by integrating the normalised autocovariance function until its first zero crossing (e.g. Kaimal and Finnigan, 1994). This ratio describes the ratio of attenuation and turbulent transport scales. It is also worth noting that this ratio can alternatively be approximated by the ratio of the cospectral peak frequency (${f}_{m}={n}_{\mathrm{m}}U/(zd)$, where n_{m} is the normalised cospectral peak frequency) to the cutoff frequency (${f}_{\text{co}}=\mathrm{1}/\left(\mathrm{2}\mathit{\pi}\mathit{\tau}\right)$) (Rannik et al., 2016). High values for $\mathit{\tau}/\text{ITS}$ describe periods when the attenuation timescale is large relative to the timescale of turbulent scalar transport, whereas low values correspond to cases when the attenuation timescale is small relative to this turbulence scale. When CF was estimated with methods 2 or 3, the relative bias in CF was typically within ± 3 % and increased linearly with $\mathit{\tau}/\text{ITS}$ at Siikaneva. The different behaviour of the two sites in Fig. 6b and c may be related to differences in spectral characteristics of turbulent transport at these two sites since measurements in Hyytiälä were made above forest (roughness sublayer flow) and in Siikaneva above short vegetation (boundary layer flow). For Method 4 the bias was within ± 2 % at both sites and independent of the $\mathit{\tau}/\text{ITS}$ ratio. The remaining small bias obtained with Method 4 can be speculated to be related to the quadrature spectra, which was assumed to be zero, although strictly speaking it is not (not shown). For Method 1 the relative bias in CF scaled with $\mathit{\tau}/\text{ITS}$ and, after normalisation with ITS, the data from both field sites collapsed into one common relationship (Fig. 6a). Horst (1997) derived a simple formula to estimate CF for different situations as a function of $(\mathit{\tau}/\text{ITS}{)}^{\mathit{\alpha}}$, where $\mathit{\alpha}=\mathrm{7}/\mathrm{8}$ in neutral and unstable cases and α=1 in stable cases. Similar dependencies could be derived for the bias in CF derived with Method 1. However, these dependencies should not be used to rectify the biased CF values, but rather CF should be estimated with methods that do not result in biased CF values in the first place.
Since the relative bias in CF when estimated with Method 1 scales with $\mathit{\tau}/\text{ITS}$, this dependence can be used to assess how much scalar fluxes are biased with different EC setups when Method 1 is utilised in the EC data processing. For the following example calculations we approximate a linear dependence between ITS and $(zd)/U$ and used the fit to estimate ITS for different z−d and U combinations. For a short tower (z−d = 1.5 m) with moderate wind speed (U = 2 m s^{−1}) and attenuation (τ = 0.2 s), fluxes are biased by 4 % based on the dependency shown in Fig. 6a. For the same tower, higher attenuation (τ = 0.8 s) results in larger bias (9 %) if Method 1 is utilised. For a tall tower the biases are not as significant since larger eddies dominate the transport, and hence $\mathit{\tau}/\text{ITS}$ is smaller for a given τ than in the case of a shorter tower. For a tall tower (z−d = 20 m) with moderate wind speed (U = 3 m s^{−1}) and attenuation (τ = 0.2 s) the bias is 1 %–2 %, and with higher attenuation (τ = 0.8 s) the bias is 3 %–4 %. These example calculations demonstrate the magnitude of the flux bias resulting from biased spectral corrections at contrasting EC sites.
4.3 CO_{2} and H_{2}O fluxes corrected with the four methods
The applicability of the results acquired with attenuated T was analysed by processing CO_{2} and H_{2}O fluxes from summer months at the two sites (Sect. 3.1). The gas fluxes were corrected with the four spectral correction methods in order to assess the systematic biases stemming from biased spectral corrections. For CO_{2} the fitting of H to H_{e,ps}, of H to H_{e,cs} and of HH_{p} to H_{e,cs} estimated response times of 0.05, 0.10 and 0.13 s at Hyytiälä and of 0.13, 0.12 and 0.18 s at Siikaneva, respectively. The lowpassinduced time lag (t_{lpf}) estimated as a fitting parameter when fitting HH_{p} to H_{e,cs} yielded close to 0.1 s at both sites for CO_{2}. For reference, the time lags estimated via crosscovariance maximisation which are combinations of signal travel times in the sampling line (t_{phys}) and t_{lpf} yielded on average 1.4 and 0.2 s at Siikaneva and Hyytiälä, respectively. For Hyytiälä the signal travel time calculated from the sampling tube dimensions and flow rate was 0.06 s. However, this calculation neglects the additional time lag caused by the LI7550 interface unit and hence can be assumed to slightly underestimate t_{phys}.
Based on theoretical considerations (Sect. 2) and results obtained with attenuated T time series, H fit to H_{e,cs} (Method 2) was projected to estimate smaller response times than the other two methods that were projected to agree (see, for example, Fig. 5b). Instead, here for CO_{2} the HH_{p} fit to H_{e,cs} gave longer response times than the other two methods (Fig. 7). For signals with small attenuation both τ and t_{lpf} are likely difficult to estimate accurately by fitting HH_{p} to H_{e,cs}. The lowpassfilteringinduced lag (t_{lpf}) can attain values only with the temporal resolution of the data being processed (here 0.1 s), and hence for time series exhibiting small attenuation t_{lpf} may well be zero. This would cause H_{e,ps} and H_{e,cs} to be similar, and hence they would result in similar estimates for the response time, as observed here for CO_{2}. This granular estimation of t_{lpf} might explain why the response times for CO_{2} did not follow the expectations. It should be also noted that the methods to estimate the response time have a limited accuracy of estimation stemming from the algorithms used in the estimation, instrument limitations and meteorological conditions under which the observations are made (Aslan et al., 2021). Furthermore, the empirical transfer function derived from cospectra (H_{e,cs}) contains the attenuation related also to sensor separation (Horst and Lenschow, 2009), whereas the empirical transfer function estimated from power spectra (H_{e,ps}) is related only to gasanalyserinduced lowpass filtering.
For H_{2}O, τ depends on RH (Mammarella et al., 2009; Ibrom et al., 2007a; Runkle et al., 2012), and hence the analysis was done in RH bins. Figure 8 shows an example of empirically derived transfer functions for H_{2}O at Siikaneva in moderate relative humidity conditions. H_{e,cs} showed less highfrequency attenuation than H_{e,ps} due to crosscovariance maximisation as predicted by the theory (Sect. 2). HH_{p} fit to H_{e,cs} resulted in a similar value for τ as the one obtained by fitting H to H_{e,ps}. The other fitting parameter in HH_{p} fit (t_{lpf}) was 0.14 s, indicating the lowpass filtering of H_{2}O time series caused an additional lag to the H_{2}O time series. Note that H does not perfectly describe the shape of H_{e,ps}. This is due to adsorption and desorption of H_{2}O to sampling tube walls, a process which cannot be described accurately with Eq. (3) (Nordbo et al., 2013). Hence different transfer functions have been proposed for H_{2}O (De Ligne et al., 2010; Nordbo et al., 2014).
Figure 9 shows results obtained for H_{2}O in different RH bins at the two sites. The H fit to H_{e,cs} (Method 2) gave smaller values for τ in agreement with the results obtained with attenuated T (Fig. 5b) and expectations based on the theory (Sect. 2). The HH_{p} fit to H_{e,cs} (Method 4) was able to reproduce the H_{2}O timelag dependency on RH, yet it estimated systematically smaller values for the timelag difference between H_{2}O and CO_{2} than those obtained with crosscovariance maximisation. The timelag difference between H_{2}O and CO_{2} can be assumed to be related to t_{lpf} of H_{2}O since H_{2}O is more attenuated than CO_{2}. The difference between H_{2}O and CO_{2} time lags obtained with crosscovariance maximisation showed linear dependence on estimated response times in relative humidity bins with a similar dependency at both sites ($y=ax+b$, where a=0.47 and b = −0.05 s for Siikaneva and a=0.44 and b = −0.07 s for Hyytiälä, and where y is the timelag difference and x the response time. Fits were based on orthogonal linear regression). Difference between H_{2}O and CO_{2} time lags resulted from stronger lowpass filtering of H_{2}O than CO_{2} signal, and this analysis assumes that the mere travel times through the sampling lines (i.e. t_{phys}), which do not depend on the lowpass filtering effects, were the same for H_{2}O and CO_{2}. These results are in line with Fig. 2 where it was shown that the lowpassfilteringinduced time lag was primarily determined by a linear dependence on τ with a secondary dependence on turbulence timescale. However, the slopes obtained above were different from the ones derived in Sect. 4.1.
The CF values estimated with the four methods agreed at Hyytiälä for low RH periods, but they diverged at high RH periods. For instance when RH > 80 % the four methods gave on average 1.39, 1.29, 1.23 and 1.37 at the Hyytiälä site, respectively. These results are in line with the ones obtained with attenuated T time series (Sect. 4.2 and Fig. 5a): when τ was large, Method 1 overestimated and Method 3 underestimated CF, similar as here for H_{2}O at high RH (and hence large τ). In contrast, for Siikaneva Method 1 gave systematically higher values for the correction factor regardless of RH, whereas the other three methods gave more similar CF RH dependency. This difference to Hyytiälä was likely due to lower τ values also at high RH (Fig. 9). Note also that the EC measurement setup (measurement height, gas sampling system) was not identical at the two sites.
On average, the differences between the four different methods to calculate CF were small, typically within ± 3 % at both sites for both gases (CO_{2} and H_{2}O) (Table 3). The biggest mean relative difference (4.1 %) was found for H_{2}O measurements at the Siikaneva site for Method 1. This was likely due to the combination of low measurement height and high attenuation which resulted in biased CF values, in accordance with the findings obtained with the attenuated T time series (Fig. 6a). Interestingly, for CO_{2} the biggest difference at both sites was obtained with Method 3, yielding smaller fluxes than the reference (Method 4). The difference was bigger than the one obtained with Method 1 which contradicts the findings presented above with attenuated T (Sect. 4.2 and Table 2) and expectations based on theory (Sect. 2), although the differences found between the methods were generally small. Note that if t_{lpf}≈0, then both power spectra and cospectra are attenuated similarly (i.e. with H) resulting in methods 1, 2 and 4 being similar. Infrared gas analysers for CO_{2} and H_{2}O can often be mounted close to the inlet, and in both setups considered here inlet lines are fairly short. The measurement of other gases, including CH_{4} and N_{2}O, usually requires larger equipment which often has to be operated on the ground and/or some distance away from the mast, requiring longer inlet lines with slow time response. In these situations the correct shape of the cospectral transfer function matters a lot more than for the measurement setups included here, and hence the variability between methods is expected to be larger.
4.4 Correct form for cospectral transfer function
There has been a longstanding debate about what the correct form for cospectral transfer function is when the scalar measurements are done with a firstorder linear sensor and vertical wind speed is not attenuated. The seminal paper on EC frequencyresponse corrections by Moore (1986) described the flux attenuation related to the scalar sensor with $\sqrt{H}$ (see Eq. 5), which was later deemed erroneous, for example, by Eugster and Senn (1995), Horst (1997, 2000), and others. Partly based on similar derivations as shown here in Sect. 2, the correct form for the cospectral transfer function was argued by Horst (1997, 2000) and Massman (2000) to be H instead of $\sqrt{H}$ . Lately, both forms (H and $\sqrt{H}$) have been utilised in the literature (e.g. Fratini et al., 2012; Foken et al., 2012b; Mammarella et al., 2016; Hunt et al., 2016; Wintjen et al., 2020) without clear consensus on which one is correct. Part of this confusion is likely related to the fact that some of the past studies have been using incorrect terminology by using the term cospectral transfer function (Eq. 5) even though they have derived/utilised a transfer function for the amplitude spectrum (Eq. 10).
Briefly summarising the findings in this study, it was shown in Sect. 2 in accordance with, for example, Horst (2000) that in the case the travel time in the gas sampling system could be accurately estimated and accounted for, the attenuation of cospectrum (and hence flux) could be described with H. However, crosscovariance maximisation inadvertently accounts also for the scalar lowpassfilteringrelated time shift (t_{lpf}) and consequently induces a phase shift between the lagcorrected attenuated scalar and vertical wind time series. Hence in the case that crosscovariance maximisation is used to align the two time series, the correct form for cospectral transfer function is HH_{p}, in which the part related to the phase shift (H_{p}) is calculated with t_{lpf}. It was shown above that HH_{p} can be approximated with $\sqrt{H}$. Therefore the debate about which is the correct form for cospectral transfer function (H or $\sqrt{H}$) relates to the lowpassfilteringinduced phase shift as discussed already by Eugster and Senn (1995) and Horst (2000), albeit with a different conclusion. Note, however, that $\sqrt{H}$ is an approximation and not an exact representation of the cospectral transfer function.
The influence of lowpassfilteringinduced phase shift on estimation of highfrequency response of an EC setup was analysed. The analysis assumed that the EC setup consisted of a fastresponse anemometer and a linear, firstorderresponse scalar sensor. Spectral corrections aiming at correcting the EC fluxes for lowpass filtering were estimated with four methods: three widely used methods and one newly proposed method which specifically accounts for the interaction between the lowpassfilteringinduced phase shift and attenuation. Based on theoretical considerations and experimental results we come to the following conclusions:

Crosscovariance maximisation overestimates the signal travel time in the EC sampling line since it inadvertently accounts also for the lag caused by lowpass filtering of scalar time series caused by the nonideal measurement system. The bias in the estimated time lag depends linearly on lowpassfilter response time (τ) with a small additional dependence on turbulence timescale. Hence, other means for estimating the physical time lag might be warranted, especially in the case of noisy measurements (Langford et al., 2015). Further research on this topic is required.

Both power spectra and cospectra are attenuated with the same transfer function (H, as noted also by Horst, 2000) in the case that the travel time of the scalar signal in the sampling line can be accurately estimated. However, if crosscovariance maximisation is used, then attenuation of cospectra follows HH_{p} in which H_{p} accounts for the bias in the estimated travel time caused by crosscovariance maximisation, and it was shown that HH_{p} can be roughly approximated with $\sqrt{H}$. Therefore, all flux calculation algorithms which rely on crosscovariance maximisation and at the same time use H to describe the attenuation of the cospectra (e.g. Horst, 1997; Massman, 2000) can be projected to be biased in the same way as Method 1 used here if the response time (i.e. τ) is accurately estimated (e.g. from power spectra). In these cases the bias in EC fluxes can be 10 % or above for short tower measurements with moderate attenuation. Hence, data reprocessing is warranted in order to rectify this bias.

In order to estimate and correct for the flux attenuation correctly, it is vital to accurately describe the attenuation of the cospectra in the correction procedure. Hence, while fitting H to cospectra (Method 2) biased the firstorder response time due to the use of the incorrect cospectral transfer function, it resulted only in a small bias for the flux correction factor since the method describes the attenuation of cospectra accurately once the crosscovariance maximisation has been applied. By contrast, fitting of H to power spectra (Method 1) correctly estimates the response time, but it nevertheless yields biased flux corrections since it utilised a transfer function that is incorrect when estimating the flux using crosscovariance maximisation.

The theoretical framework proposed in this paper was able to describe the changes in H_{2}O time lag as H_{2}O response time increased with relative humidity. This suggests that the findings derived here with attenuated temperature time series are applicable also in real world situations for EC gas flux measurements. However, results for CO_{2} deviated from the expectations based on the theory. Reasons for this remained unclear; however, lowpass filtering of inert gases (such as CO_{2}) in the gas sampling line may deviate from what was described here in such a way that the lowpass filter acting on CO_{2} may not induce a phase shift. This should be investigated by using CO_{2} flux data from an EC setup with pronounced attenuation of CO_{2} fluctuations in the sampling line or alternatively in a controlled environment in a laboratory (Metzger et al., 2016; Aubinet et al., 2016).
In summary, it is suggested that the spectral correction methods implemented in EC data processing software are revised so that the influence of lowpassfilteringinduced phase shift is recognised following the findings presented above. In particular, the usage of Method 1 is discouraged as it resulted in clearly biased flux values.
Here we derive Eq. (4) starting from the crossspectrum (Eq. 2):
where h_{χ} can be described with Eq. (3), which gives
where ${Z}_{\mathrm{w}}{Z}_{\mathit{\chi}}^{*}=\text{Co}+jQ$ was utilised. Next the quadrature spectrum (Q) was assumed to be zero (Horst, 1997; Massman, 2000) and hence
Now we can use Euler formula (${e}^{j{\mathit{\varphi}}_{\text{phys}}}=\mathrm{cos}{\mathit{\varphi}}_{\text{phys}}j\mathrm{sin}{\mathit{\varphi}}_{\text{phys}}$):
Data and MATLAB codes to reproduce Figs. 3, 4 and 5, in addition to time series processed with the four methods (Table 1) at the two sites, have been uploaded to the open repository Zenodo (https://doi.org/10.5281/zenodo.4769420, Peltola et al., 2021).
OP devised the original concept for the study, and all the authors provided further ideas. OP did the data analysis, with input also from ÜR. OP wrote the first draft of the manuscript, with contributions also from TA, ÜR and IM. All the authors commented on the first draft and made improvements.
The authors declare that they have no conflict of interest.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This study was supported by ICOS and 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 funded by the Academy of Finland, and Eiko Nemitz acknowledges support by the Natural Environment Research Council as part of the UKSCAPE programme delivering national capability.
This research has been supported by the European Commission, H2020 Research Infrastructures (RINGO; grant no. 730944), the Väisälän Rahasto (grant to Toprak Aslan), the Academy of Finland (grant no. 315424) and the Natural Environment Research Council (award no. NE/R016429/1).
This paper was edited by Glenn Wolfe and reviewed by Johannes Laubach, Marc Aubinet, and one anonymous referee.
Aslan, T., Peltola, O., Ibrom, A., Nemitz, E., Rannik, Ü., and Mammarella, I.: The highfrequency response correction of eddy covariance fluxes – Part 2: An experimental approach for analysing noisy measurements of small fluxes, Atmos. Meas. Tech., 14, 5089–5106, https://doi.org/10.5194/amt1450892021, 2021. a, b, c, d, e
Aubinet, M., Grelle, A., Ibrom, A., Rannik, U., Moncrieff, J., Foken, T., Kowalski, A. S., Martin, P. H., Berbigier, P., Bernhofer, C., Clement, R., Elbers, J., Granier, A., Grunwald, 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, 2000. a, b, c
Aubinet, M., Vesala, T., and Papale, D.: Eddy covariance: a practical guide to measurement and data analysis, Springer Atmospheric Sciences, Springer, Dordrecht, 2012. a, b, c
Aubinet, M., Joly, L., Loustau, D., De Ligne, A., Chopin, H., Cousin, J., Chauvin, N., Decarpenterie, T., and Gross, P.: Dimensioning IRGA gas sampling systems: laboratory and field experiments, Atmos. Meas. Tech., 9, 1361–1367, https://doi.org/10.5194/amt913612016, 2016. a, b
Baldocchi, D. D.: Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: past, present and future, Glob. Change Biol., 9, 479–492, https://doi.org/10.1046/j.13652486.2003.00629.x, 2003. a
Brock, F. V.: A Nonlinear Filter to Remove Impulse Noise from Meteorological Data, J. Atmos. Ocean. Tech., 3, 51–58, https://doi.org/10.1175/15200426(1986)003<0051:ANFTRI>2.0.CO;2, 1986. a
De Ligne, A., Heinesch, B., and Aubinet, M.: New Transfer Functions for Correcting Turbulent Water Vapour Fluxes, Bound.Lay. Meteorol., 137, 205–221, https://doi.org/10.1007/s1054601095259, 2010. a
Eugster, W. and Senn, W.: A Cospectral Correction Model for Measurement of Turbulent NO_{2} Flux, Bound.Lay. Meteorol., 74, 321–340, https://doi.org/10.1007/BF00712375, 1995. a, b
Foken, T. and Wichura, B.: Tools for quality assessment of surfacebased flux measurements, Agr. Forest Meteorol., 78, 83–105, 1996. a, b
Foken, T., Aubinet, M., and Leuning, R.: The Eddy Covariance Method, chap. 1, in: Eddy Covariance, edited by: Aubinet, M., Vesala, T., and Papale, D., Springer Atmospheric Sciences, Springer Netherlands, https://doi.org/10.1007/9789400723511_1, 1–19, 2012a. a, b
Foken, T., Leuning, R., Oncley, S. R., Mauder, M., and Aubinet, M.: Corrections and Data Quality Control BT – Eddy Covariance: A Practical Guide to Measurement and Data Analysis, Springer Netherlands, Dordrecht, https://doi.org/10.1007/9789400723511_4, 85–131, 2012b. a
Franz, D., Acosta, M., Altimir, N., Arriga, N., Arrouays, D., Aubinet, M., Aurela, M., Ayres, E., LópezBallesteros, A., Barbaste, M., Berveiller, D., Biraud, S., Boukir, H., Brown, T., Brömmer, C., Buchmann, N., Burba, G., Carrara, A., Cescatti, A., Ceschia, E., Clement, R., Cremonese, E., Crill, P., Darenova, E., Dengel, S., D'Odorico, P., Filippa, G., Fleck, S., Fratini, G., Fuß, R., Gielen, B., Gogo, S., Grace, J., Graf, A., Grelle, A., Gross, P., Grönwald, T., Haapanala, S., Hehn, M., Heinesch, B., Heiskanen, J., Herbst, M., Herschlein, C., Hörtnagl, L., Hufkens, K., Ibrom, A., Jolivet, C., Joly, L., Jones, M., Kiese, R., Klemedtsson, L., Kljun, N., Klumpp, K., Kolari, P., Kolle, O., Kowalski, A., Kutsch, W., Laurila, T., De Ligne, A., Linder, S., Lindroth, A., Lohila, A., Longdoz, B., Mammarella, I., Manise, T., Jiménez, S., Matteucci, G., Mauder, M., Meier, P., Merbold, L., Mereu, S., Metzger, S., Migliavacca, M., Mölder, M., Montagnani, L., Moureaux, C., Nelson, D., Nemitz, E., Nicolini, G., Nilsson, M., De Beeck, M., Osborne, B., Löfvenius, M., Pavelka, M., Peichl, M., Peltola, O., Pihlatie, M., Pitacco, A., Pokorný, R., Pumpanen, J., Ratié, C., Rebmann, C., Roland, M., Sabbatini, S., Saby, N., Saunders, M., Schmid, H., Schrumpf, M., Sedlák, P., Ortiz, P., Siebicke, L., Šigut, L., Silvennoinen, H., Simioni, G., Skiba, U., Sonnentag, O., Soudani, K., Soulé, P., Steinbrecher, R., Tallec, T., Thimonier, A., Tuittila, E.S., Tuovinen, J.P., Vestin, P., Vincent, G., Vincke, C., Vitale, D., Waldner, P., Weslien, P., Wingate, L., Wohlfahrt, G., Zahniser, M., and Vesala, T.: Towards longTerm standardised carbon and greenhouse gas observations for monitoring Europe's terrestrial ecosystems: A review, Int. Agrophys., 32, 439–455, https://doi.org/10.1515/intag20170039, 2018. a
Fratini, G., Ibrom, A., Arriga, N., Burba, G., and Papale, D.: Relative humidity effects on water vapour fluxes measured with closedpath eddycovariance systems with short sampling lines, Agr. Forest Meteorol., 165, 53–63, https://doi.org/10.1016/j.agrformet.2012.05.018, 2012. a, b, c, d, e
Horst, T. W.: A simple formula for attenuation of eddy fluxes measured with firstorderresponse scalar sensors, Bound.Lay. Meteorol., 82, 219–233, 1997. a, b, c, d, e, f, g, h
Horst, T. W.: On Frequency Response Corrections for Eddy Covariance Flux Measurements, Bound.Lay. Meteorol., 94, 517–520, https://doi.org/10.1023/A:1002427517744, 2000. a, b, c, d, e
Horst, T. W. and Lenschow, D. H.: Attenuation of Scalar Fluxes Measured with Spatiallydisplaced Sensors, Bound.Lay. Meteorol., 130, 275–300, https://doi.org/10.1007/s1054600893480, 2009. a, b
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 wintergrazed pasture, Biogeosciences, 13, 2927–2944, https://doi.org/10.5194/bg1329272016, 2016. a, b, c, d
Ibrom, A., Dellwik, E., Flyvbjerg, H., Jensen, N. O., and Pilegaard, K.: Strong lowpass filtering effects on water vapour flux measurements with closedpath eddy correlation systems, Agr. Forest Meteorol., 147, 140–156, https://doi.org/10.1016/j.agrformet.2007.07.007, 2007a. a, b, c, d, e
Ibrom, A., Dellwik, E., Larsen, S. E., and Pilegaard, K.: On the use of the WebbPearmanLeuning theory for closedpath eddy correlation measurements, Tellus B, 59, 937–946, https://doi.org/10.1111/j.16000889.2007.00311.x, 2007b. a
Kaimal, J. C. and Finnigan, J. J.: Atmospheric boundary layer flows: their structure and measurement, Oxford University Press, New York, 1994. a
Kristensen, L., Mann, J., Oncley, S. P., and Wyngaard, J. C.: How Close is Close Enough When Measuring Scalar Fluxes with Displaced Sensors?, J. Atmos. Ocean. Tech., 14, 814–821, https://doi.org/10.1175/15200426(1997)014<0814:HCICEW>2.0.CO;2, 1997. a
Langford, B., Acton, W., Ammann, C., Valach, A., and Nemitz, E.: Eddycovariance data with low signaltonoise ratio: timelag determination, uncertainties and limit of detection, Atmos. Meas. Tech., 8, 4197–4213, https://doi.org/10.5194/amt841972015, 2015. a, b
Mammarella, I., Launiainen, S., Gronholm, T., Keronen, P., Pumpanen, J., Rannik, U., and Vesala, T.: Relative Humidity Effect on the HighFrequency Attenuation of Water Vapor Flux Measured by a ClosedPath Eddy Covariance System, J. Atmos. Ocean. Tech., 26, 1856–1866, https://doi.org/10.1175/2009JTECHA1179.1, 2009. a, b, c, d
Mammarella, I., Nordbo, A., Rannik, U., 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, https://doi.org/10.1002/2014JG002873, 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/amt949152016, 2016. a, b, c
Massman, W. J.: A simple method for estimating frequency response corrections for eddy covariance systems, Agr. Forest Meteorol., 104, 185–198, https://doi.org/10.1016/S01681923(00)001647, 2000. a, b, c, d, e, f, g, h, i, j, k, l
Massman, W. J. and Lee, X.: Eddy covariance flux corrections and uncertainties in longterm studies of carbon and energy exchanges, Agr. Forest Meteorol., 113, 121–144, 2002. a, b
Metzger, S., Burba, G., Burns, S. P., Blanken, P. D., Li, J., Luo, H., and Zulueta, R. C.: Optimization of an enclosed gas analyzer sampling system for measuring eddy covariance fluxes of H_{2}O and CO_{2}, Atmos. Meas. Tech., 9, 1341–1359, https://doi.org/10.5194/amt913412016, 2016. a, b
Moncrieff, J. B., Massheder, J. M., deBruin, H., Elbers, J., Friborg, T., Heusinkveld, B., Kabat, P., Scott, S., Soegaard, H., and Verhoef, A.: A system to measure surface fluxes of momentum, sensible heat, water vapour and carbon dioxide, J. Hydrol., 189, 589–611, 1997. a, b
Moore, C. J.: FrequencyResponse Corrections for EddyCorrelation Systems, Bound.Lay. Meteorol., 37, 17–35, 1986. a, b
Nemitz, E., Mammarella, I., Ibrom, A., Aurela, M., Burba, G., Dengel, S., Gielen, B., Grelle, A., Heinesch, B., Herbst, M., Hörtnagl, L., Klemedtsson, L., Lindroth, A., Lohila, A., McDermitt, D., Meier, P., Merbold, L., Nelson, D., Nicolini, G., Nilsson, M., Peltola, O., Rinne, J., and Zahniser, M.: Standardisation of eddycovariance flux measurements of methane and nitrous oxide, Int. Agrophys., 32, 517–549, https://doi.org/10.1515/intag20170042, 2018. a
Nordbo, A., Kekäläinen, P., Siivola, E., Lehto, R., Vesala, T., and Timonen, J.: Tube transport of water vapor with condensation and desorption, Appl. Phys. Lett., 102, 194101, https://doi.org/10.1063/1.4804639, 2013. a
Nordbo, A., Kekäläinen, P., Siivola, E., Mammarella, I., Timonen, J., and Vesala, T.: SorptionCaused Attenuation and Delay of Water Vapor Signals in EddyCovariance Sampling Tubes and Filters, J. Atmos. Ocean. Tech., 31, 2629–2649, https://doi.org/10.1175/JTECHD1400056.1, 2014. a, b, c
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/bg1037492013, 2013. a
Peltola, O., Aslan, T., Ibrom, A., Nemitz, E., Rannik, Ü., and Mammarella, I.: Dataset for “The high frequency response correction of eddy covariance fluxes. Part 1: an experimental approach and its interdependence with the timelag estimation”, Zenodo [data set], https://doi.org/10.5281/zenodo.4769420, 2021. a
Rannik, Ü. and Vesala, T.: Autoregressive filtering versus linear detrending in estimation of fluxes by the eddy covariance method, Bound.Lay. Meteorol., 91, 259–280, 1999. a
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/amt951632016, 2016. a
Rebmann, C., Aubinet, M., Schmid, H., Arriga, N., Aurela, M., Burba, G., Clement, R., De Ligne, A., Fratini, G., Gielen, B., Grace, J., Graf, A., Gross, P., Haapanala, S., Herbst, M., Hörtnagl, L., Ibrom, A., Joly, L., Kljun, N., Kolle, O., Kowalski, A., Lindroth, A., Loustau, D., Mammarella, I., Mauder, M., Merbold, L., Metzger, S., Mölder, M., Montagnani, L., Papale, D., Pavelka, M., Peichl, M., Roland, M., SerranoOrtiz, P., Siebicke, L., Steinbrecher, R., Tuovinen, J.P., Vesala, T., Wohlfahrt, G., and Franz, D.: ICOS eddy covariance fluxstation site setup: a review, Int. Agrophys., 32, 471–494, https://doi.org/10.1515/intag20170044, 2018. a
Rinne, J., Tuittila, E.S., Peltola, O., Li, X., Raivonen, M., Alekseychik, P., Haapanala, S., Pihlatie, M., Aurela, M., Mammarella, I., and Vesala, T.: Temporal Variation of Ecosystem Scale Methane Emission From a Boreal Fen in Relation to Temperature, Water Table Position, and Carbon Dioxide Fluxes, Global Biogeochem. Cy., 32, 1087–1106, https://doi.org/10.1029/2017GB005747, 2018. a
Riutta, T., Laine, J., Aurela, M., Rinne, J., Vesala, T., Laurila, T., Haapanala, S., Pihlatie, M., and Tuittila, E. S.: Spatial variation in plant community functions regulates carbon gas dynamics in a boreal fen ecosystem, Tellus B, 59, 838–852, https://doi.org/10.1111/j.16000889.2007.00302.x, 2007. a
Runkle, B. R. K., Wille, C., Gazovic, M., and Kutzbach, L.: Attenuation Correction Procedures for Water Vapour Fluxes from ClosedPath EddyCovariance Systems, Bound.Lay. Meteorol., 142, 401–423, https://doi.org/10.1007/s105460119689y, 2012. a, b
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, https://doi.org/10.1515/intag20170043, 2018. a, b, c, d
Starkenburg, D., Metzger, S., Fochesatto, G. J., Alfieri, J. G., Gens, R., Prakash, A., and Cristóbal, J.: Assessment of Despiking Methods for Turbulence Data in Micrometeorology, J. Atmos. Ocean. Tech., 33, 2001–2013, https://doi.org/10.1175/JTECHD150154.1, 2016. a
Taipale, R., Ruuskanen, T. M., and Rinne, J.: Lag time determination in DEC measurements with PTRMS, Atmos. Meas. Tech., 3, 853–862, https://doi.org/10.5194/amt38532010, 2010. a, b
Vesala, T., Kljun, N., Rannik, Ü., Rinne, J., Sogachev, A., Markkanen, T., Sabelfeld, K., Foken, T., and Leclerc, M. Y.: Flux and concentration footprint modelling: State of the art, Environ. Pollut., 152, 653–666, https://doi.org/10.1016/j.envpol.2007.06.070, 2008. a
Wintjen, P., Ammann, C., Schrader, F., and Brümmer, C.: Correcting highfrequency losses of reactive nitrogen flux measurements, Atmos. Meas. Tech., 13, 2923–2948, https://doi.org/10.5194/amt1329232020, 2020. a, b, c