The high frequency response correction of eddy covariance ﬂuxes. Part

. The eddy covariance (EC) technique has emerged as the prevailing method to observe ecosystem - atmosphere exchange of gases, heat and momentum. EC measurements require rigorous data processing to derive the ﬂuxes that can be used to analyse exchange processes at the ecosystem - atmosphere interface. Here we show that two common post-processing steps (time-lag estimation via cross-covariance maximisation, and correction for limited frequency response of the EC measurement system) are interrelated and this should be accounted for when processing EC gas ﬂux data. These ﬁndings are applicable 5 to EC systems employing closed-or enclosed-path gas analysers which can be approximated to be linear ﬁrst-order sensors. These EC measurement systems act as a low-pass ﬁlters 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, ﬁlters). Time-lag estimation via cross-covariance 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 χ , 10 which distorts the measured cospectra between w and χ and hence has an effect on the correction for dampening of EC ﬂux 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 ﬂux data. Based on analyses using EC data from two contrasting measurement sites, we show that the low-pass ﬁltering induced time-lag increases approximately linearly with the time constant of the low-pass ﬁlter, and hence the importance of H p in describing the high frequency ﬂux loss increases as well. Incomplete description of

(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 carbon budget of different type of ecosystems (Baldocchi, 2003;Vesala et al., 2008;Mammarella et al., 2015).
Gas flux measurements are made using a three-dimensional sonic anemometer and a gas analyser, which are able to provide 5 fast-response measurements of turbulent fluctuations of vertical wind velocity and gas concentration  :::::::::::::::::: . All EC systems, including open-path, enclosed-path or closed-path gas analysers, act as low pass 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 small-scale variation. The main reasons for co-spectral attenuation are inadequate frequency response of the sensor, sensor separation and line averaging, as well as, in closed-path systems, the sampling of air trough :::::: through : tubes and 10 filters.
The frequency response correction is usually performed based on a priori knowledge of the system transfer function and the unattenuated cospectrum: (1) Here CF is the estimated spectral correction factor, Co w,χ the normalised unattenuated cospectrum between scalar (χ) and 15 vertical wind speed (w), f the frequency, H lpf the total transfer function describing the low-pass filtering, and f 1 and f 2 are 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 low-frequency :::::::: high-pass transfer function (Rannik and Vesala, 1999) is not included in Eq. (1), as here we focus only on the low-pass filtering effect of the EC system. The high-frequency ::::::: low-pass ::::::: filtering : transfer function H lpf can be derived 20 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 co-spectral 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 25 open-path 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, where the model cospectra and H lpf are estimated using in situ measurements. This is the preferable approach for closed-path systems, where 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 absortion/desorption processes inside the sampling tube (Ibrom et al., 2007a;Nordbo et al., 2014;Runkle et al., 2012) may contributed :::::::: contribute to the EC system cut-off 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 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., 2009Mammarella et al., , 2016 for the empirical 5 retrieval of H lpf . The companion paper (Aslan et al., 2020) has investigated methodological issues and flux uncertainty related to these two methods under different attenuation conditions and signal-to-noise ratio scenarios.
In EC systems, particularly those using closed-path gas analysers, the measurements of vertical wind velocity and gas concentration are not co-located, and a time delay between the two signals exists  ::::::::::::::::::  . The standard procedure is to determine this time delay for each averaging period (typically 30 minutes) by maximising the 10 related cross-covariance 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 IRGA sampling cell, separation distance between the inlet and the centre of the sonic anemometer path), but it can additionally reflect the low-pass filtering of the measured signal (Massman, 2000;Ibrom et al., 2007b, a).
In this study, we investigate how the low-pass filtering induced phase shift affects the estimation of the high frequency flux 15 loss, and we show the implications that occur when H lpf time constants are empirically derived from measured power spectra or cospectra, and related CFs values are applied to covariances estimated from the cross-covariance 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 time scale characteristics. 20 2 Theory

Transfer functions
The measured cross-spectrum (Cr m ) between the fluctuations of vertical wind speed (w ) and the attenuated and lagged scalar (χ ; in this paper attenuated variables are denoted withˆ) can be described by (e.g. Massman (2000)): 25 where Z w and Z χ are Fourier transforms of the time series of :::::: vertical ::::: wind ::::: speed : (w and : ) :::: 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 = √ −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 cross-spectra measured with a typical EC system employing closed-path gas analyser 30 where 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 high-pass filtering on flux attenuation.
Following, e.g., Horst (1997) and Massman (2000), the scalar sensor is approximated to be a linear first-order sensor for which h χ can be written as where τ is the time constant (s) of the filter, also called response time. Note that Z w Z * χ is the true cross-spectrum (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 w Z * χ = Co + jQ). Recall also that the integral of Co equals the unattenuated covariance w χ , i.e. the vertical turbulent flux of scalar χ. Following previous studies (Horst, 1997;Massman, 2000), the quadrature spectrum (Q) can be 10 assumed to be zero for stationary turbulent flow and hence Cr=Z w Z * χ ≈ 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 15 the travel time of scalar χ signal in the sampling system is zero. This means that w and the attenuatedχ (ˆdenotes attenuation throughout the text), depending on frequency, can be out of phase due to low-pass 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 first-order system (H) as 20 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 cross-spectrum of the attenuated scalar with itself yields: where Z χ Z * χ 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.
The centre of the sonic anemometer was displaced 20 cm horizontally and 1 cm vertically from the intake of the gas analyzer.
The measurement setup was designed to closely follow the ICOS EC measurement protocols (Franz et al., 2018;Rebmann 25 et al., 2018). The measurements were conducted between May and August 2019.
The Siikaneva fen site is located in Southern Finland (61 o 49.9610 N, 24 o 11.5670 E; 160 m a.s.l), consisting mainly of sedges (Eriophorum vaginatum, Carex rostrata, C. limosa) and Sphagnum-species :::::::: Sphagnum ::::::: species, namely S. balticum, S. majus and S. papillosum with the : 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 30 used in this study were measured with 10 Hz sampling frequency using a 3-D sonic anemometer (Metek USA-1, GmbH, Elmshorn, Germany) and a closed-path analyzer (LI-7000, LI-COR Biosciences, :::::: Lincoln, :::: NE, : USA). The sonic anemometer and the gas inlet were situated at 2.75 m above the peat surface, and the air was drawn to the analyzer through a 16.8 m long heated inlet sampling line. The center of the sonic anemometer was displaced 25 cm vertically above the intake of gas analyzer.
A short dataset (hereafter D S1 ) measured between 9:30 and 11:30 on 16-Jun 2013 at the Siikaneva site was used in the method validation (see Sect. 4.2), whilst long-term datasets from both Siikaneva (hereafter D S2 ) 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.

Data processing
High frequency 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., 10 2016), 2) coordinates were rotated using double-rotation 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 point-by-point to be relative to dry air, if not already done internally by the gas analyser (LI-7200), 4) turbulent fluctuations were extracted from the measurements using block-averaging and 5) time lags between gas or filtered T (see Sect. 3.2.2) and vertical wind time series were accounted for using cross-covariance maximisation. Power spectra and cospectra were 15 calculated after these processing steps. Spectral corrections and the related correction factors (CF ) 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 non-stationary fluxes (Foken and Wichura, 1996). Site specific friction velocity thresholds (0.17 and 0.3 m/s for Siikaneva and Hyytiälä, respectively) were used to discard gas flux data during low turbulence periods. Furthermore, at Hyytiälä wind directions between 150 • and 230 • were discarded since here the mast 20 construction disturbed the air flow. At Siikaneva an omnidirectional sonic anemometer was used mounted to 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): 25 where u * is the friction velocity, v k = 0.4 the von Karman constant, g the acceleration due to gravity and θ v is virtual potential temperature. Overline denotes temporal averaging and primes deviations from the mean. From the Obukhov length, the stability parameter was calculated as ζ = z−d L , where z is the measurement height and d is the zero plane displacement height.
3.2.1 Empirical estimation of the high-frequency attenuation ::::::: low-pass :::::::: filtering :: of ::: the ::::: signal The total transfer function of an EC setup can be estimated empirically from measured w-χ and w-T cospectra: or, alternatively, from the power spectra of the measured scalar and T : where Co w,χ denotes the cospectral density between w andχ, and w χ represents the covariance between w andχ. Sχ ,χ and σ 2 χ denote the power spectral density and variance ofχ, 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 Hz 5 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 high frequency losses (Ibrom et al., 2007a). Here, w and T were approximated to be free from any low-pass 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 ensemble-averaged spectra from measurement periods that fulfilled the following criteria: flux stationarity (Foken and Wichura, 1996) was below 0.3, 10 u * > 0.1m s −1 , w T > 0.02K m s −1 . For CO 2 it was also required that its turbulent flux was below -0.05 µmol mol −1 m s −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 (similarly 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 15 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 high frequency attenuation in our companion paper (Aslan et al., 2020).
Often both H e,cs and H e,ps are approximated by H with a value for τ that describes the high-frequency attenuation ::::::: low-pass :::::: filtering : of the EC setup according to Eq. (5). Then the measured flux is corrected for high-frequency attenuation ::::::: low-pass :::::: filtering : by multiplying it with CF (see Eq. (1)). There are different ways to estimate Co w,χ for Eq. (1). Here we limited the 20 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).
See the definitions for H and Hp in Eqs. (5) and (6), respectively. Values for τ (and t lpf in method :::::: Method 4) ) were estimated with nonlinear least squares fit to He,ps or He,cs where Fn and τ (and t lpf in method :::::: Method 4) ) were used as fitting parameters. The additional normalisation factor Fn took into account any inaccuracies in estimation of variances and covariances in calculation of He,ps or He,cs (Eqs. (13) and (12)), respectively. The fits were weighted with temperature cospectra, in order to give more weight to frequencies where the scalar fluxes were high.

Low-pass filtering of temperature data
In order to evaluate the performance of the different spectral correction methods presented in Sect. 3.2.1, high frequency T 5 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. (2020), 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 s, 0.3 s and 0.7 s. The correct value for CF was obtained from the ratio between w T and w T , whereT was the time series attenuated with the methodology described above. Note that here w T was calculated using cross-covariance maximisation in order to 10 mimic regular EC gas flux processing, meaning that the lag caused by low-pass 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. numerically for the constant C as a function of frequency f . Fig. 1 (upper plot) illustrates that such solution is dependent on frequency and varies over a relatively large interval up to the frequency 2 3 π τ (the relationship holds well up to such 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 one. The ::: For :::::::::: illustration, ::: the lower plot of Fig. 1 presents the left and right side terms of the approximate equality using the constant 5 coefficient C = 0.63 :::::::::::: (H p ≈ 1/ √ H) ::::::::: calculated ::::: using : a ::::: value :::: 0.63 ::: for ::: the ::::::::: coefficient :: C. In spite of the large variation of the exact coefficient with frequency ::::: shown :: in ::: the ::::: upper :::: plot :: in ::: Fig. :: 1, the selected constant value produces very good correspondence for the whole frequency range . A :::::::: suggesting :::: that :: C ::: can ::: be ::::::::::: approximated ::: to :: be ::::::: constant ::::: albeit :::::: 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 τ = 0.01 π τ to 2 3 π τ for a range of time constants from 0.05 to 1 s. We observed that the 10 optimum coefficient was independent of τ and equal to 0.63 (with 2-digit 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 ).

15
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 of 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 and hence it can be used :::::::: supporting ::: the ::::::::: estimation :: of :::: t lpf as a fitting parameter when evaluating the high frequency response of the EC system :: in ::: the :::::: Method :: 4 ::::: when ::::: fitting ::::: HH p :: to ::::: H e,cs (see Sect. 3.2.1 , Method 4 ::: and ::::: Table : 1).

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 cross-covariance functions from the Siikaneva site with three different values of τ applied to the same 2-hour time series. The peak of the w-T cross-covariance shifted to longer positive times (i.e.T lags w) and the cross-covariance maximum decreased as τ increased. This is related to the phase shift and flux attenuation caused by 25 the low-pass filter, respectively. The low-pass 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 where the time lag derived from cross-covariance 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 low-pass 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  (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 :: In ::::::::: particular, :::: HH p :::::: shows ::::::: negative :::::: values :: at :::: high ::::::::: frequencies ::::::: whereas : √ H . :::: does ::: not. : When the correction factor was calculated with Method 4, which takes into account the low-pass filtering induced phase shift, the attenuation factors (1/CF ) agree with the normalised cross-covariance maxima (see Fig. 3). This indicates that 5 the attenuation of cross-covariance 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 cross-covariance values at zero lag as per predicted by theory (Sect. 2). Note that the time lag discussed here only represents the low-pass filter effect, not other effects that also alter the time lag. This is only possible because we work with degraded temperature time series, where phase shifts through low-pass filtering are the only cause for time delays.
In order to evaluate the four methods to estimate CF further, turbulent T data were low-pass filtered with three values of 10 τ (0.1 s, 0.3 s 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 :::::: albeit ::::: slight ::::::::: differences ::: can :: be :::::: noted. The performance of all four methods decreases with increasing τ at both sites, and Method 1 showed 15 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 Relative difference between the fluxes based on the different correction methods applied to :::: mean : sensible heat fluxes attenuated :::::: obtained : with ::: the different values for τ :::::::: correction :::::: 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 τ /ITS (Fig. 6) where ITS is the integral time scale of w T , calculated from the autocovariance of w T time series by integrating the normalised autocovariance function until its first zero crossing (e.g. Kaimal
Hyytiälä, respectively. ::: For ::::::: Hyytiälä ::: the ::::: signal ::::: travel :::: time ::::::::: calculated :::: from ::: the :::::::: sampling :::: tube ::::::::: dimensions :::: and :::: flow ::: rate :::: was :::: 0.06 : s. : 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 time than the other two methods that were projected to agree (see e.g. 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 5 with small attenuation both τ and t lpf are likely difficult to estimate accurately by fitting HH p to H e,cs . The low-pass filtering induced lag (t lpf ) can attain values only with the temporal resolution of the underlying data itself ::: 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., 2020). 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 gas analyser induced low-pass 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   (Nordbo et al., 2013). Hence different   The time lag 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 cross-covariance maximisation showed linear dependence on estimated response times in relative humidity bins with a similar dependency at both sites (y=0.47x-0.05 s :::::::::: y = ax + b, ::::: where ::::::: a = 0.47 :::: and :::::::::: b = −0.05 s for Siikaneva and y=0.44x-0.07 s ::::::: a = 0.44 ::: and :::::::::: b = −0.07 s : for Hyytiälä, where y is the time lag difference and x response time. Fits were based on orthogonal linear regression). Difference between H 2 O and CO 2 time 15 lags resulted from stronger low-pass 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 low-pass filtering effects, were the same for H 2 O and CO 2 .
All values are in %. The CF values estimated with the four methods agreed at Hyytiälä for low RH periods, but they diverged at high RH periods.

Site
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 :::::: Method 5 1 ) overestimated and method ::::::::::: overestimated ::: and ::::::: Method : 3 ) underestimated CF similarly as here for H 2 O at high RH (and hence large τ ). In contrast, for Siikaneva Method 1 gave systematically higher value 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.

10
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 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).

Correct form for cospectral transfer function
There has been a long-standing debate about what is the correct form for cospectral transfer function when the scalar measurements are done with first-order 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 √ H (see Eq. (5)), which was later deemed erroneous, e.g. by Eugster and Senn (1995); Horst (1997Horst ( , 2000 and others. Partly based on similar deriva-5 tions as shown here in Sect. 2, the correct form for the cospectral transfer function was argued by Horst (1997Horst ( , 2000 and Massman (2000) 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 transfer function 10 for the amplitude spectrum (Eq. (10)).
Briefly summarising the findings in this study, it was shown in Sect. 2 in accordance with e.g. 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, cross-covariance maximisation inadvertently accounts also for the scalar lowpass filtering related time shift (t lpf ) and consequently induces a phase shift between the lag corrected attenuated scalar and 15 vertical wind time series. Hence in the case that cross-covariance maximisation is used to align the two time series, the correct form for cospectral transfer function is HH p , where 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 √ H. Therefore the debate about which is the correct form for cospectral transfer function (H or √ H) relates to the low-pass filtering induced phase shift as discussed already by Eugster and Senn (1995) and Horst (2000), albeit with a different conclusion. :::: Note :::::::: however, :::: that :::: √ H :: is :: an ::::::::::::: approximation, ::: not ::: an ::::: exact ::::::::::: representation ::: of 20 :: the ::::::::: cospectral ::::::: transfer ::::::: function. :

Summary and conclusions
The influence of low-pass filtering induced phase shift on estimation of high frequency response of an EC setup was analysed.
The analysis assumed that the EC setup consisted of a fast-response anemometer and a linear, first-order-response scalar sensor.
Spectral corrections aiming at correcting the EC fluxes for high frequency attenuation ::::::: low-pass ::::::: filtering were estimated with 25 four methods: three widely used methods and one newly proposed method which specifically accounts for the interaction between the low-pass filtering induced phase shift and high frequency attenuation. Based on theoretical considerations and experimental results we come to the following conclusions: -Cross-covariance maximisation overestimates the signal travel time in the EC sampling line since it inadvertently accounts also for the lag caused by low-pass filtering of scalar time series caused by non-ideal measurement system. The 30 bias in the estimated time lag depends linearly on low-pass filter response time (τ ) with a small additional dependence on turbulence time scale. :::::: Hence, ::::::::::: investigation :: of ::::: other ::::: means ::: for ::::::::: estimating ::: the ::::: signal ::::: travel :::: time ::::: might ::: be ::::::::: warranted.
-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 first-order response time due to the use of the incorrect cospectral transfer function, it resulted only in a small bias for the flux correction factor, 15 since the method describes the attenuation of cospectra accurately once the cross-covariance 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 bias in EC fluxes can be 10% or above if method 1) is used for short tower measurements with moderate attenuation.

20
-All flux calculation algorithms which rely on cross-covariance maximisation and at the same time use H to describe the attenuation of the cospectra 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). Hence, for instance the analytical formulas of Horst (1997) and Massman (2000) will result in biased fluxes when an accurate estimate of scalar sensor response time is used in these algorithms.

25
-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. Similarly, Wintjen et al. (2020) showed that the transfer function related to the phase shift (H p ) needed to be taken into account when processing reactive nitrogen fluxes.
Code and data availability. Data and matlab codes to reproduce Figures 3, 4 and 5 will be uploaded to open data repository Zenodo upon acceptance of this manuscript along with flux time series processed with the four methods (Table 1) at the two sites.

5
Author contributions. 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 the first draft and made improvements.
Competing interests. The authors declare that they have no conflict of interest.