Advances in the True Eddy Accumulation Method: New theory, implementation, and field results

The true eddy accumulation method (TEA) provides direct measurements of ecosystem-level fluxes for a wide range of atmospheric constituents. TEA utilizes conditional sampling to overcome the requirement for a fast sensor response usually demanded by the state-of-the-art eddy covariance method (EC). However, the assumptions and conditions required for the TEA method are often not met. Here we explore the limitations caused by the assumption of zero mean vertical wind velocity during the averaging interval and by the fixed accumulation 5 interval. We extend the theory of TEA method to non-zero vertical wind velocity by employing information about the scalar transport. We further derive a new method with adaptive time varying accumulation intervals. The new method, termed short-time eddy accumulation (STEA), was successfully implemented and deployed to measure CO2 fluxes over an agricultural field in Braunschweig, Germany. The measured fluxes matched very well against a conventional EC system (slope of 1.05, R of 0.87). We 10 provide a detailed description of the setup and operation of the STEA system in the flow-through mode, devise an empirical correction for the effect of buffer volumes, and describe the important considerations for the successful operation of the STEA method. The new theory developments reduce the bias and uncertainty in the measured fluxes and create new ways to design eddy accumulation systems with finer control over sampling and accumulation. The results encourage the application of TEA and 15 STEA for measuring fluxes of more challenging atmospheric constituents such as reactive species as well as other constituents where no fast gas analyzers are available.


Introduction
Micrometeorological methods provide non-invasive, in situ, integrated, and continuous point measurement for ecosystem fluxes on a scale ideal for ecosystem study (Baldocchi et al., 1988;Baldocchi, 2014). Among micrometeorological methods, eddy 20 covariance (EC) has become the de-facto method for measuring ecosystem fluxes for the past 40 years. The EC method is the most direct micrometeorological method. It is also relatively easy to set up and operate. These features have led to the wide use and adoption of the EC method at hundreds of sites worldwide, including several regional and global flux measurements networks such as ICOS and FLUXNET (Hicks and Baldocchi, 2020).
The EC method depends on the fast measurement of vertical wind velocity and the scalar concentration (such as an atmospheric constituent). The requirement for fast measurement frequency (10 to 20 Hz) limits the application of the method to a handful of atmospheric constituents where fast gas analyzers are available. Alternative methods that work for slow gas analyzers include: (i) signal downsampling methods (Lenschow et al., 1994) such as disjunct eddy accumulation (Rinne et al., 2000a;Turnipseed et al., 2009) and disjunct eddy covariance (Rinne and Ammann, 2012), and (ii) indirect methods such as flux gradient methods (e.g Rinne et al. (2000b)) which depend the Monin-Obukhov similarity theory (Monin and Obukhov, 1959) 30 and the relaxed eddy accumulation (REA) based on the assumption of flux-variance similarity (Businger and Oncley, 1990) The true eddy accumulation method (TEA) (Desjardins, 1977) is the most direct and mathematically equivalent alternative to eddy covariance. Unlike EC, the TEA method requires the concentration measurements to be carried out once every averaging interval (30 minutes). For a long time, the development of the TEA method was hindered by the difficulty of fast air flow rate control and the strict operational requirements (Businger and Oncley, 1990;Hicks and McMillen, 1984). A recent improvement 35 to the TEA method used a new type of mass flow controllers, online coordinates rotation, and several online treatments of the signal to overcome important limitations of the method's applicability (Siebicke and Emad, 2019). The new system showed a good match with a reference eddy covariance system with coefficients of determination of up to 86% and a slope of 0.98.
The non-necessity of high-frequency measurements of the scalar in TEA, although a major advantage, introduces multiple difficulties. First and foremost, exact equality between EC fluxes and TEA fluxes is only possible when the mean vertical wind 40 velocity is assumed to be zero during the flux averaging interval. This assumption is almost never met under field conditions and the residual vertical mean velocity contributes to the error in the flux. Nonzero mean vertical wind velocity is a source of error for all eddy accumulation methods, including TEA (Hicks and McMillen, 1984), relaxed eddy accumulation (REA) (Pattey et al., 1993;Businger and Oncley, 1990;Bowling et al., 1998), and disjunct eddy accumulation (DEA) (Turnipseed et al., 2009). The reported bias on the flux due to nonzerow varied with different studies and accumulation methods. For TEA, Hicks 45 and McMillen (1984) recommended thatw should not exceed 0.0005 σ w if accumulated mass is measured and 0.02 σ w when concentrations are measured directly. Turnipseed et al. (2009) reported that a mean vertical wind bias of ± 0.25 σ w lead to ± 15% errors in the flux using the disjunct eddy accumulation method. Values reported for the REA method showed that a loss of approximately 5% of the flux due to aw of 0.20 σ w (Pattey et al., 1993), which agrees with the recommendations of Businger and Oncley (1990). Additionally, the absence of high-frequency information means any decision on air sampling 50 such as flow rate and sampling direction is final. The lack of high-frequency information also implies that sample accumulation should happen on a time scale similar to the flux averaging interval (30 to 60 minutes). These limitations impose restrictive design considerations related to the size and function of sample accumulation reservoirs. They also dictate that the sampling apparatus needs to accommodate a large dynamic range to cover the range of wind velocities during the flux averaging interval.
In this paper, we address these limitations of the TEA method. First, we revise the theory of true eddy accumulation and 55 extend the TEA equation to non-ideal conditions. The new generalized TEA equation allows obtaining TEA fluxes equivalent to fluxes measured with EC when the vertical wind velocity is nonzero. This is achieved by incorporating knowledge about the scalar transport represented by the transport asymmetry coefficient. We show analytical and empirical ways to obtain the transport asymmetry coefficient and provide an interpretation of its value. Then, we describe the sensitivity of the flux to values 2 Theory

Eddy covariance
The net ecosystem exchange (NEE) of a scalar c (such as an atmospheric constituent), N c is the total vertical flux wc across the measurement plane at a height h and the change of storage below that height (Gu et al., 2012).
Where w is the vertical wind velocity (m s −1 ), c is the molar density (mol m −3 ) of the scalar of interest (such as CO 2 ). The previous equation can be reached either from a holistic mass balance approach or by averaging the continuity equation for the scalar c and integrating from the surface to measurement height h. In both cases, horizontal advection is ignored as a virtue of the assumption of horizontal homogeneity and molecular diffusion is ignored due to its small magnitude (Gu et al., 2012). For a full discussion on the equations of surface flux, see, for example (Finnigan et al., 2003;Foken et al., 2012a). 75 The storage term measurements and value are beyond the scope of this study, therefore we ignore it. Consequently, the total vertical flux is represented by the first term on the right hand side of Eq. (1) which can be further decomposed into turbulent and mean advective parts.
The overlines denote ensemble averages that obey Reynolds averaging rules. Primes represent departures from the mean.

80
The ensemble averages are estimated experimentally by the time averages, thus for a stationary time series drawn from an ensemble, the flux for the averaging period T avg can be written as where w(t) and c(t) are realizations of the vertical wind velocity and the scalar quantity such as CO 2 concentration, respectively.

85
The first term on the right-hand side of Eq.
(2) is the covariance between vertical wind velocity and the scalar concentration.
It is often referred to as the turbulent flux or the first-order approximation of the eddy flux, whereas the mean vertical advective term (second term on the right-hand side of Eq. (2)) is known as "Webb correction" or Webb-Pearman-Leuning (WPL) correction (Webb et al., 1980), and exists due to fluctuations in air density (Fuehrer and Friehe, 2002). The mean vertical wind velocity induced by density fluctuations of air parcels, hereafter w d , is often nonzero and is needed to account for the 90 mean advective term in the flux. It is, however, not possible to directly measure it. One reason is that w d caused by air density fluctuations is rather small (less than 1 mm s −1 ). But more importantly, any offset in the vertical wind velocity will appear in the measured w, consequently obscuring w d . Several reasons can contribute to a nonzero mean vertical wind velocity. This includes: tilted coordinates, biases in instruments, flow perturbations, and meteorological reasons induced by local circulation or topographical drainage Paw U et al., 2000;Heinesch et al., 2007). As a result, the measured slowly varying 95 termwc need to be discarded and the correct w d need to be estimated. One way to estimate w d is by utilizing the knowledge of the NEE of another scalar (Gu et al., 2012) or, in the case of WPL theory, by assuming that the net mean vertical mass flux of dry air is zero (stationarity of dry air) and calculating w d from sensible and latent heat fluxes.

True Eddy Accumulation
The true eddy accumulation method circumvents the need to measure individual realizations of the scalar concentration. In-100 stead, it is sufficient to measure the mean product wc for updraft and downdraft once at the end of each averaging interval T avg .
The product of w and c is realized by physically collecting air samples with a flow rate proportional to the vertical wind velocity w. The method is formulated assuming ideal conditions, where the mean vertical wind velocity during the averaging period is assumed to be zero. Whenw = 0, the second term on the right hand side of Eq.
(2) will be zero and the turbulent flux 105 w c will equal the mean product wc. By separating wc depending on the direction of the vertical wind velocity we can write Hence, by sampling air with a flow rate proportional to vertical wind velocity and accumulating it according to its direction in updraft and downdraft reservoirs, one can measure the quantity wc and consequently the flux without having to measure the 110 high-frequency fluctuations of the scalar, c (Desjardins, 1977;Hicks and McMillen, 1984).
A simpler alternative formulation can be reached using the law of total expectation. We write the flux as the expected value of the random variable wc conditional on the direction of the vertical wind velocity sign(w) Sampling air proportional to the magnitude of vertical wind velocity requires a scaling parameter, A that maps vertical wind 115 velocity to the flow rate. The scaling parameter is the product of the pump calibration coefficients and other coefficients used to adjust the system's dynamic range. For a short interval of time dt a sample of the volume V sample = A|w| dt will be collected in the system.
The accumulated sample volume in each of the two reservoirs during a long enough averaging period T avg (30 to 60 minutes) will be By the end of the averaging period T avg , the flux will be equal to the difference in the scalar accumulated mass between updraft and downdraft reservoirs.
If it is desired to formulate the flux in terms of the accumulated scalar concentration (mol m −3 ) instead of the accumulated mass, the average density of accumulated samples in each of the reservoirs will equal the accumulated mass of the scalar 125 divided by the accumulated volume Where C acc is the accumulated scalar density and the arrows indicate the reservoir. The measured concentration in Eq. (8) is the weighted mean of the scalar concentration and the magnitude of the vertical wind velocity. We can simply rewrite the accumulated concentration in terms of the wind and the scalar concentration as When w is assumed to be zero, |w ↑ | = |w ↓ | = |w|/2, and we can write the flux in terms of concentrations of accumulated samples, similar to Hicks and McMillen (1984)  to guarantee a zero mean vertical velocity. A common approach to nullify mean vertical wind velocity in EC measurements is to rotate the wind coordinates in post processing to force w to zero for each averaging interval, this method -commonly referred to as double rotation -is not feasible in eddy accumulation methods. The planar fit method (Wilczak et al., 2001) is better suited for the online application in the TEA method (Siebicke and Emad, 2019). The planar fit method aligns the sonic coordinates to the long-term streamline coordinates by aligning the wind vector to the plane that minimizes the sum of squares 145 of the vertical wind velocity means for a long period of time (weeks to months). This approach, while minimizes the vertical wind velocity means of the individual averaging intervals, does not force them to be zero. Considerable spread of w values around zero can still be observed after applying the planar fit method.
The key to extending the TEA equation to non-ideal case is to obtain an estimate of the scalar meanc, and consequently remove the advective termwc. We achieve this by using the weighted mean of c and |w| and correcting for the correlation 150 between them.
The weighted mean (c W ) of the scalar, c and wind magnitude |w| can be written similar to Eq. (9) as By decomposing into mean and fluctuating parts, we can write It follows that Substitutingc in Eq.
(2) we can write the flux as w c = wc −w |w| |w|c +w |w| |w |c We can obtain all the terms in Eq. (14) from our measurements except for the covariance term |w |c . 160 We define the "transport asymmetry coefficient" for the scalar c, (α c ) as the ratio of the covariance between the wind magnitude and the scalar to the covariance between the wind and the scalar.
where ρ c|w| , ρ cw are the correlation coefficients between c and |w|, c and w, respectively. σ |w| and σ w are the standard deviations of |w| and w, respectively. After substitution, we write the flux as where ρ is the correlation coefficient between the vertical wind velocity w and the scalar concentration c.
We can express the analytical value of α c in terms of the moments of the joint probability density function. We evaluate The term |w| c is obtained using After solving the integrals in Eq. (20) and Eq. (21) and substituting in Eq. (19), we find the value of α c can be written using the mean vertical wind velocity and standard deviation as where erf is the error function. Therefore, the error in the flux when when failing to account for the correlation between the scalar and the magnitude of vertical wind velocity will lead to a flux biased by the last term of Eq. (16)w/|w|α c . We can further substitute the expected value of |w| by the mean of the folded normal distribution (Leone et al., 1961) and obtain an analytical expression for the expectation of the flux error due to a nonzero vertical wind velocity 190 Experimental evidence has shown that different scalars behave similarly (Ohtaki, 1985;Wesely, 1988). We expect the value of α to be similar for different scalars due to the similar transfer mechanism. Therefore, an empirical approach for estimating α for one scalar is to use α value from another available scalar, e.g. sonic temperature.
We can express the value of α in terms of the updraft and downdraft contributions to the flux (flux ↑ and flux ↓) as the difference of updraft and downdraft fluxes to the total flux. 195 α c can be written to a good approximation as where the updraft and the downdraft contributions to the flux (flux ↑ , and flux ↓ ) are defined as

200
where P(w ↑ ) is the probability of a vertical wind velocity with a positive direction, P(w ↓ ) is the probability of a vertical wind velocity with a negative direction.
The expression of α in terms of updraft and downdraft flux contributions is closely related to quadrant analysis. We include it here for completeness. Quadrant analysis is commonly used to inspect the contributions from different quadrants in the (w , c ) plane by sorting the instantaneous values into four categories according to the sign of the two fluctuating components e.g. (Katul et al., 1997;Raupach, 1981;Katsouvas et al., 2007).

205
We find that α can be written in terms of quadrants as Where S i is the fraction of the flux transported by contributions in quadrant i, given as wc i is the conditional average defined as The indicator function I i obeys Following the definition of Thomas and Foken (2007), the pairs S 2 and S 4 are ejections and sweeps for downward directed 215 net flux (negative ρ wc ) and S 1 and S 3 for upward directed net flux (positive ρ wc )

Calculating the corrected TEA flux
The new general equation for TEA (Eq. 17) extends the validity of the method to conditions where the mean vertical wind velocity is nonzero. We show here how the correct TEA flux can be calculated from the measured physical quantities.
The weighted mean over an averaging period T avg can be written as which in terms of the quantities we are measuring, translates to Similarly After substitution and simplification, we obtain the correct TEA flux in terms of the measured quantities Where F TEA is the kinematic flux density (mol m s −1 ). C ↑ acc and C ↓ acc are the mean concentrations (mol m −3 ) of the scalar c in updraft and downdraft reservoirs at the end of the accumulation averaging period T avg . V ↑ and V ↓ are the accumulated sample volume (m 3 ) in updraft and downdraft reservoirs at the end of the averaging period. |w| is the mean of the absolute 230 vertical wind velocity (m s −1 ) during the averaging period.w is the mean of the vertical wind velocity. α c is the transport asymmetry coefficient for the scalar c (dimensionless).

Short-time eddy accumulation
The original formulation of the true eddy accumulation method requires the samples to be accumulated for the entire averaging interval T avg before the measurement can take place. This can pose limitations on the operation and the applicability of the 235 method.
We propose a modification for the TEA method, where samples can be accumulated for a sequence of shorter intervals τ i that add up to the averaging period T avg . This formulation can be achieved by applying the law of total expectation to the random variable cw with respect to a partitioning variable Y that divides the averaging period T avg into multiple non-overlapping partitions with the length τ i . It follows that the expectation of cw is the conditional expected value of cw given Y and the flux 240 is equal to This allows to write Eq. (8) as a sum of j intervals with the length of τ i and a scaling factor A i each. The concentration for either updraft or downdraft reservoirs can then be calculated from The concentration in either updraft or downdraft reservoirs at the end of the averaging interval is the mean of the short interval concentration measurements C i weighted by the sample volume during the short interval V i .
We call this new variety of eddy accumulation, the short-time eddy accumulation method (STEA).

Effect of buffer volumes 250
The short-time eddy accumulation method can be achieved in two ways, either using expandable buffer volumes (e.g. bags), which are emptied after each short interval measurement C i or using a flow-through system with rigid buffer volumes. The flow-through system has a practical operational benefit but requires additional correction to reverse the effect of buffer volumes on the concentration signal. Buffer volumes act as low pass filters (Cescatti et al., 2016). They attenuate the magnitude of the high-frequency part and shift the phase of the signal. The buffer concentration at time step n is dependent on the new input 255 sample concentration and the buffer concentration from the previous step y[n − 1]. Thus, the buffer volume concentration y n response to an input C i can be described with the difference equation whereq is a dimensionless flow rate that is defined as the ratio between the sample mass to the total mass of air in the buffer volume, at each time step n Where V i and ρ i are the volume and density of the short accumulation sample, respectively, while V b and ρ b are the volume and the air density of the buffer, respectively. Equation (38) characterizes a first-order linear filter. The treatment as a discretetime process aligns with the discrete operation of the STEA method.
The system response is characterized by the dimensionless flow rate. The time constant of the system is defined as the 265 required time for the system to reach 1/e from a step increase and relates toq by (Taylor et al., 2013).
where ∆t is the sampling interval.

Experiment period
The experiment started in September 2019. The first period (September 2019 to April 2020) was used for the development of the method which coincided with winter. The measurements during this period showed a lower quality and were excluded from 275 the analysis.
The winter period was followed by a period of stable operation, running from May 2020 until October 2020. During this period, we had a stable and continuous operation, however, intermittent discontinuities still existed due to blocked inlets, rain, or technical failures. From this period, we selected six weeks in summer, spanning from 18 June 2020 to 31 July 2020, to compare the different methods.

Instruments
EC and STEA measurement complexes were mounted at 5 m height above the ground (Fig. 1). The instruments used in the experiment for flux measurements and data analysis are listed in Table 1. Meteorological variables were logged using a Sutron 9210 XLite logger (Sterling, USA). All the raw data needed for flux processing were synchronized on the STEA computer and remote servers for real-time processing.

285
The EC system comprised a dedicated sonic anemometer (uSonic-3 Omni H) and an open-path infra-red gas analyzer (IRGA). Wind and scalar density data were acquired at 20 Hz frequency. Relative to the Class-A sonic anemometer used for STEA, the northward, eastward, and vertical separation of the IRGA was −17 cm, 26 cm, and −15 cm, respectively. The 2.7 TEA system description The TEA system used in the experiment is based on an earlier system of Siebicke and Emad (2019). The new system used the same mass flow controllers and shared most of the operating software. It has, however, several differences and improvements.
One major difference is the use of fixed stainless steel buffer volumes instead of expandable bags. The system was developed initially as a hybrid TEA-EC method to run the TEA method in a flow-through mode (Siebicke, 2016). The system was set up 295 to operate in the STEA flow-through mode described earlier in the theory section. A constant duration for the short intervals (τ i ) equal to one minute was used. The STEA system is comprised of two identical sampling lines, one for updrafts and one for downdrafts. Each of the sampling lines has two rigid buffers in a sequence connected using 6 mm Teflon tube (Fig. 2).
The STEA sampling inlets were installed in close proximity to the sonic's center of measurement volume. The horizontal separation was 22 cm, while the vertical separation between the two inlets was 2 cm. Upon sampling, the collected samples 300 were carried using 6 mm Teflon tubes to the first set of buffers. The sampling can be summarized in the following steps (see a detailed description of the system operation and sampling in (Siebicke and Emad, 2019)):  The colored bottom bar below shows the range of pressure values at each stage.
1. 3D wind measurements are acquired from the sonic anemometer (uSonic-3 Class A) with a 10 Hz sampling frequency.
2. Wind coordinates are rotated into the streamline coordinates using the planar fit method without an intercept (Dijk et al., 2004). The fit is performed online as a running window operation with a window width of 2 days and an update frequency 305 of once every 30 minutes.
3. The mean vertical wind from the previous window width (30 min.) is removed to minimize w. This is equivalent to applying a high-pass filter to the vertical wind velocity measurements. This section describes the steps followed to obtain the final and corrected STEA flux. Firstly, we discuss the effect of water vapor on the measured concentrations of other scalars and how we corrected that remaining water cross-sensitivity. Then, we present the procedure of data quality screening. Next, we detail the steps of calculating the final STEA flux. Finally, we present the buffer volume empirical correction we applied.
where r c is the dry mole fraction of the species c, χ c is the wet mole fraction measured by the instrument, and χ w is the water mole fraction measured by the instrument. For the LGR gas analyzer, these coefficients were estimated as a = 335 −1.219 × 10 −06 (±2.169 × 10 −09 ), b = 1.229 × 10 −12 (±1.073 × 10 −13 ) (Hiller et al., 2012) for CO 2 where the unit is ppm.
We found that using the same parameters could not control for all the effects of water. A linear slope different from zero was still found when supplying the gas analyzer with air of varying water signal and constant CO 2 . This suggests a remaining cross-sensitivity on the presence water vapor. To control for this small cross-sensitivity we used a linear fit to obtain the slope and corrected for it. 340 We were not able to supply the gas analyzer with air of known CO 2 signal. Instead, we used the systems buffer volume to collect air from the atmosphere, closed the inlets, and supplied the gas analyzer with enough sample flow rate for measurement.
The accumulated sample was enough to supply the gas analyzer for ca. 10 min. We repeated the measurements several times and used the obtained dataset for correcting the renaming cross-sensitivity.

345
Raw data were processed to ensure the removal of outliers due to measurement errors and instruments malfunction. This included the following steps Vickers and Mahrt (1997).
-Dropouts removal: some sensors would get stuck on one value, the first value is kept and the rest are discarded. -Deadband removal: measurement of short interval events involve regularly switching the sampling line coming to the gas analyzer from updraft to downdraft reservoirs. This will cause contamination from subsequent samples. We experimentally chose a 25-second deadband at the start of each short interval event. The measurements falling within the deadband 355 were removed. Figure 3 shows an example of deadband removal at the start of each averaging interval.
-Detection of sample contamination: periods where the flow rate to the gas analyzer is smaller than 400 mL min −1 are flagged. Under these conditions ambient air might enter the system and contaminate the collected samples. When the number of flagged data points exceeds 10% of the total points in the sampling interval, data in the sampling interval are discarded. After measurements are quality checked and erroneous data points are excluded, the final STEA flux is calculated as follows -Short interval statistics: for each short interval sample, the gas analyzer will have several repeated measurements for the concentrations C i , however, only one value is needed for the flux calculation. We use the median to obtain the representative value in order to minimize uncertainty and exclude outliers. Figure 3 shows an example of data quality 365 checking and choice.
-Calculate air molar volume: the molar volume of air is needed to express the flux in dynamic flux units. The molar volume is calculated using sonic temperature, pressure, and humidity measurements.
-Calculate short intervals weights: following Eq. (37), the measured short interval concentration should be weighted by the ratio of the accumulated volume during that interval to the total buffer volume V i /V tot .

370
-Calculate values of α θ : values of α θ are calculated using vertical wind velocity and sonic temperature measurements.
Values of α θ larger than 1 are discarded as they indicate a problem with the measurement.    We simulated the effect of buffer volumes on the high-frequency sonic temperature signal. The loss in the fluxes was parameterized by artificially degrading the sonic temperature in a procedure similar to Goulden et al. (1996) and Berger et al.

EC reference flux measurements and computations
The raw data from the two sonics and the high-frequency gas density measurements from the IRGA were used to compute eddy 385 covariance fluxes for water vapor and CO 2 in the period from 1 April 2020 to 1 November 2020 using EddyPro® software (LI-COR Env. Inc. USA) version 7.0.4. The flux processing steps were chosen to be as similar as possible to the TEA processing scheme. The calculation involved the following steps: -Statistical screening for the data quality issues following (Vickers and Mahrt, 1997), including tests for spikes, amplitude resolution, dropouts, absolute limits, and higher moment statistics.

390
-De-trending of raw time series by block averaging.
-Compensation of the time lag between the wind and the scalar time series using covariance maximization.
-Tilt correction using the planar fit method without an intercept (Dijk et al., 2004). The planar fit procedure was performed in two ways. First, for the entire experiment period, and second, as a running window operation with a width of 2 days, similar to the procedure performed online by the TEA system. -Analytical high and low-frequency corrections to correct for the spectral attenuation of the IRGA (Moncrieff et al., 2005(Moncrieff et al., , 1997

Data selection for method comparison
For comparing the fluxes calculated from both methods, we selected averaging intervals according to the following criteria: -Spike removal: following Vickers and Mahrt (1997) using a window width of 6 hours and a threshold of 2 standard 400 deviations.
-Rainy periods exclusion: data records where rain was recorded were excluded.
-Flux quality flags: periods where the flux quality flag is 1 or 2 according to Foken et al. (2005) were excluded.
-STEA low flow rate: averaging intervals flagged with the low flow rate flag described earlier were discarded.
After applying the above criteria, 1971 averaging intervals remained. They accounted for 54.1% of the whole comparison 405 period. Table 2 shows a summary of data quality checks results. To compare the overall difference between the two methods, we used the coefficient of determination R 2 and the slope of the orthogonal distance regression (ODR) (also known as major-axis regression and model II regression). ODR considers the errors in x and y as opposed to OLS regression which assumes the error in x is negligible (Wehr and Saleska, 2017).

Numerical simulations 410
In addition to the experiment described earlier, we set up a numerical simulation to test the magnitude of the error due to nonzerow on the flux. We used high-frequency measurements obtained from the IRGA and the sonic anemometer during one We applied minimal quality checks on the resulting fluxes before the comparison. We applied the steady state test following Foken et al. (2005), which removed 24% of the averaging intervals. We limited the values of |α θ | to less than 0.7, which 420 removed an additional 11% of the averaging intervals. The excluded averaging intervals occurred almost exclusively during low developed turbulence and night-time conditions.

Results and Discussion
We first discuss the problem of nonzero mean vertical wind velocity and the new generalized TEA equation. Then, we discuss the value and interpretation of the transport asymmetry coefficient α. Next, we discuss the newly proposed short-time eddy 425 accumulation method. Then, we discuss some results and aspects of the STEA flux calculations. Afterwards, we will present the flux intercomparison between STEA and EC. Finally, we discuss the effect of using fixed buffer volumes on the fluxes and the proposed empirical correction.

Nonzero mean vertical wind velocity
We presented a new formulation for the TEA equation. The new equation (Eq. (16)) employs information about the scalar 430 transport to allow the estimation ofc from available TEA measurements and, consequently, constraining the bias termwc.
Besides the correction of the nonzero w bias, the estimation of the scalar meanc is essential for the WPL correction and the calculation of storage fluxes.
The terms of Eq. 16 account for different contributions to the flux. The first term on the right hand side is equivalent to calculating the flux as the difference in accumulated mass between updraft and downdraft. Whenw = 0, the equation is 435 reduced to this term only. The second term accounts for the bias introduced by the advective termwc by using the weighted mean of the scalar and the magnitude of wind as an estimate forc. We show that whenw = 0, the first two terms are equivalent to using the concentration formula of Hicks and McMillen (1984) shown in Eq. (10), with the unequal volume correction of Turnipseed et al. (2009) that accounts for the small difference between the weighted meanc W and average of concentrations (C ↑ acc + C ↓ acc )/2. Refer to appendix A for details about this equality. The new third term c |w |/|w| corrects for the correlation 440 between the scalar and the magnitude of the wind. Ignoring the third term will result in a flux biased with the ratiow/|w|α c .
The new TEA equation reveals an important insight. When using the new equation to calculate the flux, the error in the flux when w = 0 is independent of the scalar concentration and is governed by the characteristics of the turbulent transport. This gives higher confidence in using the TEA method for measuring atmospheric constituents with high background concentration and small flux (low deposition velocity).

445
To quantify the effect of nonzero w on the fluxes, we used the results of the numerical simulation explained in the methods section. The results of the comparison are presented in Fig. 4. Additionally, we used the slope and the coefficient of determination, R 2 obtained from a linear fit of the calculated fluxes against the reference EC flux and the mean absolute difference, MAD to compare the different formulations.  . This assumption, although used in the literature e.g. (Wyngaard and Moeng, 1992), is not adequate. While the wind might be normally distributed for most stability classes (Chu et al., 1996), the scalar can depart significantly from normality (Berg and Stull, 2004). Other distributions might be more suited for approximating the joint probability distribution (Frenkiel and Klebanoff, 1973). For example, Katsouvas et al. (2007) showed using experimental data that a third-order Gram-Charlier distribution was necessary and sufficient in most of the cases for describing the quadrant time and flux contributions.
We found using high-frequency measurements that the value of α for CO 2 correlates moderately with the skewness of the measured scalar (r = 0.61). On average, updrafts have larger contribution to the flux. The mean of α for CO 2 and sonic temperature calculated from high-frequency measurements for periods with negligiblew was approximately 0.2 under unstable 470 and good turbulent mixing conditions (|ρ wc | > 0.25) with a standard error, SE = 0.01. Under stable stratification (ζ > 0), the mean of α was approximately equal to −0.18 but with a higher spread around the mean, SE = 0.09. These values generally agree with values found from studies using conditional sampling (Greenhut and Khalsa, 1982) and LES simulations (Wyngaard and Moeng, 1992) which found that updrafts contribution to the flux is 2 to 3 times larger than downdrafts.
We showed that the bias in TEA flux whenw = 0 is dependent on the value of α c . For TEA, α c is not readily available, since 475 its calculation requires the high-frequency information of the scalar. Similarity of scalar transport suggests that α values for different scalars should be similar, allowing the use of α θ calculated from sonic temperature as a substitute for α c . The similarity was confirmed empirically by calculating the values of α θ and α c from high-frequency measurements. A linear fit with a slope of 0.98 and R 2 of 0.962 was obtained during steady-state and well-developed turbulence conditions. During such conditions, α θ can substitute α c to calculate the flux correction ratio. However, the correction becomes large and unreliable in periods 480 where σ w and ρ cw are small, associated with small fluxes during night-time and stable conditions. Additionally, temperature is considered a bad proxy during near-neutral conditions (McBean, 1973;Hicks et al., 1980) due to its contribution to buoyancy.
We noticed that the variance in α values is higher under weakly developed turbulence. We experimentally determined the threshold for the safe use of α for correction as |ρ cw | = 0.2. Below this threshold, values of α larger than 0.5 are observed, making the correction unreliable. This threshold can be seen as an indicator for the violation of assumptions of homogeneity 485 and stationarity or other problematic conditions. Similar uses for the correlation coefficient are common in the literature e.g. (Foken and Wichura, 1996).
We determined experimentally that the error in the flux due to nonzerow becomes significant (larger than 10% of the flux) whenw exceeds 0.21σ w for periods with good turbulent mixing conditions (|ρ w,CO2 | > 0.2). This threshold is close to the analytical value of 0.323 σ w obtained from the Gaussian joint probability distribution. To push this threshold further, α θ 490 calculated from sonic temperature can be used during good turbulent mixing conditions (|ρ w,CO2 | > 0.2). Simulations indicate that the average relative confidence interval for the predicted value of α θ from α CO2 is 0.17% of the fit value. In summary, to keep the error in the flux below 10%, α θ can be safely used to correct for biasedw as long asw < 0.7 σ w . This limit is considered forgiving and easy to achieve with online coordinates rotation and other simple online treatments. The only times where this limit is expected to be reached is when σ w is very small (e.g. during night-time conditions) where other problems 495 such as low turbulence mixing and violations of the method's assumptions are expected to occur. These periods largely overlap with periods considered of low quality and are usually excluded from the analysis (Foken et al., 2012b).

Short-time Eddy Accumulation
We proposed the short-time eddy accumulation method (STEA), a modification of the TEA method where the accumulation time is divided into shorter intervals of time that add up to the flux averaging interval. The accumulation on shorter time 500 scales brings many advantages. First, it allows adapting to the local range of vertical wind velocity values which improves the resolution and dynamic range of the system. This can be achieved by exploiting the autocorrelation of the wind velocity signal to predict a scaling parameter A i better adapted to the local velocity field for each interval. For a short interval, the range that the sampling apparatus need to cover will be on average smaller than the range of the whole averaging interval. We found for a short averaging interval of one minute, the range was on average 60% smaller than that of the whole flux averaging interval.

505
As a result, the upper bound of the required dynamic range for w reported by Hicks and McMillen (1984) as 5 σ w is lowered to 3.33 σ w .
Additionally, the accumulation on varying intervals means the measurement frequency can be adjusted to match that of the gas analyzer or the precision requirements. This can be useful for reactive species and other traces gases, where relatively fast gas analyzers are available but not fast enough for EC. Finally, the STEA method facilitates using the STEA system in flow-through mode using rigid reservoirs. The operation in flow-through mode requires two sets of buffer volumes in a series. The ideal operation of such a system can be achieved as  It is important for this scheme to keep the mass flow rate to the second set of buffers constant so that the assumption of time 525 invariance of the linear filter used to model the buffers is not violated.

STEA fluxes computations
In this section, we will discuss some aspects related to the calculation of the STEA fluxes. We first discuss the effects of water vapor on CO 2 concentration measurements. Then, we discuss the effect of coordinates rotation on the fluxes. Finally, we discuss the effect of density fluctuations on eddy accumulation methods.

Water correction
Treatment of the residual cross-sensitivity of CO 2 on water signal using a linear fit produced a small slope of −1.17 × 10 −4 shown in Fig. (7). Thus, a difference in water concentration of 4000 ppm between updraft and downdraft reservoirs, typically observed at extreme conditions, will lead to a difference on the order of 0.5 ppm for CO 2 .
Applying the water correction using the quadratic fit and the slope correction reduced the STEA fluxes in comparison to the 535 direct calculation of mixing ratios. However, it improved the fit between the STEA and the reference EC flux ( Correction none polynomial slope Figure 7. Effect of water correction on the measured CO2 concentration using the LGR FGGA-24r-EP instrument. Points represent measured CO2 by the gas analyzer when air with constant CO2 concentration and varying H2O concentration was supplied. Lines represent linear regression fits. Red colored points and line represent CO2 measurements after applying the polynomial correction (Hiller et al., 2012;Rella, 2010). In blue are the CO2 measurements after applying our slope adjustment correction to remove additional cross-sensitivity on water.

Coordinates rotation
The online coordinates rotation produced stable rotation angles over the experiment period. The eddy covariance fluxes calculated using the Class-A sonic using a two-month long dataset (1 June 2020 to 1 August 2020) produced the rotation angles: The use of online rotation with a moving window of two days minimized the residual vertical mean wind in comparison to using the whole period of the experiment. This is likely due to a better adaptation to the local wind field. The mean vertical wind 545 values ranged from −0.05 σ w to 0.38 σ w (mean magnitude of 0.251 σ w ) compared to −0.04 σ w to 0.04 σ w (mean magnitude of 0.06 σ w ) when using the online rotation with a short moving window.
To estimate the effect of the online rotation method on the fluxes, we calculated EC fluxes using the two different rotation approaches while keeping other treatments constant. The comparison revealed that the online rotation with a moving window had minimal effect on the fluxes: a slope of 1 and R 2 of 0.98 were obtained when using a linear fit. Nevertheless, this compar-550 ison only included data of good quality. A comprehensive comparison might be needed to identify the effect under non-ideal conditions.

Effect of density fluctuations
Changes in air density due to temperature, pressure, and dilution of water vapor and other gases bias the measured flux. If the density of a scalar is measured, a correction is needed to account for these effects (Webb et al., 1980). However, If the mixing 555 ratio is measured, no correction is required since it is a conserved quantity. In TEA and STEA, after samples are collected and mixed in buffer volumes, the mean mixing ratio is measured. Therefore, no correction for density effects is needed. The resulting flux is equivalent to the flux measured with mixing ratios r c w . However, the density fluctuations might affect TEA and STEA differently. Since, the mass flow rate of air is dependent on air density. The more dense the air is, the higher the mass flow rate is. If such an effect is not taken into account by using mass flow sampling, the resulting flux will be biased. The 560 bias is equivalent to having the wind speed measurement dependent on air density in EC.

STEA/EC flux intercomparison
The measured CO 2 fluxes using the STEA method in flow-through mode showed a good match with the reference EC fluxes ( Fig. (8)).
The time series of measured CO 2 fluxes in Fig. (8 -a) shows that the STEA method was able to reproduce the daily dynamics 565 of CO 2 flux. The estimated fluxes using the STEA method appear to have less spikes and smoother in general, this is likely due to the smoothing effect of buffer volumes and the lower sensitivity of the closed path gas analyzer to rain and high humidity.
The mean diurnal cycle estimates from the two methods shown in Fig.(8 -b) match very well. However, a small time shift can be observed on the mean diurnal cycle a result of the phase shift introduced by the low-pass filtering effect of the buffer volumes.

570
The regression results shown in Fig. (8 -c) show that the measured CO 2 fluxes using the STEA method in flow-through mode have a very good agreement with the reference EC fluxes. The magnitude of STEA fluxes was comparable to EC fluxes (ODR slope = 1.05). This indicates that the STEA method does not introduce systematic error to the fluxes. However, the coefficient of determination R 2 was 0.87, which indicates a 13% unexplained variance contributed by the uncertainty in the two estimates. We suggest three different mechanisms contributing to the observed uncertainty leading to the unexplained 575 variance. First, the random sampling error arising from the stochasticity of turbulence (Hollinger and Richardson, 2005). The observed uncertainty from the two methods calculated as the standard deviation of the difference is 4.29 µmol m −2 s −1 this is comparable to reported value of 2.7 µmol m −2 s −1 using two tower estimates (Hollinger and Richardson, 2005). The errors also show heteroscedasticity with the error increasing along with the absolute magnitude of the flux, a similar behavior was observed by Hollinger and Richardson (2005)  Finally, the different processing steps between the two methods can contribute to the uncertainty. In particular, the effects of time-lag compensation, spectral corrections, and statistical screening. We determined the combined effect of these treatments by calculating the EC flux with and without the treatments and found that the effect on the flux was negligible.

Effect of buffer volumes
Using fixed buffer volumes attenuates the signal. The effect of buffer volumes can be described as a low-pass first-order linear 590 filter. Figure (9) shows the filter's magnitude and phase responses. The magnitude response |H| plot shows how the magnitudes of different frequencies are attenuated, the smaller the dimensionless flow rate is, the larger the time constant is. Consequently, the attenuation is stronger. To understand the effect of buffer volumes use on the measured scalar concentration. We carried out a simulation on a surrogate signal generated from sonic temperature. The simulation showed the decline can reach up to 10% of the fluxes under 595 operation ranges similar to those of our experiment (for τ = 11 minutes) (Fig. 10). The empirical correction was consistently able to mitigate most of the attenuation when the filter properties are assumed to be constant, (i.e the flow rate needs to be constant for each short interval). This assumption was difficult to maintain using the 1-minute switching regime. The simulation showed the empirical correction for the buffer volumes worked best when the correction factor was obtained using a linear fit, as opposed to taking a ratio of the attenuated flux to the true flux for each averaging interval. The correction factor, 600 in this case, is the reciprocal of the slope of the linear regression between the attenuated flux and the true flux. The correction factor calculated using Eq. (10) shows a good agreement between sensible heat flux and CO 2 . However, the uncertainty of the correction factor increased with increasing buffer volume time constant. For our experiment, the average time constant for the first-order linear filter used to model the buffer volume was estimated to be τ = 700 seconds. This value was used to simulate the loss on the fluxes using the sensible heat flux calculated from the sonic anemometer. The correction factor was obtained 605 from the slope of the attenuated flux and was equal to 1.18 In this paper, we revised the theory of the true eddy accumulation method and extended the applicability of the method to non-ideal conditions where the mean vertical wind velocity during the averaging interval is not zero. The new generalized equation allows estimating the scalar mean during the flux averaging interval and define conditions where the error in the 610 flux is significant. We found that the error in the TEA flux is a function of the disparity of atmospheric transport and defined ways and conditions to constrain that error. We found that it is possible to achieve minimum bias in the TEA flux under most atmospheric conditions. We believe these results increase the confidence in using the TEA method for different atmospheric constituents and under a variety of atmospheric conditions.
Additionally, we proposed an alternative method for the measurement of ecosystem-level fluxes. The new method, referred to 615 as short-time eddy accumulation (STEA), allows the sample accumulation to be carried out on shorter varying-length intervals.
The method offers more flexibility than TEA and has many potential benefits including a better dynamic range and higher accuracy than the TEA method, and the ability to operate under a flow-through scheme using fixed buffer volumes. The flexibility introduced by the new method offers new ways to design eddy accumulation systems particularly suited for a specific atmospheric constituent or a specific gas analyzer. For example, the accumulation time can be tailored to measure reactive 620 species or to distribute the gas analyzer time to measure the fluxes at different heights.
Furthermore, we presented a prototype evaluation of the STEA method under the flow-through regime. We described the details of the system design and operation. We compared flux measurements from our new system against a reference EC system over a flat agricultural field. The fluxes from the two methods were in very good agreement. We highlighted the importance of different processing and design aspects between the two methods and their potential effect on the fluxes.

625
Finally, we analyzed the effect of buffer volumes in the flow-through operational mode on the fluxes and proposed an empirical correction to correct for the underestimation resulting from the low-pass filtering behavior of the buffer volumes.
In summary, the generalized TEA equation and the new STEA method provide direct flux measurement methods that complement the state-of-the-art EC method. They extend the coverage of micrometeorological methods to new trace gases and atmospheric constituents beyond the scope of the EC method.
is equivalent to using (C ↑ acc + C ↓ acc )/2 as an estimate forc in the second term on the right hand side of Eq. (2). We write the conditional expectation ofw as w = w|sign(w) = |w ↑ | P(w ↑ ) − |w ↓ | P(w ↓ ) where sign(w) is the sign of vertical wind velocity. P(w ↑ ) and P(w ↓ ) are the observed probabilities of the sign of w, which equals the ratio of the time the wind is positive or negative to the total integration interval time.
and similarly by substituting |w|/2 with After rearrangement and simplification we get to When Eq. (A6) is compared with Eq. (2), it is clear that the term C + acc +C − acc 2 is used as an estimate forc.