Optimal Frequency-Response Corrections for Eddy Covariance Flux Measurements Using the Wiener Deconvolution Method

We describe a new direct correction approach to accurately restore frequency attenuated eddy covariance (EC) measurements. The new approach utilizes the Wiener deconvolution method to optimally estimate the original signal from noisy atmospheric measurements. Key features over conventional EC spectral correction methods include (i) the use of physics-based response functions, (ii) the ability to account for the non-linear phase contributions, and (iii) the direct restoration of the original signal rather than simulating the effect on an ideal reference spectrum. The new correction approach is compared to conventional spectral correction methods in a numerical simulation where the magnitude of the key limitations of conventional methods is explored under conditions relevant to common EC set-ups. The simulation results showed that the spectral correction methods commonly used for calculating EC fluxes introduced systematic error up to 10% to the restored fluxes and substantially increased their random uncertainty. The errors are attributed to the effect of using inappropriate response functions, failing to account for the contribution of the non-linear phase, and due to the assumption of spectral similarity on the scale of averaging intervals. The Wiener deconvolution method is versatile, can be applied under non-ideal conditions, and provides an opportunity to unify analytical and “in-situ” spectral correction methods by applying existing transfer functions to directly restore attenuated spectra. Furthermore, the Wiener deconvolution approach is adaptable for use with various micrometeorological measurement techniques such as eddy accumulation and flux profile measurements.


Introduction
The eddy covariance (EC) method has become the de facto method for measuring the atmospheric exchange of energy, momentum, and atmospheric trace constituents. EC enables B Anas Emad anas.emad@uni-goettingen.de 1 Bioclimatology, University of Göttingen, Büsgenweg 2, 37077 Göttingen, Germany 30 A. Emad direct and continuous measurements of a range of atmospheric fluxes, such as CO 2 , H 2 O, sensible and latent heat fluxes, on a scale suitable for ecosystem studies (Baldocchi 2014;Hicks and Baldocchi 2020). The turbulent flux of a scalar in EC is calculated as the covariance between the scalar concentration c and the coincident vertical wind velocity w. In ecological applications, the covariance between w and c is commonly calculated for a long-enough averaging interval (30 to 60 min) using the high-frequency measurements of the scalar concentration and the wind. Typically, measurement frequencies of 10 Hz or more are required to resolve relevant scales of flux-carrying turbulent eddies (Aubinet et al. 1999).
Fluxes measured using the EC method suffer from unavoidable biases caused by the frequency response limitations of the EC system and data processing methods. Measured signals required for flux calculation are often attenuated due to interactions with the measurement apparatus. This attenuation can result from limitations of the instrument design, data processing methods, or due to physical processes such as water adsorption and desorption on the walls of tubes. Signal attenuation is notably more significant in closed-path EC systems and when sorbing scalars like water vapor are involved, with the underestimation of measured fluxes under these configurations reaching up to 40 (Shaw et al. 1998;Ibrom et al. 2007).
Attenuated signals can be more comfortably identified and treated in the frequency domain by inspecting their spectra and cospectra. Therefore, depending on the part of the spectrum affected, the attenuation can be classified into high-frequency losses equivalent to the effect of a low-pass filter and low-frequency losses equivalent to the effect of a high-pass filter (Massman and Clement 2005;Foken et al. 2012).
Methods for correcting spectral attenuation can be broadly put into two classes based on their assumptions and how they are applied. The analytical method, also called the transfer function method or the theoretical method (Moore 1986), and the empirical (or in-situ) method (Goulden et al. 1996a). The analytical method aims to derive transfer functions that describe the response of the measurement system to individual causes of attenuation, which is then multiplied by a model cospectrum to estimate the required correction. Early examples of analytical methods dealt with distortions introduced by line averaging of sonic anemometers (Kaimal et al. 1968;Silverman 1968). The net frequency response of the EC system in the analytical method is considered as the product of the several independent response functions that account individually for different causes of spectral attenuation. The resulting net frequency response function is multiplied by a model of the unattenuated cospectra of the scalar and vertical wind velocity to obtain a correction factor (Moore 1986). Several phenomenological models of turbulence spectra and cospectra have been utilized in the analytical method, most notably the models of Kaimal et al. (1972Kaimal et al. ( , 1976. Empirical methods have different assumptions about the net transfer function and the ideal unattenuated spectrum. A different scalar undergoing little to no attenuation is assumed to be a reference. The spectrum of the reference scalar is considered to be the ideal spectrum by utilizing the hypothesis of spectral similarity (Ohtaki 1985). Most commonly, the spectrum of the sonic temperature measured using a sonic anemometer is considered the unattenuated ideal spectrum. The net system response is typically assumed to be a first-order low pass filter (RC filter) with one parameter to estimate, the time constant τ (Goulden et al. 1996b). The filter's time constant is estimated from the ratio of the attenuated scalar spectrum to the ideal spectrum (Ibrom et al. 2007). But also the use of the cospectra is common (Mammarella et al. 2009). Less commonly, τ can be experimentally estimated by adding CO 2 to the analyzer (Goulden et al. 1996b) or by utilizing the observed phase shift (Shaw et al. 1998). The estimated filter is used to simulate the attenuation on the reference scalar (sonic temperature) and obtain a correction factor that can be used directly to correct the attenuated fluxes or parameterized using wind speed and atmospheric stability classes (Goulden et al. 1997;Berger et al. 2001;Ibrom et al. 2007;Fratini et al. 2012). Other methods that do not fit exactly in this classification scheme are the ogive method, which involves comparing the cumulative frequency curves (ogives) of sensible heat flux and the scalar to determine the flux loss (Ammann et al. 2006) and the wavelet approach, in which the variance of the attenuated scalars is adjusted to match the theoretical power law in the inertial sub-range in the wavelet domain (Nordbo and Katul 2013). Both of these methods make similar assumptions regarding the ideal shape of the spectra or the cospectra.
The analytical and the empirical methods both have known limitations. Firstly, the use of ideal reference spectra and cospectral to estimate the attenuation is not well-justified, regardless of whether these ideal spectra and cospectra are obtained from turbulence spectral models or based on the assumption of spectral similarity. Experimental data point at a large degree of variability and significant deviations from the ideal spectra and cospectra (Mamadou et al. 2016). This variability often results in varying correction estimates, depending on the chosen correction method (Reitz et al. 2022;Zhu et al. 2015). Furthermore, although the spectral similarity might hold for a large ensemble of averaged spectra under favorable conditions, there is no reason why the spectra should be similar in each averaging interval, as implied by empirical methods. In addition to the measurement height and meteorological conditions, variations in spectra can be caused by the differences in spatial and vertical distribution of sources and sinks of different scalars and by the contribution of sonic temperature to buoyancy (McBean 1973). Secondly, conventional methods lack a proper treatment for the contribution of the phase spectrum to the attenuated covariance. The contribution of the phase spectrum, also known as the contribution of quadrature spectrum and the out-of-phase contribution to the covariance (Kaimal and Finnigan 1994) is caused by the changes of the phase angle between w and c due to the phase shifting effect of the attenuating filter. This contribution has been acknowledged before (Hicks 1972;Kristensen and Lenschow 1988;Eugster and Senn 1995;Horst 2000). However, it was either ignored (Moore 1986), replaced by an approximation based on simplifying assumptions (Eugster and Senn 1995;Peltola et al. 2021), assumed negligible (Massman 2000;Horst 2000;Mammarella et al. 2009), or, assumed that digitally shifting the time series would compensate for the contribution of phase spectra (Ibrom et al. 2007;Fratini et al. 2012). Finally, the consequences of some assumptions and treatments of empirical methods are subject to ongoing discussion such as model assumptions, noise effect on the spectra, interactions with the time lag correction, and how the spectral transfer functions are applied Aslan et al. 2021).
A direct spectral correction approach for EC measurements remains a desired goal. Few studies have attempted to directly restore the original signal or the (co)spectra. Massman (2000) mentioned that a direct correction is possible but too computationally expensive. Wilczak and Bedard (2004) used simple inverse filtering to correct pressure measurements but eliminated frequencies below a critical limit to avoid noise amplification. (Nordbo and Katul 2013) proposed an empirical wavelet-based approach to restore the original signals by compensating for the estimated variance loss in the wavelet space. More recently, Polonik et al. (2019) have evaluated directly correcting the attenuated cospectra with empirically estimated transfer function. Despite these efforts, a general direct spectral correction approach that optimally restores the attenuated signals in the presence of noise is currently unavailable. The goal here is to develop such an approach.
This paper aims to revise the theoretical basis of the problem of signal attenuation in EC flux measurements and show that the Wiener deconvolution is an optimal solution to estimate the original signal. In order to estimate the size of the problems addressed here, a numerical simulation is utilized to examine the extent of the key limitations of common spectral correction methods that can be resolved by the Wiener deconvolution method. Future 32 A. Emad work will present an extended evaluation of the Wiener deconvolution method using field measurements from multiple sensors.
We start in Sect. 2 by examining the theoretical basis of the spectral attenuation of EC measurements in the framework of linear system theory. We derive an accurate cospectral transfer function which forms the basis for subsequent analysis of the role of phase spectrum and measurement noise in signal restoration. Based on this analysis, the Wiener deconvolution method is proposed as an optimal solution for restoring attenuated noisy atmospheric measurements in Sect. 2.2. Subsequently, we discuss the theoretical aspects of modeling the system response and estimating its parameters within this framework. We then proceed to explore the key limitations of conventional methods in a numerical simulation. The simulation set-up and evaluation are presented in Sect. 3 where the effects of three interrelated issues are investigated: (i) the inappropriate assumptions of the attenuation model, (ii) the contribution of the non-linear phase to the attenuation, and (iii) the consequences of the spectral similarity assumption. We present and discuss the simulation results in Sect. 4. Finally, we discuss the requirements and practical limitations of applying the Wiener deconvolution method to restore attenuated atmospheric scalars.

Theory
We begin by providing a short overview of signal attenuation in EC measurements, followed by the construction of a simple accurate cospectral transfer function. Next, we present the Wiener deconvolution as an optimal solution to the problem of signal attenuation. We then discuss the models of attenuation and the problem of model parameter estimation.
It is worth noting that different concepts in the EC literature are interchangeably referred to as "transfer function". This study uses terminology consistent with linear system theory and signal processing (e.g. Smith (2007)). An overview of the notation used is provided in Table 1.

Spectra and Cospectra of Atmospheric Signals
For an infinite stationary time series c(t) of a scalar, c, (such as an atmospheric constituent) observed over a region of finite-time T , the power spectrum of the scalar c, S cc is defined as (Miller and Childers 2012): where ω = 2π f is the angular frequency, and C(ω) denotes the Fourier transform of c(t). It is customary in EC literature to define the power spectrum in terms of the autocorrelation function of the scalar c over a distance r 1 in the x direction, R(r 1 ). Such a definition can be shown to be equal to the previous definition since R(r 1 ) and the spectrum S(ω) form a Fourier transform pair. Time and space in R(r 1 ) can be interchanged in the previous definition by making use of Taylor's frozen turbulence hypothesis (Kaimal and Finnigan 1994).
The phase spectrum of C is defined as ∠C(ω) ≡ arg(C(ω)) where arg is the argument of C. The area under the spectrum equals the biased variance of the signal, σ 2 c (Parseval's theorem): For two variables such as the scalar c(t) and the vertical wind speed w(t), the crossspectrum Cr cw is similarly defined as: where [C(ω)] * is the complex conjugate of C(ω). The cross-spectrum has a real part and an imaginary part. The real part of the cross-spectrum is the cospectrum C wc and the imaginary part is the quadrature spectrum Q wc . It can be shown that the quadrature spectrum is the Fourier transform of the odd part of the cross-covariance function R wc (r 1 ) and that it makes no contribution to the cross-covariance at (r 1 = 0), R wc (0) (Kaimal and Finnigan 1994). Similar to the spectrum, integrating the cospectrum of w and c, C wc gives the covariance wc: The spectra and cospectra for atmospheric measurements are usually binned, normalized, averaged, or weighted by the frequency to emphasize certain features (Stull 1988). They are also often defined as one-sided (0 to +∞) since the signals we deal with are real-valued. Several methods exist to estimate the true spectral density from a limited number of samples of a random process (Percival and Walden 1993). In EC measurements, a variant of the periodogram method is used to calculate the spectra for half-hourly intervals (Bartlett 1948).

Signal Attenuation
For a system measuring a scalar quantity, the observed signal at time step t, c m (t) can be described as the convolution of the true signal c(t) and the impulse response function of the measurement system h(t): where * denotes convolution. The response of the system is assumed to be Linear and Time-Invariant (LTI). Consequently, the impulse response function h(t) completely characterizes the behavior of the system. The impulse response function h(t) represents the output of the dynamic system in response to a sudden change (impulse) in the time domain. The assumption of linearity implies that the system satisfies the superposition principle, which means that a linear change in the system's input will induce a linear change in the output. Time invariance implies that the system output does not depend on the time of the measurement.
The LTI system can also be characterized by its representation in the frequency domain H (ω). The (complex) frequency response function of the system H (ω) is defined as the Laplace transform of the impulse response function h(t). The convolution from Eq. 5 is equivalent to multiplication in the frequency domain (the convolution theorem), therefore: where C m (ω) and C(ω) are the Fourier transforms of the attenuated signal c m (t) and the true signal c(t), respectively, H (ω) is the frequency response of the system, and · denotes point-wise multiplication.
In the case of EC measurements and assuming the attenuation is only affecting the scalar c, the attenuated (measured) covariance, wc attn is the integral of the real part of the crossspectrum: Here, wc is considered in the general case, but w c can be easily obtained by removing the zero-frequency component (mean of the time series). For notational simplicity, we drop the frequency from the notation (e.g., H · W instead of H (ω) · W (ω)).
We can write Eq. 7 using the amplitude and the phase shift as: where the phase cross-spectrum, ∠Cr cw , is the phase shift between c and w at each frequency, ∠H is the phase spectrum of the filter, |Cr cw | is the amplitude cross-spectrum, and |H | is the magnitude response of the filter. The previous equation can be further decomposed to relate the attenuation to the spectra of c, and w as: where ∠C and ∠W are the phase spectra of c(t) and w(t), respectively. A cospectral transfer function is often defined in the EC literature as the transfer function that will produce the attenuated cospectra when multiplied by the true cospectra. The earliest use of such a function is the formulation of Moore (1986) where he used the magnitude response of the filter, |H | (in his notation √ H ) to express the attenuation on the cospectra. Eugster and Senn (1995) approximated the phase induced attenuation for an RC filter as |H |, thereby using |H | 2 as a cospectral transfer function which was further supported by Horst (1997Horst ( , 2000. To evaluate these formulations, we construct an accurate cospectral transfer function. We rearrange Eq. 8 to isolate the cospectra: Therefore, the function Λ(ω) defined when ∠Cr cw = π 2 + n π as: can be considered an accurate cospectral transfer function that will produce the desired attenuation when multiplied by the cospectrum. The function Λ(ω) highlights the importance of considering the phase shift introduced to individual frequency components if it is to be used for correcting attenuated cospectra. This raises the question of whether simply shifting the time series in time would correct for the phase shift. To answer this question, it is important to consider the concept of group delay. Group delay represents the time delay each frequency component undergoes due to phase shift and is defined as the derivative of the phase spectra with respect to frequency (Mandal and Asif 2007). If the phase shift is a linear function of the frequency, all frequency components would be shifted equally and the net effect on the time series would be similar to shifting the time series by a constant time lag. Therefore, the time lag correction procedure, which is commonly applied in EC data processing and involves shifting the scalar time series to maximize its correlation with the vertical wind velocity, can only correct for linear phase shifts at best. Any remaining effects on the cospectra due to the nonlinearity of the phase will continue to bias the fluxes. We refer to these remaining effects on the flux as "nonlinear phase contribution" to the flux attenuation.

Signal Reconstruction Using the Wiener Deconvolution Method
Using the properties of convolution, recovering the original signal seems conceptually simple. The original signal can be restored by dividing the Fourier transform of the attenuated signal by the frequency response function of the filter and then taking the inverse Fourier transform. Rearranging Eq. 6, we can write: and the original signal will be: However, the problem with this approach, commonly called "naive deconvolution", is that the noise in the signal will be amplified, particularly in regions of the spectrum where the attenuation is high and the signal is low. This results in a recovered signal dominated by noise and is often unrecognizable. In most real-world applications, the measured signal is contaminated by noise. Therefore, the observed signal c m (t) can be more accurately described as a sum of the attenuated signal and an additive noise term ε(t) independent of c(t): In the presence of noise, we wish to get the best estimate of the original signal c(t) denoted asĉ(t). The goal is then to find the optimal filter G that, when applied to the measured signal, results in the best estimate of the original signal such as: This objective can be expressed as minimizing the mean squared error between the original signal and the estimated signal. The error that needs to be minimized can be written as: where E denotes mathematical expectation. We substitute the estimated signalĈ with its expression from Eq. 15 and the frequency domain representation of Eq. 14. We can then write the error as: where E is the Fourier transform of the noise ε(t). Expanding the previous equation, we find: Considering that the noise is assumed independent of the signal, we have E C · E * = E C * · E = 0. Substituting E|C| 2 and E|E| 2 with their corresponding mean power spectral densities,S(ω) andS N (ω), we obtain: This optimization problem was solved by Wiener (1964). Further details of the different aspects of the solution are also provided by Meditch (1969) and Kay (2006). The minimum value of the error in Eq. 19 is found by taking the Wirtinger derivative with respect to G and setting it equal to 0. The solution that minimizes the error is found to be: Rearranging the previous equation, we find the formulation to the optimal deconvolution filter as: where SNR is the mean signal-to-noise ratio spectrum defined as: The Wiener deconvolution method combines knowledge about the noise spectra and the system response to produce the best estimate for the original signal c(t). The estimated signalĉ minimizes the least-squares error between the observed signal c m (t) and the original signal c(t).
After constructing the deconvolution filter G(ω), the estimated signalĉ is obtained by multiplying the deconvolution filter by the Fourier transform of the measured signal and taking the inverse Fourier transform of the product:

Modeling the System Response
The knowledge of the response function of the linear time-invariance (LTI) system is essential for the application of the Wiener deconvolution method and conventional spectral correction methods. The response function can be derived from the mathematical description of the system, usually by finding the differential equations that govern the system and applying the Laplace transform to them. However, a description of the physical dynamics of the system is not always possible. In such a case, the system can be assumed to follow a certain simplified response approximating the true system behavior. Therefore, we distinguish between models developed from the underlying physical dynamics of the system and empirical models approximating the behavior of the system.

Empirical Models
Empirical models for the system response function are assumed based on the observation of the magnitude response |H | of the system. The simplest form of these models is the first order low-pass filter, also known as the RC filter and the infinite impulse response filter. The first order system (such as an RC circuit) is described using the following differential equation: where τ is the characteristic time constant of the system. The frequency response function of this system in the Laplace domain is given by: where s is a complex variable known as the Laplace variable. The magnitude squared response of the RC filter is given by: The inverse Laplace transform of Eq. 25 gives the impulse response function of the RC filter in the continuous time domain: The low-pass RC filter is well-suited to model certain sensors with simple responses such as the gain of a resistor-capacitor circuit. However, the overall system response of an EC system is typically more complex and involves multiple sources of attenuation which makes it difficult to be adequately described by such a simple response function. Despite this, the magnitude response of the RC filter has been extensively used to model the attenuation caused by EC systems in numerous studies (Moore 1986;Eugster and Senn 1995;Goulden et al. 1996b;Horst 1997;Berger et al. 2001;Ibrom et al. 2007;Mammarella et al. 2009;Fratini et al. 2012;Peltola et al. 2021). However, the RC filter was found to be a poor approximation for the observed attenuation on many EC systems (de Ligne et al. 2010) which motivated the use of modified functions to model the observed spectral attenuation ratio such as the Gaussian function (Aubinet et al. 2001), or the Gaussian function with an adjustable exponent (de Ligne et al. 2010). These functions, while effective for modeling the spectral shape, do not qualify as valid filters for application as convolution as they do not include information about the phase shift.

Process-Based Models
In process-based models, the system's response function is developed by describing the physical mechanism causing the attenuation. A good example relevant to this study is the process-based physical model for the behavior of water vapor in tubes. This model describes the attenuation of scalars experiencing absorption and desorption as a diffusion process in addition to Taylor dispersion (Kekäläinen et al. 2011; The model describes the mechanism in which water drops adsorb and condense on tube walls in periods of high humidity until an equilibrium is reached. Then, condensed liquid water evaporates in less humid conditions enhancing the content of air with water vapor. This mechanism is formulated as an advection-diffusion equation and its solution under appropriate boundary conditions for the breakthrough curve gives the frequency response function of the system. The sorption model was tested in laboratory and field conditions and was shown to be a good representation of the attenuation of water signals in EC systems (Nordbo et al. , 2014. The frequency response function of the sorption model is given by the solution in the Laplace domain: where μ, λ, and κ are dimensionless variables that characterize the response function. The parameter μ is related to molecular diffusion in air and Taylor dispersion and can be estimated from the system dimensions using Taylor theory (Taylor 1954). The parameter λ is related to the relative area available for condensation and evaporation. The parameter κ is related to the length scale attached to the diffusive evaporation process. The model parameters λ and κ are not easy to know a priori and need to be estimated from the data. They were found to correlate with the material, the age, and the dirtiness of tubes and particle filters (Nordbo et al. 2014). The inverse Laplace transform for Eq. 28 gives the impulse response to the system h(τ ): The function F and a detailed account of this solution are provided in Kekäläinen et al. (2011, Eq. 65)

System Response Parameter Estimation
As discussed earlier, the parameters for the system frequency response are often unavailable and need to be estimated from the data. The magnitude squared of the response function can be found from Eq. 6 to be the ratio of the attenuated spectrum to the true spectrum: In practice, the true spectrum is not known and is substituted with that of another scalar, based on the assumption of spectral similarity. The spectra typically used for parameter estimation are long-term averaged spectra chosen during favorable atmospheric conditions. Therefore, this use is rather valid and different from using the spectral similarity hypothesis for estimating the loss on the short scale of averaging intervals.
Once |H | 2 is known, the filter parameters can be estimated using non-linear regression (Ibrom et al. 2007;Fratini et al. 2012;Nordbo et al. 2014). The use of the cospectra to estimate |H | or |H | 2 e.g. (Aubinet et al. 1999(Aubinet et al. , 2001Mammarella et al. 2009) will create biased estimates due to the effect of the phase shift as can be seen from Eq. 10.
In the presence of noise, the measured spectra contain the spectra of the noise. The measured spectrum of noisy signal |C m | 2 can be written as: Here, the measured spectrum |C m | 2 is the sum of the true spectrum |C m | 2 , the noise spectrum |E| 2 , and a correlation term between the noise and the signal 2 [CE * ]. The correlation term will converge to its expected value of zero given an additive uncorrelated noise. Therefore, the noise spectra can be subtracted from the measured spectra prior to fitting. A discussion of the effect of noise on model parameter estimation is provided in Aslan et al. (2021).

Numerical Simulation
The simulation aims to estimate the extent of the problems that can be potentially addressed with the Wiener deconvolution method. This goal is realized by analyzing the sensitivity of the corrected flux to the individual and the combined effects of three related issues that are intrinsic to conventional methods (i) inappropriate system response, (ii) the contribution of the non-linear phase to the attenuated flux, and (iii) the consequences of the spectral similarity assumption. The simulation is limited to the common practices that are believed to contribute most to biases and uncertainties in the corrected fluxes. It does not aim to provide a comprehensive evaluation of all variations of EC spectral correction methods. Starting from a high-frequency dataset of EC measurements, two models for attenuation were considered: (i) a first-order low pass filter similar to Eq. 25 commonly used in EC literature (hereafter referred to as RC), and (ii) a physical process-based model for sorbing scalars shown in Eq. 28. Both filters were discretized using the bilinear transform and applied to each averaging interval in the frequency domain, similar to Eq. 7. The process was repeated by varying the filter parameters to cover a range of values that are similar to field conditions relevant to eddy covariance setups. For each combination of the sorption model parameters, a matching RC filter that approximates its magnitude response was considered. The time constant τ of the RC filter was estimated using non-linear least squares, such that it minimizes the mean square error between the magnitude response of the RC filter and the sorption model. Both of these models are shown in Fig. 1 for a few selected parameters.
In addition to direct convolution, the attenuation was also applied by multiplying the cospectra by the magnitude response of the filter without considering the phase contribution similar to the approach used in Moore (1986), Fratini et al. (2012), Nordbo et al. (2014). The attenuated flux in this approach is calculated as: The results of the two different attenuation methods were evaluated by computing correction factors for all the scalars using three strategies that are analogous to the ways conventional 40 A. Emad methods are applied: (i) F comp , which is the ratio of the true flux to the attenuated flux, (ii) F tl , which is the correction factor obtained after removing the introduced time lag in the attenuated flux by covariance maximization, and (iii) F no.φ , which is the correction factor of fluxes that were attenuated without the phase response using Eq. 32. The investigation of the three issues proceeded as follows 1. Inappropriate response function: The use of the RC filter to model the response of the attenuation process is prevalent in the EC literature. The consequences of this simplifying assumption were investigated for the range of parameters of the sorption model. Scalars were attenuated using both the sorption model and the RC filter. The ratios of the attenuated fluxes were compared with and without the correction for time lag, using F comp and F tl , respectively. 2. Contribution of the non-linear phase: As discussed in the theory section, time lag correction can only account for the contribution of the linear phase to the attenuation. Here, the contribution of the non-linear phase is investigated for both the RC and the sorption filters by comparing the attenuated flux with no phase F no.φ to the time lag corrected attenuated flux F tl . The difference between the two is considered the non-linear phase contribution. 3. Consequences of spectral similarity assumption: The assumption of spectral similarity is used in all empirical methods. The systematic bias and uncertainty due to this assumption were investigated by comparing the attenuated fluxes of sonic temperature and the scalar under examination using both the RC and the sorption filters.
The interactions between these three issues were explored using three restoration strategies closely imitating conventional empirical methods.
-R no.φ restoration strategy: This strategy involves removing the time lag of the attenuated scalar through covariance maximization, followed by correcting the flux using the correction factor obtained from simulating the attenuation on sensible heat flux while ignoring the phase F no.φ . This strategy assumes that the effect of phase shift is eliminated by time lag correction, and therefore simulates the attenuation using Eq. 32. A similar approach can be found in Fratini et al. (2012), Nordbo et al. (2014). The strategy R no.φ is prone to biases arising from spectral similarity assumption and the non-linear phase contributions on the restored fluxes. -R comp : In this strategy, the correction factor F comp of sensible heat flux is used to correct the attenuated flux without correcting for phase shift. The R comp strategy is prone to biases from the combination of spectral similarity, linear phase, and non-linear phase to scalar attenuation. This strategy is not commonly used and is included only to facilitate the comparison of the combined attenuation effect. -R tl : In this strategy, time lag removal by covariance maximization is applied to both the attenuated scalar and the simulated attenuation on sensible heat flux. Then, the correction factor F tl of sensible heat flux is used to correct the time-lag corrected attenuated flux. This approach is similar to the one employed by Goulden et al. (1997) These correction strategies cover most of the varieties of empirical spectral correction methods. However, certain variants were not included, such as the evaluation with |H | 2 , as well as methods that depend entirely on parametrization (e.g. (Ibrom et al. 2007)). This is because these variants could introduce additional uncertainty to the evaluation.

High-Frequency Dataset
A high-frequency EC dataset was used as an input for the simulation. Measurements were taken over an ideal flat agricultural field of the Thünen Institute, located at 52.297 N, 10.449 E in Braunschweig, Germany. A complete description of the site and instrumentation is available elsewhere (Emad and Siebicke 2023). Data used for the simulation spanned 10 days from 15 June 2020 until 26 June 2020 and comprised 512 30-min averaging intervals. The variables used in the simulation were the three-dimensional wind velocity components measured using a sonic anemometer (uSonic-3 Class A, Metek GmbH, Germany) and gas concentration measurements of CO 2 and H 2 O obtained from an open-path infra-red gas analyzer (LI-7500A, LI-COR Biosciences GmbH, Germany).
The raw data were prepared for the simulation by applying common statistical screening for data quality issues (Vickers and Mahrt 1997). The initial time lag between the wind and the scalar time series was compensated using covariance maximization. The wind coordinates were rotated using the planar fit method (Wilczak et al. 2001), and any remaining residualw was removed.

Simulation Parameters
The parameters adopted for the sorption model were based on empirically determined parameters from experimental data and are considered typical for closed-path eddy covariance setups (Nordbo et al. 2014). The parameter μ was set at a fixed value, μ = 0.001 since its value was found to have a minor effect on the shape of the spectrum within its reasonable range compared to λ and κ. The range of values considered for λ and κ was from 10 −5 to 2.5 for both parameters. The tube length considered was fixed at L = 1 m, and the flow velocity was set to v = 1 m s −1 . The simulation covered a total of 5000 treatments.

Simulation Results Quality Control
Fluxes were calculated for 30-min averaging intervals. Then, quality checks were applied to computed fluxes and simulation results. Averaging intervals with low developed turbulence, low u * , and low fluxes of CO 2 , sensible and latent heat were excluded from the analysis similar to typical quality checks of EC fluxes and spectra (Sabbatini et al. 2018). After quality checks, 34% of high-quality averaging intervals remained for the final evaluation. Additionally, for simulation treatments that required time lag optimization, unrealistic estimated time lags that produced spurious correlations were eliminated.

Simulation Results
The simulation results are limited to fairly ideal conditions and high-quality fluxes. The input data were taken from a flat site where the variability of spectra is not expected to be very high. Furthermore, moderate simulation parameters were chosen for the sorption model. In particular, the ratio of the tube length to flow velocity was set to 1. Finally, the quality checks employed left only very well-behaved periods to be considered in the evaluation. Therefore, we note that the estimates obtained from the simulation strictly relate to the current conditions and are expected to vary for other scalars and sites, with more complex sites likely giving more variable results.

Inappropriate Response Function
Using an RC filter to model the response of the sorption process introduces systematic bias to the fluxes. If the filter-induced time lag is not corrected, the sorption filter will always attenuate the signal more than its approximated RC counterpart (Fig. 2a). This is not surprising as the time lag caused by the strong phase shift of the sorption filter will lead to the degradation of the covariance. On average the attenuation caused by the sorption model was 1.3 to 1.7 times higher compared to the RC low pass filter.
Correcting the linear phase shift by time lag optimization reduces the attenuation ratio considerably (Fig. 2b). However, the remaining effects of the non-linear phase shift and the differences in the filter's magnitude response lead to a systematic bias reaching up to 15% of the recovered fluxes in either direction (underestimation or overestimation). The magnitude and direction of the bias are mainly controlled by the shape of the filter's magnitude response. The inability of the RC filter to accurately approximate the shape of the sorption filter leads to uneven attenuation at different frequencies. The mismatch between the two filters is manifested by two features (Fig. 3). The first feature is the different steepness of the two filters (roll-off) as seen from the "over-est" scenario, here the RC filter is causing stronger attenuation due to having a slower roll-off (less steep). The importance of the filter's roll-off in the attenuation is due to the unequal distribution of energy in the spectra of turbulence. The second feature is the dip in the magnitude response of the sorption filter observed around the 0.1 Hz frequency in the treatment "under-est". The dip is not possible to emulate using the simple RC filter and its position is controlled by the model parameters κ and λ (Nordbo et al. 2014).
The observed spread for individual averaging intervals around the mean is shown for one set of parameters in (Fig. 2c). The spread translates to increased random uncertainty in the restored fluxes and is mainly driven by the variability in the distribution of spectral energy within different averaging intervals. The spread shows dependence on the discrepancy between the RC and the sorption filter, therefore the largest bias is accompanied by the largest variability. Fig. 3 The effect of the filter magnitude response on the attenuation ratio outcome after correcting the effect of linear phase shift. a Magnitude responses of the sorption model (solid lines) and the RC filter approximating its response (dashed) for three values of λ and a fixed value of κ, κ = 1.02 that represent three possible outcomes: over-est: the RC filter is causing stronger attenuation to the flux (λ = 0.36). under-est: the sorption filter is causing stronger attenuation to the flux (λ = 2.3). min-bias: minimum bias between the sorption and the RC filter (λ = 0.86). b Ratio of attenuated CO 2 fluxes using a process-based sorption model to fluxes attenuated using a matching RC filter

Phase Non-linearity
The contribution of the non-linear phase to the attenuation of the fluxes can reach up to 5% of the flux for both the RC filter and sorption filters (Fig. 4). The average non-linear phase contribution shows an increase with increased attenuation. This non-linear contribution was considered as the difference between the attenuated fluxes using cospectra adjustment without a phase to the attenuated flux after removing time lag using covariance maximization. This procedure assumes that the time lag correction will remove the effects of the linear phase. The mechanism by which the non-linear phase shift contributes to the attenuation is by causing different phase shifts to different frequency components of the signal. The net effect on the time series of the scalar is a nonuniform distortion that is more likely to reduce the correlation between w and c. The correlation loss can not be corrected by digitally shifting the time series. Conversely, it will likely reduce the effectiveness of the time lag correction procedure.

44
A. Emad   Fig. 5 Effects of the assumption of spectral similarity on restored fluxes in conventional method. Ratio of recovered to attenuated flux when a correction factor obtained by simulating the attenuation effect on sonic temperature is used to correct CO 2 fluxes attenuated with the same filter a if an RC filter is used. b if the sorption filter is used (κ = 1.83). Gray lines show the recovery ratio for individual averaging periods, red line is the mean for the whole dataset, and the dashed lines represent the mean ± 2 SD Fig. 6 a Spectra of three scalars: CO 2 , H 2 O, and Ts demonstrating spectral similarity of averaged spectra. Gray lines show spectra from individual averaging periods, dashed lines represent the mean ± 2 SD. b Variability of the spectra across frequencies, calculated as the ratio of the standard deviation to the mean for binned frequencies (n = 176)

Limitations of Spectral Similarity
The assumption of spectral similarity introduced a small averaged systematic bias of up to 3% of the corrected fluxes if other aspects of the correction were treated optimally. However, a larger increase was observed in the uncertainty of the corrected fluxes. The standard deviation of restored fluxes increased up to 10% of the flux value (Fig. 5). The observed systematic bias is likely caused by the slight differences between the averaged spectra shapes as well as the differences in individual averaging intervals spectra (Fig. 6). The variability in the spectra from one averaging interval to another is not uniform across frequencies and different scalars (Fig. 6b). The observed variability in the outcome of the applied filter to different scalars showed a very weak dependence on wind speed and stability classes. While the wind speed and stability explain part of the variability of the spectra, large variability remained after controlling for wind speed and stability (not shown). This variability is potentially due to the distributions of sources and sinks in addition to buoyancy in the case of sonic temperature. Fig. 7 Combined effects of phase contributions and the similarity assumption on the recovered fluxes when using an RC model for attenuation and recovery. Three conventional correction strategies are shown where the attenuated scalar is corrected by simulating the attenuation on sensible heat flux. R comp the attenuation is simulated as convolution (magnitude and phase). R no.φ assumes the phase shift is corrected by time lag correction and simulates the attenuation on the cospectra without a phase shift. R tl attenuation is simulated as convolution but linear phase shift is corrected separately by covariance maximization. a Solid lines represent the mean systematic bias averaged for the whole dataset. Dashed lines show the variability in restored fluxes represented as the mean ± 2 SD. b Variability in restored fluxes for individual averaging intervals shown as the standard deviation of the recovery ratios

Combined Effects on the Corrected Fluxes
All three conventional methods for restoring the fluxes introduced varying degrees of systematic bias and increased uncertainty to the recovered fluxes (Figs. 7, 8). With regard to systematic bias, the strategy R no.φ was by far the worst performer when using either the RC or the sorption filter, followed by R comp , then R tl . The same dynamics are observed when using an RC filter to model the sorption process but with much larger biases.
These results can be understood by considering the three attenuation sources: the damping of the spectra due to the filter's magnitude, the attenuation due to the linear phase shift, and the attenuation due to the non-linear phase shift. The different responses to these attenuation sources between the sonic temperature and the scalar as well as whether a correction method accounts for one or more attenuation sources are what drive the observed biases and uncertainties. R no. phi is applied to the time lag corrected fluxes and only corrects the attenuation attributed to the magnitude of the cospectra while ignoring the contribution of the non-linear phase. Thus, correcting the attenuated flux (that still has the additional attenuation of the non-linear phase) using R no. phi will introduce the largest bias. This is still valid, even though no phase similarity between Ts, and CO 2 was observed, potentially because accounting for the non-linear phase reduces the systematic bias by randomizing the systematic error leading to increased variability in the result, as seen in Fig. 7.
However, including the linear phase in the correction factor does not improve the recovery as evident by the better performance of the R tl strategy. Accounting for both linear and nonlinear components of the phase with R comp has inferior performance compared to accounting for the non-linear phase and removing the effect of the linear phase using time lag correction in the R tl strategy. The reason is likely the difference in the contribution of the linear phase 46 A. Emad

Fig. 8
Combined effects of phase contributions, spectral similarity assumption, and the use of a simplifying model on the recovery of fluxes attenuated with a sorption model. Two conventional correction strategies are shown where the attenuated scalar is corrected by simulating the attenuation on sensible heat flux. Each of the two methods is evaluated in the case of using the correct sorption model for simulating the loss (first row) and when a simplifying RC filter is used (second row). R no.φ assumes the phase shift is corrected by time lag correction and simulates the attenuation on the cospectra without a phase shift. R tl attenuation is simulated as convolution but linear phase shift is corrected separately by covariance maximization. a Using a sorption filter for correction. b Using a simplifying RC filter for correction. Lines represent the mean systematic bias averaged for the whole dataset. Line shapes represent the method used, solid lines for R no.φ and dashed for R tl . b and c Variability in restored fluxes for individual averaging intervals for R no.φ and R tl respectively. e and f Variability in restored fluxes for individual averaging intervals for R no.φ and R tl respectively between sonic temperature and the scalar, thus it is better to remove them separately rather than including them in the correction.

Limitations of Conventional Spectral Correction Methods
The analysis presented in this paper demonstrates that both analytical and empirical conventional spectral correction methods have inherent shortcomings due to the inadequate representation of the true cospectra, the inappropriate simplified models, and ignoring the nonlinear phase effects. Consequently, these methods are likely to produce unreliable approximations for flux corrections.
To examine the limits of conventional analytical correction methods, we formulated in Eq. 11 an accurate cospectral transfer function that does not assume a specific shape of the attenuation filter. The equation shows that the use of a cospectral transfer function for estimating the attenuation requires the knowledge of the phase shift between c and w represented by the phase spectrum, ∠Cr cw and the phase shift introduced by the attenuation filter. The phase spectra have no available models, unlike the spectra, and are not expected to be similar among different scalars. The equation suggests that the commonly used |H | is only correct if the filter has no phase shift, an unachievable condition for causal physical processes (Smith 2007). Alternatively, |H | can be used if time lag correction eliminates the impact of the phase shift. However, as explained in the theory section, shifting the time series can only effectively correct constant group delay that is caused by a linear phase shift. In nearly all practical scenarios, the phase is non-linear, and the remaining non-linear phase distortion will lead to biased corrected flux. This occurs because a non-linear phase will shift the individual frequency components of the cospectrum at varying rates and thus will change the signal in a way that cannot be predicted solely from the spectrum. This generally supports the recommendation of using |H | after time lag correction ). However, this approach is a first-order approximation and is not exact. For similar reasoning, estimating |H | as the ratio of attenuated to ideal cospectra is considered an approximation.
In addition to the incomplete representation the phase distortion, the use of modeled spectra and cospectra in the analytical method is likely to increase the uncertainty and systematic biases in corrected fluxes. This is because these models are generally simplistic and fail to capture the variability of the cospectral shapes across time and site-specific conditions.
Empirical methods face similar challenges related to the accuracy of the ideal spectrum and the models of the attenuation filter. However, a key difference is that instead of simulating the attenuation on a model spectrum, empirical methods assume that the ideal spectrum originates from another scalar, based on the assumption of spectral similarity.
The extent of these issues on empirical methods was investigated in a numerical simulation. The results revealed that the largest systematic biases arose from the use of simplified response functions, followed by ignoring the non-linear contribution of the phase shift, and finally, from assuming spectral similarity on the scale of averaging intervals. When combined, these effects can lead to more than 10% systematic bias in the corrected fluxes even for small filter time constants, as demonstrated by the correction strategy R no.φ (Fig. 8d). Many EC flux correction methods, including the widely used in-situ methods (Ibrom et al. 2007;Fratini et al. 2014), adopt a similar approach and are expected to produce similar systematic biases when applied to correct water vapor fluxes in closed-path EC systems. The use of spectral similarity assumption has been found to increase the uncertainty in restored fluxes, in general this increase depended on the magnitude of the attenuation and reached values of more than 10% (1 SD).
In summary, the limitations of conventional methods result from the incomplete representation of the transfer function, the use of inappropriate response functions, and the reliance on indirect simulation of the attenuation on a proxy scalar or model cospectra rather than directly correcting the signal.

Signal Restoration Using the Wiener Deconvolution Method
The current study demonstrated that the Wiener deconvolution method is an optimal solution to restore attenuated signals of atmospheric scalars. The Wiener deconvolution method offers the best estimate of the original signal from noisy measurements by minimizing the mean squared error between the true and the measured signals. As a direct method, the Wiener deconvolution avoids many of the shortcomings in conventional methods that arise from the assumption of an ideal cospectrum or spectral similarity and brings many advantages for flux calculation. Unlike conventional methods that are often limited by atmospheric conditions that influence the ideal cospectrum, the performance of the Wiener deconvolution method is not dependent on the properties of the input signal once the attenuation model parameters have been estimated. Additionally, the direct restoration of the scalar signal in Wiener deconvolution allows for obtaining corrected statistics of the original scalar, such as variance and the integral time scale. Furthermore, the Wiener deconvolution completely corrects the phase distortions, making it independent from time lag correction procedures that may still be applied after signal estimation to correct for time lag introduced by sources other than the system response.
The successful application of Wiener deconvolution requires meeting three criteria: (i) knowledge of the system frequency response function, (ii) the response function must be linear and time-invariant, and (iii) the knowledge of the mean signal-to-noise ratio spectrum. The knowledge of the system's frequency response function requires modeling the physical process leading to the attenuation. While models for certain processes are available, such as the sorption model shown earlier, other causes of attenuation might not have readily available models. The analytical models developed for conventional methods, such as line averaging and sensor separation, are good candidates for use with the Wiener deconvolution method. In this regard, the Wiener deconvolution method offers a chance to unify analytical and in-situ methods where physics-based models are applied to the observed spectrum instead of the modeled one.
Once the frequency response function is known, it needs to satisfy the conditions of linearity and time invariance. These assumptions are not specific to the Wiener deconvolution method, although they are not stated explicitly for conventional methods, they are implied by expressing the spectral attenuation as convolution e.g. (Moore 1986;Massman and Clement 2005). Linearity implies that for a linear combination of inputs, the system produces a linear combination of outputs. This idealization is rarely met by physical systems. However, in most cases, the response function can be assumed to be linear or approximated as such without significant error. The assumption of time invariance means that the system's output does not depend on the time the filter is applied. Time invariance is not always guaranteed, as demonstrated by the sorption model shown earlier, where the model parameters λ and κ are dependent on relative humidity, which can vary during the averaging interval when the correction is applied. The violation of time invariance will lead to biases in the restored signal. This issue can be mitigated by applying the correction over shorter intervals where the model parameters are assumed to be constant. However, there is a limit on how fine the time resolution can be as the signal can not be sharply defined in both time and frequency domains (Gabor limit) (Downey 2016). The application of the Wiener deconvolution method also requires finding the signal-to-noise (SNR) spectrum, which requires the knowledge of the power spectral density of the signal before attenuation. This information is usually unavailable. However, knowledge of the approximate power spectral density of the true signal is usually sufficient. For atmospheric scalars, the ensemble average of the sonic temperature spectra or the modeled spectra of Kaimal et al. (1976) should be sufficient for most cases. The SNR spectrum is used to suppress the noise in the restored signal by balancing the trade-off between amplifying the signal and suppressing the noise. An imperfect SNR spectrum may result in noisier estimates, but is unlikely to introduce systematic errors to the restored fluxes. The noise behavior of the sensor can be statistically estimated from instrument measurements if not known beforehand. For instance, experiments can be performed to measure the variance of a measured signal of known concentration in the case of gas analyzers. In cases where experimentation is not feasible and the noise is assumed to be white, the noise spectrum can be estimated by extrapolating the slope in the inertial subrange to zero and attributing the remaining measured power at the higher range of the spectrum to the noise.

Limits of Signal Recovery
The recovery of attenuated frequencies is only possible if attenuated frequencies do not fall below the detection limit of the measurement apparatus. If any frequency components are attenuated to zero or below the noise level, the information is lost and can not be recovered. In such a case, the contribution of the lost frequencies to flux must be estimated using a different source of information, which will likely increase the flux uncertainty. Therefore, the conditions that result in the elimination of frequencies should be avoided by adjusting the system design and operation. For example, in the case of water vapor in tubes, the response function can be solved for the lowest magnitude resolvable by the measuring device to determine safe thresholds for controllable parameters such as tube length and flow rate.

Conclusions
Reliable estimates of the atmospheric exchange require an optimal recovery of attenuated fluxes due to signal frequency losses. The analysis in this paper showed that conventional spectral correction methods introduce systematic biases and increased random uncertainty to the corrected fluxes. These biases and increased uncertainty are due to using simplified response functions, ignoring the contribution of non-linear phase shift, and the assumptions of spectral similarity and model spectra. A new approach based on Wiener deconvolution was shown to be an optimal solution to directly restore attenuated signals without relying on the assumptions of spectral similarity for individual averaging intervals. The Wiener deconvolution approach has wide applications for micrometeorological methods where an ideal spectrum can not be assumed or when the spectrum is distorted due to aliasing, such as in the case of eddy accumulation methods (Emad and Siebicke 2023), flux profile measurements (Kaimal 1986), and disjunct eddy covariance method (Rinne and Ammann 2012;Hill et al. 2017).
More research is currently underway to evaluate the impact of using the Wiener deconvolution method for a larger dataset covering multiple sites on yearly flux estimates and energy balance budget, and to compare it against conventional spectral correction methods.

A. Emad
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A: List of Symbols
See Table 1.