The high frequency response correction of eddy covariance fluxes. Part 1: an experimental approach for analysing noisy measurements of small fluxes

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 and 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 5 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), where 10 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 (PSAA20). For CSA, we use three different methods: (1) a plain version of Lorentzian equation describing the H (CSAH ), (2) a square-root of the H (CSAH ), and (3) a square-root of the H with shifted vertical wind velocity time series via cross-covariance maximisation (CSAH,sync). PSAI07 tends to overestimate the time constant when low-pass filtering is low, whilst the new PSAA20 successfully estimates 15 the expected time constant regardless of the degree of attenuation and SNR. CSAH underestimates the time constant with decreasing accuracy as attenuation increases due to the omission of the quadrature spectrum. CSAH overestimates, but its accuracy increases with time-lag correction in the CSAH,sync. 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 of ±2%. 20 PSAA20 showed a similar variation, yet slightly better accuracy. CSAH underestimated fluxes by up to 4%, while CSAH overestimated them by up to 3%, a bias which was mostly eliminated with time-lag correction in the CSAH,sync (−2% to 1%). The accuracies of both PSA and CSA methods were not affected by SNR level, instilling confidence in EC flux 1 https://doi.org/10.5194/amt-2020-478 Preprint. Discussion started: 23 December 2020 c © Author(s) 2020. CC BY 4.0 License.

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 (PSA A20 ). For CSA, we use three different methods: (1) a plain version of Lorentzian equation describing the H (CSA H ), (2) a square-root of the H (CSA √ H ), and (3) a square-root of the H with shifted vertical wind velocity time series via cross-covariance maximisation (CSA √ H,sync ). PSA I07 tends to overestimate the time constant when low-pass filtering is low, whilst the new PSA A20 successfully estimates 15 the expected time constant regardless of the degree of attenuation and SNR. CSA H underestimates the time constant with decreasing accuracy as attenuation increases due to the omission of the quadrature spectrum. CSA √ H overestimates, but its accuracy increases with time-lag correction in the CSA √ H,sync . 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 P SA I07 show a bias of ±2%. 20 P SA A20 showed a similar variation, yet slightly better accuracy. CSA H underestimated fluxes by up to 4%, while CSA √ H overestimated them by up to 3%, a bias which was mostly eliminated with time-lag correction in the CSA √ H,sync (−2% to 1%). The accuracies of both PSA and CSA methods were not affected by SNR level, instilling confidence in EC flux 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).
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. 30 (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 to describe sensor related attenuation independently. This should in principle provides a better estimation of the time constant of the gas analysis system, because the spectral data is not mingled with other components, such as, e.g., sensor separation. In this approach, the effect of sensor separation is then treated explicitly with an additional correction step (Horst and Lenschow, 2009). 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 5 uncertainty in the PSA approach. This paper explores this uncertainty and proposes a novel, more robust approach to account for the noise in PSA.
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 be disappear from the measured signal. On the other hand, the use of the CSA relies on the correct determi-10 nation of the time lag between w and χ, which may be difficult in 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., 2020).
EC measurements conducted under low flux conditions result in relatively higher noise in signal, i.e. low signal-to-noise 15 ratio (SNR) (Smeets et al., 2009). EC fluxes with low SNR are typically found in many ecosystems especially for methane (CH 4 ) and nitrous oxide (N 2 O) as well as for other gases and aerosol particles, in which any implementation of the FRC becomes uncertain Nemitz et al., 2018;Oosterwijk et al., 2018). Low SNR can also be observed for CO 2 at specific ecosystems (e.g. lakes), for CH 4 over well-drained soils or peatlands during winter, for N 2 O in long periods outside the high emission periods (e.g. fertilizer applications, freeze-thaw cycles or rain events), all of which, although small, 20 ultimately contribute to the final flux estimates, and, hence must be corrected from systematic bias. Thus, investigation of uncertainties in commonly used FRC methods is of great importance in order to obtain 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 have 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 25 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 attenuated gas concentration dataset (e.g. CH 4 , N 2 O) with known characteristics. 30 We use a first order low-pass filter, which solely depends on a single time constant (τ LP F ), 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. Firstly, we try to retrieve the system time constants using the PSA and CSA, then compare those with original values (i.e., τ LP F ). Secondly, in order to demonstrate how variation in time constant estimation further affects the cumulative fluxes, we low-pass filter a one month T time series, and correct the attenuation with 35 3 https://doi.org/10.5194/amt-2020-478 Preprint. 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 (F corr ): w χ corr = w χ meas F corr . (1) One way to calculate F corr is to estimate the ratio of the integrated unbiased and biased cospectra as a function of frequency. In order to define the co-spectrum of the unattenuated scalar under the assumption that the normalized cospectrum of all scalars 10 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 Moore, 1986). The attenuated cospectra are obtained by convolution of the reference cospectra with the transfer function (H), which characterizes the filtering of the EC system.
Another way to calculate F corr is to simulate the attenuation with a recursive filter in time, rather than in frequency space, 15 in which F corr is defined as the ratio of the unattenuated and attenuated co-variances (see Goulden et al., 1997). In addition, based on the same approach, an experimental method was proposed by Ibrom et al. (2007a) to parameterise 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 F corr can be calculated. 20 In 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 on the low frequency part of the spectra, and to fulfill 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 25 et al., 2007a).
For the PSA, H can be calculated using the power spectra of χ and T (Eq. (2)), where the effect of sensor separation should additionally be treated via the method proposed by Horst and Lenschow (2009). For PSA, H is derived as: where S χ indicates the power spectrum of measured target gas mixing ratio, S T is the power spectrum of T and variances 30 (σ 2 χ and σ 2 T ) 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 σ s . Where instrumental noise is large, the power spectrum asymptotically approaches a line of slope +1 at the high frequency end (Fig. 1) which represents the fact that white noise affects all frequencies equally, but is multiplied by f in this representation (see Sect. 2.2). Thus, prior to calculation of Eq.
(2), the noise contribution to the power spectra of S χ should be removed and this is typically done by fitting a line to the high frequency end (Ibrom et al., 2007a). Finally, the frequency dependence of H P SA can be be described through a sigmoidal 5 curve which is characterised by the time constant (τ ) of the measurement system (see for more details Peltola et al., 2020, and references therein): For PSA, Ibrom et al. (2007a) updated Eq. (2) by introducing a normalisation factor, which is used to secure the spectral similarity especially for small fluxes. As described in their Eq. (6), the time constant and the normalisation factor (F n ) are 10 obtained via fitting the following equation to the dampened and noise free χ data: In our study, we follow the same procedure (hereafter PSA I07 ) for the time constant calculation for the PSA, which is summarized in Sect. 3.2. In addition to PSA I07 , we used a new comprehensive method for PSA, which is summarized in Sect.

15
Alternatively, for the CSA, H is calculated as where f is the natural frequency, Co wχ indicates the cospectrum of measured w and target gas mixing ratio χ; Co wT is the cospectrum of measured kinematic heat flux, w χ and w T are covariances calculated across a frequency range where the cospectra are not attenuated. 20 For CSA, τ is obtained similar to PSA via fitting H emp to Eq. (5), however, there is an ongoing debate on whether the correct transfer function for the cospectra would be H emp , instead of H emp (Moore, 1986;Eugster and Senn, 1995;Horst, 1997Horst, , 2000Fratini et al., 2012;Hunt et al., 2016) as discussed throughouly by Peltola et al. (2020) who also provide a mathematical derivation on the subject. Since both forms exist widely in the literature, here we opt to apply both approaches to CSA to demonstrate discrepancies between the approaches. 2.2 A new approach to estimate the time constant from a noise-contaminated power spectrum In the PSA I07 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 the noise part is fitted with a line of an unconstrained slope (blue dashed line) and then the transfer function is fitted (black solid line). This approach can be problematic as we will show in Sec. 4.1 below, because if the resulting line has a slope less than unity, its extrapolation 30 to lower frequencies erroneously removes true signal. By contrast, a line with a predefined slope of 1 can only be fit for very noisy spectra where the 1:1 is fully reached. Here we introduce a new alternative approach (hereafter PSA A20 ), which performs the estimation of the time constant and accounts for instrument noise in a single non-linear comprehensive fitting step.
Equation (4) can be extended to the following equation, which includes also the noise component of the scalar: 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 5 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 where b is the y-axis intercept of the power spectra of the white noise multiplied with f , which is shown linearly in log-log scale. All terms in the equation have been multiplied by f because this is the standard normalisation used to depict the spectral 10 density functions (Fig. 1). The detailed derivation of Eq. (7) can be found in Appendix A.
In this approach, the frequency ranges used for fitting both for transfer function and noise are not separated, hence the fitting is performed across the entire frequency domain.

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, two-dimensional coordinate rotation of the wind velocity vector, and linear de-trending to D 1 . We did not apply For PSA, we calculated power spectra of T , following Sabbatini et al. (2018), normalised 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. 20 Then we took an ensemble average of the 70 power spectra. For PSA I07 , we first removed the noise from the power spectra (see Sect. 3.2.3), then retrieve the time constant (i.e. τ P SA I07 ) 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 PSA A20 , we obtained the time constant (i.e. τ P SA A20 ) via fitting Eq. (4) using the entire frequency domain.
Low-pass filtering introduces a time-lag in the scalar of interest in addition to any physical time-lag that may be caused by 25 transport through sampling lines and/or sensor separation (Massman, 2000;Ibrom et al., 2007b). In our case, since the scalar of interest is represented by the attenuated T , this time-lag occurs between T and w, two entities that are not separated by a further physical time-lag. In order to demonstrate the effect of low-pass induced time-lag on the time constant estimation in CSA, we generated a dataset in which the time-lag was adjusted via the maximisation of the crosscovariance, which was explored in addition to the original not time-lag adjusted dataset. Next, we calculated cospectra of T and w of both datasets, (3). All fits were estimated for the frequency range from 0.01 to 2 Hz. For a more complete analysis on the effect of low-pass induced time-lag on time constant estimation see the companion paper (Peltola et al., 2020).
In summary, for the two methods of PSA and three methods of 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

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 Lee, 2002): where χ O is the output of a scalar sensor, χ I (t) is the true scalar concentration (i.e., input) and τ LP F is the characteristic time constant of the sensor response (Massman and Lee, 2002). The spectral response (h LP F (ω)) of such system can be obtained by the Fourier transform of the ratio of output signal to the input signal, i.e. χ O /χ I (Horst, 1997): where ω = 2πf and j = √ −1.
The desired (i.e. low-pass filtered) output data in frequency domain can be estimated as where Z I 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 Z O and taking the real part. In practice, the complex conjugate of h LP F is used to derive the correct 5 temporal lag (scalar lag with respect to w).
We followed this procedure to filter T with five different τ

Noise superimposition
In this study we use Gaussian white noise, which has equally distributed spectral densities across all frequencies (Stull, 2012), 10 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.

Noise removal
The instrumental noise is often removed via the approach introduced by (Ibrom et al., 2007a), in which, to detect the noise, a line with unconstrained slope is fitted to the power spectra in logarithmic frequency space at a certain frequency range, which is assumed to be solely dominated by noise. Then, the fitted line is extrapolated towards lower frequencies and finally subtracted from the original spectra in logarithmic frequency space as described in Ibrom et al. (2007a) (cf. dashed lines in Fig. 4). The 20 boundaries of the frequency domain of fitting are defined by visual inspection of the power spectra. In our study, the optimal frequency domains used for noise fitting were detected as 3, 2.6 and 2.3 to 5 Hz for the different attenuation level 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.

Data processing when estimating long-term budgets
Data processing steps when estimating the long-term budgets using the D 2 dataset are summarized in Fig. 3. 25 We applied regular EC data processing, which included de-spiking, coordinate rotation, de-trending. Next, the T time series were deteriorated with values of τ that varied between 0.1 and 0.5 s to mimic a realistic range of scalar attenuations (mimicking the conditions, e.g., of CH 4 , N 2 O). Later, the low-pass induced time-lag was accounted for via maximisation of the cross-covariance, which was followed by the calculation of the covariances. Then we applied the frequency response correction using Eq. (1), where F corr was estimated using the method proposed by Fratini et al. (2012) 1 : 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 D) for CO instead of parameterisation of F corr proposed by Ibrom et al. 5 (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 ms −1 ), unrealistic sensible heat fluxes and non-stationary conditions (Foken and Wichura, 1996), eliminating 553 half-hourly data points out of 1440 (ca. 38%).
The time constant (τ ) in H emp is either estimated by PSA (i.e., τ P SA I07 and τ P SA A20 ) or CSA (i.e., τ CSA √ H,sync , τ CSA √ H and τ CSA H ), and this yielded cumulative fluxes F P SA I07 , F P SA A20 , F CSA √ H,sync , F CSA √ H and F CSA H , respectively. For sim- increases and therefore, goodness of fit increases (Table 1). As discussed in Sect. 2.2, white noise causes a slope of one, and thus any discrepancies from one results in overestimated noise removal. Based on Fig. 4 alone, it is evident that removing 1 Here we preferred using the square-root to describe the true transfer function as it is a good approximation when maximisation of the cross-covariance is used for the time-lag correction as shown by Peltola et al. (2020).
https://doi.org/10.5194/amt-2020-478 Preprint. Discussion started: 23 December 2020 c Author(s) 2020. CC BY 4.0 License. the noise from power spectra can be done with higher accuracy when the high frequency attenuation increases and/or SNR decreases.
The results of the time constant estimation are shown in Fig. 5a, and 5b for PSA I07 and PSA A20 , respectively. Results are presented as medians with 25 th and 75 th percentile ranges for the repeat simulations. It can be seen that for PSA I07 the interquartile range (IQR) expands as the amount of noise increases, while for PSA A20 it is small and almost constant with only 5 a slight increase. PSA I07 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 via 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 one, and thus better fitting parameters. However, especially for low attenuation 5 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 shortcoming of the linear fitting, the visual inspection of both detecting the frequency ranges for noise removal and the H fitting constitutes another uncertainty source due to its subjectivity. It requires expertise on the topic, and is hard to Figure 5. Time constants calculated using the power spectra approach, comparing (a) P SAI07, and (b) P SAA20, in several low-pass filtering condition (τLP F =0.1 -0.5 s) as a function of signal-to-noise ratio (SNR) over the range 10-1.1, which correspond to the amount of noise (e.g. 10-90 %). The solid red curve represents the median, while the shaded area represent inter-quartile ranges. The dashed line corresponds to the expected value (τLP F ), which was used for the artificial low-pass filtering.
automise 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.
PSA A20 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 5 the PSA A20 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 PSA A20 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 characterises the noise. This said, in many EC studies, the type of the 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 characterise 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 5 measurements with very low SNR and those of white and blue noise (see Appendix C). It should be noted, however, that there are situations when noise is complex and difficult to predict a priori. Figure 6 illustrates cospectra of three low-pass filtered cases (i.e. τ =0.1, 0.3 and 0.5 s) with most noisy conditions (i.e. lowest SNR of 1.1) in addition to raw (i.e. original) cospectra. It shows that noise contamination does not affect the shape of the 10 cospectra.

Time constant estimation with CSA
CSA H systematically underestimates the time constant, while, in contrast, CSA √ H overestimates, the accuracy of which is further significantly improved when the time-lag correction is applied to the time-series, in CSA √ H,sync (Fig. 7). The reason for this improvement is discussed in detail in the companion paper (Peltola et al., 2020).
The accuracy of CSA H decreases with increasing attenuation. The spectral energy of two variables comprises a real part, 15 i.e. the cospectrum (Co), and an imaginary part, i.e. the quadrature spectrum (Q) (Stull, 2012). The latter is often neglected by assuming the attenuation does not change the phase between the two variables (Horst, 2000). In order to test this assumption we further investigated the effect of Q on the time constant estimation (see Appendix B). In short, application to a couple of case studies shows that CSA H can under-or over-estimate the time constant due to neglected Q, depending on the sign of Q.
The bias is larger for higher attenuation since the central term (2τ c Q/Co term in Eq. (B1)) increases linearly with τ . We 20 showed that (Table B1 and B2), the bias can be eliminated when the ratio of Q over Co, both of which are obtained from unattenuated records, is utilized in the H (Eq. (B1)). However, the use of the ratio is not applicable in real-world data of closed-path systems where Q and Co are distorted in case of attenuated time series, hence this uncertainty can be noted as a shortcoming of the CSA.
In the case of CSA √ H,sync , the variation around the expected values can be attributed to the shortcomings in the maximi-25 sation of the crosscovariance 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.
In practice, the time-lag correction when utilized with the square-root of H can eliminate the above mentioned phase effect, the theoretical premise of which is further investigated in Peltola et al. (2020).
Lastly, as expected, the level of noise does not affect the accuracy of any of the CSA methods at any SNR level as the random 30 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 P SA A20 approach. P SA I07 based fluxes showed a bias of ±2% with negligible response to the SNR level, whereas P SA A20 showed a somewhat similar behaviour with slightly better accuracy reflecting the more accurate time constant estimation previously shown.

5
Fluxes based on CSA √ H,sync were very close to the expected value with similar biases as the PSA methods. CSA √ H overestimates the fluxes up to 3%, while the CSA H method underestimates fluxes by up to 4%. These findings are partly consistent with the observed biases on time constant estimation, meaning that where the time constant and the low-pass filtering were overestimated (e.g. with the P SA I07 especially for τ =0.1 s and CSA √ H ), the spectral correction factor and thus the fluxes were overestimated, too. In summary, the findings indicate that using PSA methods (particularly P SA A20 ) for time constant estimation provides reasonably accurate FRC. On the other hand, for CSA methods, CSA √ H,sync better approximates the time constants, hence FRC. These results are in agreement with the analysis presented in Peltola et al. (2020), showing that CSA √ H,sync well approximates the effect of phase shift on the estimation of the time constant and flux correction factor. See more details in Peltola et al. (2020).

4.4
Review of typical signal-to-noise ratios and response times encountered during closed-path flux measurements 15 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, configuration of instruments, and setup of EC system. Ibrom et al. (2007a) examined Figure 7. Time constants calculated using the cospectral approaches, i.e. CSAH (red), CSA √ H (green), and CSA √ H,sync (blue), for several low-pass filtering condition (τLP F =0.1 -0.5 s) as a function of signal-to-noise ratio (SNR) over the range 10-1.1, which correspond to the amount of noise (e.g. 10-90 %). The solid curves represent the medians, while the shaded areas represent iner-quartile ranges. The dashed line corresponds to the expected value (τLP F ), which was used for the artificial low-pass filtering. fluxes of water vapour and CO 2 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 f c of CO 2 as low as 0.325 Hz (i.e. ca. τ = 0.5 s), indicating strong attenuation, while H 2 O 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) Inc., Billerica, USA) that also measures mole fractions of CO 2 and H 2 O. They reported high attenuation for their measurement setup (time constant of 0.68 s for their EC system), and high noise disturbance in their power spectra (their Fig. 7) with a SNR of about one. 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 CO 2 , assuming that both fluxes are affected by the same attenuation.

5
In addition, N 2 O fluxes measured over urban areas (Järvi et al., 2014), and CO 2 fluxes over lake ecosystems  can also be examples of low SNR conditions.
Given the wide range of SNR and attenuation conditions summarised above, we analysed only a limited range of SNR and attenuation. Also, the impact of a 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 dif-10 ferent 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.
Finally, the key constraint of the study was that we artificially simulate the various attenuation and SNR conditions, which should be verified with in situ measurements.
15 Figure 8. Relative biases of the cumulative fluxes derived with the various approaches compared with the reference flux as a function of SNR for different attenuation time scales (0.1-0.5 s), for P SAI07 (blue) and P SAA20 (black), calculated as, e.g., 100(FP SA I07 −FREF )/FREF .

Conclusions
Here we investigated the limitations of two commonly used approaches to empirically estimated 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 setups 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. PSA I07 and PSA A20 .
The latter is newly introduced here and does not require noise removal prior to fitting of the response function. For CSA, we examined three approaches: (1) utilising the transfer function (H) suggested by Horst (1997), i.e. CSA H , (2) using √ H, i.e.

5
CSA √ H , and (3) implementing √ H with shifted w time series via maximisation of the crosscovariance, i.e. CSA √ H,sync . We generated an artificially 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. PSA I07 overestimated the noise contribution and consequently the signal loss and time constant for low attenuation conditions, but better performance was found 10 as attenuation increases. The new PSA A20 approach successfully estimated the time constants regardless of the attenuation and SNR level, identifying the noise and signal comprehensively, 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. By contrast, systematic differences were found between the results from three approaches based on cospectra, where the choice of 15 the shape of H differs in the literature. CSA H underestimated the time constants with increasing bias as attenuation increases, the accuracy of which was shown to be affected by neglecting the filtering effects on the quadrature spectrum, which cannot be avoided. CSA √ H overestimated the time constants, the bias of which was mostly eliminated when time-lag correction was ap-plied in CSA √ H,sync . See the companion paper (Peltola et al., 2020) for a thorough discussion of the interdependence between time-lag estimation via cross-covariance maximisation and frequency response corrections.
We then examined the effect of the different approaches to estimate the time constants on the cumulative fluxes: fluxes corrected using the P SA I07 based time constants showed the bias of ±2% in comparison with reference fluxes (see Sect. 3.3), where P SA A20 showed similar behaviour, yet slightly better accuracy. By contrast, fluxes corrected using CSA H based 5 time constants were underestimated by up to 4%, while CSA √ H overestimated the fluxes by up to %3, the bias of which was further reduced with CSA √ H,sync , varying within ±2%. The SNR did not affect the accuracy of either PSA or CSA approaches, alleviating concerns on EC flux measurements with low SNR level.
In summary, for the empirical estimation of parameters of H of closed-path EC systems, our findings showed that PSA A20 is the most accurate, precise and robust method when power spectra are used. This finding was independent of SNR and degree 10 of attenuation. On the other hand, using the square root of H, which is a good approximation of the attenuation of cospectra in the presence of phase shifts (Peltola et al., 2020), provided the correct estimation of τ and F corr when cospectra are used, together with time-lag quantification by cross-covariance maximisation.
Finally, given the constraints of this study, we encourage additional studies based on the real attenuation and SNR conditions, investigating also other type of noise contamination to provide a step forward in efforts to standardise the EC method, which 15 is of great importance to avoid systematic biases of fluxes and improve comparability between different datasets.
Appendix A: Derivation of the PSA A20 In this section, we derive the new approach for the PSA (i.e. PSA A20 ) 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 20 holds (i.e. power spectra of all scalars follow a similar shape). After these assumptions the power spectrum (S χ ) of χ can be written as: where σ 2 χ and σ 2 T are the variances of χ and T related to turbulent signal (no attenuation or noise), S T (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 25 that here the noise in the T measurements was neglected as being small except in very low turbulence conditions (when eddycovariance 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: where the proportionality constant F n was introduced due to the effect that noise adds to σ 2 χ and thus biases the normalisation (see Ibrom et al., 2007a). The spectral density is best shown in log-log space as wide range of frequencies and spectral densities are well displayed (Stull, 2012), with natural frequency (f ) on the x-axis, and spectral density multiplied with f on the y-axis.
After multiplication by f , Eq. (A2) becomes: = log(f ) + log(b), and in non-logarithmic space as e log(f )+log (b) , which is f b. Substitution into Eq. (A3) yields: It is worth mentioning that this method can be used to retrieve the variance of noise as well, since b equals the ensemble 10 averaged noise power spectra (S χ,n ) divided with σ 2 χ . However, we do not further examine it as it is not the main interest of this study.
Appendix B: Assessment of the effect of quadrature spectrum on time constant estimation in the CSA In Sect. 4.2, we briefly presented the effect of Q on time constant estimation in the CSA H . In this section, the details are provided. 15 The transfer function used in this study (i.e. Eq. (3)) is the simplified version of the transfer function of Horst (2000) (their Eq. (5)) which neglects Q and assumes that the sonic anemometer has a perfect frequency response. The theoretical transfer function including Q can be written as: where Q is quadrature spectrum and Co is the cospectrum of two variables and τ is the system time constant. 20 For this particular analysis we used two T datasets, the first of which is D 1 from Siikaneva and second of which is from Hyytiälä (see Peltola et al., 2020). We filtered the data by following the same procedure described in Sect. 3.2.1 with various τ LP F ranging from 0.05 to 1 s. Then, we used different transfer functions (i.e. Eq. (3) and Eq. (B1)) to obtain time constants to examine the effect of Q.    W m −2 , and the wind speed was between 2.1 and 3.7 m s −1 .

Appendix C: Identification of instrumental noise
Here we perform a short analysis to characterise the type of noise in an example high frequency CH 4 time series by comparing the spectra of 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 Relationships), Power spectra of methane were calculated using measurements from the 10 July at 12:00-14:30. The high frequency range of the normalised frequency weighted CH 4 power spectrum follows the f +1 power law scaling, which is consistent with white noise and inconsistent with, e.g., blue noise (Fig. C1). f Co(f ) w χ = 2 π n/nm 1 + (n/nm) 2 , where n is the normalised frequency, nm is the cospectral peak frequency, which is derived from in situ measurements in Siikaneva as: where z is the measurement height, d is the displacement height, and L is the Obkuhov length, a measure of atmospheric stability.

5
Author contributions. IM and OP designed the study and TA processed and analyzed the data. UR investigated the distortion caused by quadrature spectrum in CSA. UR, TA and AI investigated the limitations of noise removal procedure of PSAI07. AI developed the PSAA20 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.