True eddy accumulation trace gas flux measurements: proof of concept

Micrometeorological methods to quantify fluxes of atmospheric constituents are key to understanding and managing the impact of land surface sources and sinks on air quality and atmospheric composition.

Important greenhouse gases are water vapor, carbon dioxide, methane, and nitrous oxide. Further important atmospheric constituents are aerosols, which impact air quality and cloud formation, and volatile organic compounds. Many atmospheric constituents therefore critically affect the health of ecosystems and humans, as well as climate.

The micrometeorological eddy covariance (EC) method has evolved as the method of choice for CO₂ and water vapor flux measurements using fast-response gas analyzers. While the EC method has also been used to measure other atmospheric constituents including methane, nitrous oxide, and ozone, the often relatively small fluxes of these constituents over ecosystems are much more challenging to measure using eddy covariance than CO₂ and water vapor fluxes. For many further atmospheric constituents, eddy covariance is not an option due to the lack of sufficiently accurate and fast-response gas analyzers.

Therefore, alternative flux measurement methods are required for the observation of atmospheric constituent fluxes for which no fast-response gas analyzers exist or which require more accurate measurements. True eddy accumulation (TEA) is a direct flux measurement technique capable of using slow-response gas analyzers. Unlike its more frequently used derivative, known as the relaxed eddy accumulation (REA) method, TEA does not require the use of proxies and is therefore superior to the indirect REA method.

The true eddy accumulation method is by design ideally suited for measuring a wide range of trace gases and other conserved constituents transported with the air. This is because TEA obtains whole air samples and is, in combination with constituent-specific fast or slow analyzers, a universal method for conserved scalars.

Despite the recognized value of the method, true eddy accumulation flux measurements remain very challenging to perform as they require fast and dynamic modulation of the air sampling mass flow rate proportional to the magnitude of the instantaneous vertical wind velocity. Appropriate techniques for dynamic mass flow control have long been unavailable, preventing the unlocking of the TEA method's potential for more than 40 years.

Recently, a new dynamic and accurate mass flow controller which can resolve turbulence at a frequency of 10 Hz and higher has been developed by the first author. This study presents the proof of concept that practical true eddy accumulation trace gas flux measurements are possible today using dynamic mass flow control, advanced real-time processing of wind measurements, and fully automatic gas handling.

We describe setup and methods of the TEA and EC reference flux measurements. The experiment was conducted over grassland and comprised 7 d of continuous flux measurements at 30 min flux integration intervals. The results show that fluxes obtained by TEA compared favorably to EC reference flux measurements, with coefficients of determination of up to 86 % and a slope of 0.98.

We present a quantitative analysis of uncertainties of the mass flow control system, the gas analyzer, and gas handling system and their impact on trace gas flux uncertainty, the impact of different approaches to coordinate rotation, and uncertainties of vertical wind velocity measurements.

Challenges of TEA are highlighted and solutions presented. The current results are put into the context of previous works. Finally, based on the current successful proof of concept, we suggest specific improvements towards long-term and reliable true eddy accumulation flux measurements.

1 Introduction

The ability to observe the exchange of trace gases between the earth's surface and the atmosphere is key to understanding the functioning of ecosystems. Trace gas flux measurements allow quantification of how natural and anthropogenic systems affect atmospheric composition.

Many studies over the past decades have observed carbon dioxide (CO₂) and water vapor fluxes at ecosystem scale using micrometeorological methods (Baldocchi et al., 1988). Eddy covariance (EC) (Baldocchi, 2003, 2014) has become the most widely used method for measuring turbulent fluxes. Today the EC method is routinely being applied the world over including the major flux networks FLUXNET, ICOS, and NEON.

The EC method requires fast-response gas analyzers which only exist for a few trace gas species, above all CO₂ and water vapor but more recently also other trace gases, including methane (CH₄) and nitrous oxide (N₂O). However, for a large number of trace gases and atmospheric constituents, the applicability of the EC method is limited by a lack of fast-response gas analyzers, by the high power demand necessary for sustaining high sample flow rates in some closed-path gas analyzer systems, and by a possibly small signal-to-noise ratio of high-frequency measurements.

A number of alternative turbulent flux measurement methods exist which can use slow-response gas analyzers and might provide more accurate results than eddy covariance with fast-response analyzers. These methods are applicable to a wide range of conserved trace gases, isotopes, aerosols, volatile organic compounds, and other atmospheric constituents. An overview on selected micrometeorological methods applicable to slow-response gas analyzers follows, presenting the air sampling principles and timings and stating advantages and disadvantages of each method.

True eddy accumulation (TEA) is an alternative to the EC method. Unique properties of the TEA method are highlighted which make TEA stand out from other methods. This study is a contribution towards a practical implementation of the TEA method.

1.1 Micrometeorological methods suitable for slow-response gas analyzers

1.1.1 True eddy accumulation (TEA)

True eddy accumulation (Desjardins, 1977; Hicks and McMillen, 1984) refers to the sampling of air, separating updrafts and downdrafts on the condition of the sign of the vertical wind velocity. The mass flow rate of physical air samples needs to be proportional to the magnitude of the vertical wind velocity and controlled at 10 Hz or above to resolve flux-relevant turbulence scales. For conserved scalars, the net flux can then be determined from the difference in scalar concentration between the accumulated updraft and downdraft samples, respectively, over a certain flux integration interval, e.g., 30 min.

The idea of eddy accumulation (EA) goes back to early considerations by Desjardins (1972) who proposed the method for physically sampling trace gas fluxes. He reported a first experiment of conditionally sampling temperature and deriving sensible heat flux through mathematical accumulation (Desjardins, 1977). We use the term “true eddy accumulation” rather than just “eddy accumulation” to refer to the original formulation of eddy accumulation (Desjardins, 1977; Hicks and McMillen, 1984), specifically with vertical wind proportional air sampling, as opposed to later derivatives of eddy accumulation such as “relaxed eddy accumulation”, which is subject to constant mass flow and further limitations (see Sect. 1.1.2).

Literature on true eddy accumulation is sparse, with just over a dozen published studies. Very few studies performed actual flux measurements. Desjardins (1977), Speer et al. (1985), Neumann et al. (1989), Beier (1991), and Komori et al. (2004) presented early prototypes of true eddy accumulators and disjunct true eddy accumulators (Rinne et al., 2000). Others conducted simulations (Hicks and McMillen, 1984; Businger and Oncley, 1990) and contributed technology (Buckley et al., 1988) and reviews (Businger, 1986; Speer et al., 1986; Hicks et al., 1986). However, its practical implementation has long been difficult, particularly the accurate and dynamic control of mass flow rates. None of the experiments above produced significant long-term data sets. Correlation of TEA fluxes with EC fluxes was generally relatively low with coefficients of determination of, e.g., R²=0.07 (Speer et al., 1985), R²=0.41 (Neumann et al., 1989), and R²=0.64 (Komori et al., 2004). Until today there has been no TEA instrument commercially available.

Recently we have successfully performed a series of TEA flux experiments using a new and fully digital approach to dynamic and fast mass flow control and real-time processing of wind data. We are further working to advance TEA flux corrections and TEA simulations. Those experiments (unpublished) yielded a tight correlation between TEA and EC flux measurements, with coefficients of determination of up to R²=0.96, exceeding R² values from any of the above cited literature. The current work presents the first of the TEA and EC intercomparison experiments performed over short vegetation during spring 2015 in detail.

Figure 1 True eddy accumulation (TEA; see Sect. 1.1.1) and disjunct eddy accumulation (DEA; see Sect. 1.1.3) sampling scheme. Vertical wind, w, (a), scalar density, CO₂, (b) and vertical wind proportional mass flow rate (c). Black solid lines indicate the continuous true atmospheric signal. Sampling resolution of TEA and DEA are 10 Hz and 10 s, respectively. Note that active sampling time for DEA is only 1 % of the sampling time for TEA.

The concept of the TEA sampling scheme is illustrated in Fig. 1. The true vertical wind velocity (Fig. 1a, black line) is sampled at a frequency of, e.g., 10 Hz (Fig. 1a, blue dots), using an ultrasonic anemometer. Likewise, the air sampling device samples the true atmospheric time series of the scalar, e.g., CO₂ (Fig. 1b, black line), at the same time resolution of 10 Hz (Fig. 1b, blue dots). The time variable flow rates at which samples are being accumulated are shown in Fig. 1c. Separate accumulation of updrafts (red lines) and downdrafts (orange lines) is distinguished.

Following whole air sampling, the atmospheric constituent of interest can be trapped in a number of ways. Constituents can be accumulated as whole air samples in bags, absorbed in gas washing reservoirs, adsorbed on to chemicals using cartridges, continuously sampled with denuders, trapped as reaction products with chemicals, or retained using mechanical filters.

The true eddy accumulation principle is not limited to passive trace gases. Here, we suggest that the TEA method has the potential to measure fluxes of dust, pollen, bacteria, fungi, and other biological material carrying physical, chemical, and genetic information. The latter materials can be accumulated on appropriate filter media.

True eddy accumulation has a number of advantages over other methods. Sample accumulation over the duration of typical flux averaging intervals of 30 to 60 min allows for the use of slow-response gas analyzers. The key advantage of TEA over EC is the applicability to a much wider range of atmospheric constituents, assuming that slow-response analyzers are more readily available than fast-response analyzers, and better accuracy can be obtained through signal averaging.

The key advantage of TEA over other variants of eddy accumulation, i.e., relaxed eddy accumulation or hyperbolic relaxed eddy accumulation, is that true eddy accumulation is the only direct method in the family of accumulation methods. As a direct method it does not require the use of proxies (other scalars) and coefficients like the β coefficient in relaxed eddy accumulation and therefore does not depend on scalar similarity (Ruppert et al., 2006). This property of a direct measurement method is essential for quantifying fluxes of constituents which cannot be measured by other means (e.g., the EC method). Scalar similarity of the fluxes of the constituent of interest and the proxy cannot be assessed without first quantifying both fluxes themselves. The direct TEA method is independent of prior knowledge.

Another advantage over other types of eddy accumulation (relaxed eddy accumulation or hyperbolic relaxed eddy accumulation) or any type of disjunct eddy sampling (e.g., the disjunct eddy covariance method or the disjunct eddy accumulation method) is the continuous sampling of the air by the TEA method such that the signal is recovered in its entirety. Continuous sampling avoids noise associated with disjunct sampling (Lenschow et al., 1994). Likewise, omitting samples at times of small vertical wind velocities, which is common practice in relaxed eddy accumulation, would effectively be disjunct sampling, trading in noise for the sake of higher concentration differences between accumulated updrafts and downdrafts.

The long averaging intervals further allow for repeated measurements of the same sample, improving precision. The trace gas concentration of the accumulated samples, which is by design constant, at the time of analysis and the typically long analysis integration times are best matched with low sample flow rates through the gas analyzer. Low flow rates result in low power consumption and a low pressure drop over system components. A low pressure drop is beneficial for the stability and accuracy of the gas analyzer's reading.

1.1.2 Relaxed eddy accumulation (REA)

Given the challenges associated with the original formulation of true eddy accumulation, Businger and Oncley (1990) proposed a modified version of eddy accumulation, today known as relaxed eddy accumulation (REA). REA is based on the concept of flux–variance similarity. In order to relate the scalar flux to the variance of the vertical wind velocity, a proportionality factor, β, was introduced, so REA became an indirect method.

Figure 2 Relaxed eddy accumulation sampling scheme. Vertical wind, w (a), scalar density, CO₂ (b), and mass flow rate (c). Black solid lines indicate the continuous true atmospheric signal. Sampling resolution of REA is 10 Hz. A fraction of the CO₂ time series (gray lines, b) is not sampled by REA due to use of a vertical wind velocity dead band for small velocities (thresholds indicated by red lines, a). This dead band causes gaps in the otherwise constant mass flow rate (c).

The advantage of the relaxed eddy accumulation method is that air is sampled at a constant flow rate (Fig. 2c). This meant that the dynamic high-frequency modulation of flow rates as a function of the magnitude of the vertical wind velocity as in the TEA method was no longer required. REA still accumulates updrafts and downdraft separately controlled by the sign of the vertical wind velocity.

A second modification was introduced in REA: at times of small positive or negative vertical wind velocities, no air samples are taken. This “dead band” is illustrated in Fig. 2a. Figure 2b shows the air sampling scheme: the true scalar time series, e.g., CO₂ density (black line), is sampled at a regular frequency of, e.g., 10 Hz (blue dots), if the vertical wind velocity (Fig. 2a, black line), sampled at the same 10 Hz frequency (Fig. 2a, blue dots), is larger than the thresholds defining the dead band. A certain fraction of the scalar time series is thus omitted from sampling (Fig. 2b, gray line).

The use of a dead band has two advantages: the concentration difference between the updraft and downdraft accumulated samples increases (Pattey et al., 1993; Katul et al., 1996), improving the ratio of the flux signal to the noise of the gas analyzer. Secondly, use of a dead band leads to less frequent switching between updraft and downdraft samples, which relaxes the need for fast-response valves to some degree and would reduce material wear. However, a lack of scalar similarity can lead to flux underestimation as simulated by Ruppert et al. (2002), who also found that flux errors increased with dead-band size. Another disadvantage is the impact of the dead band on the flux itself of unknown magnitude, depending on the co-spectrum of scalar and vertical wind velocity.

The simplifications of the REA method relative to the TEA method, particularly the constant mass flow rate, have facilitated wide adoption of the REA method. More than 200 studies on REA flux measurements and simulations have been published since its description less than 30 years ago (Businger and Oncley, 1990). The significant number of REA studies suggests that there is a need for alternatives to the eddy covariance method for certain applications.

Despite being simpler to implement than TEA, REA has distinct disadvantages. Being an indirect method, the accuracy of REA remains critically dependent on the correct determination of an a priori unknown β factor. β varies with scalar and with atmospheric conditions. Typical β values obtained from measurements and simulations (Wyngaard and Moeng, 1992; Businger and Oncley, 1990; Oncley et al., 1993; Pattey et al., 1993; Baker et al., 1992; Gao, 1995; Milne et al., 1999; Katul et al., 1996; Baker, 2000; Ammann and Meixner, 2002; Held et al., 2008) are around 0.55 but range from ca. 0.4 to ca. 0.7, introducing significant uncertainty of up to several tens of percent of the measured flux.

Scalar similarity between a constituent of interest and a suitable proxy for determination of the β factor is often lacking (Ruppert et al., 2006; Cancelli et al., 2015). The alternative use of a constant β factor leads to lower accuracy of the estimated flux (Foken and Napo, 2008; Ruppert et al., 2002).

A variant of REA is hyperbolic relaxed eddy accumulation (HREA) (Shaw, 1985; Bowling et al., 1999). HREA maximizes concentration differences between accumulation reservoirs through the use of hyperbolic dead bands. Thus, HREA can resolve small fluxes such as stable isotope fluxes of ¹³C and ¹⁸O (Bowling et al., 1999; Wichura et al., 2000). However, HREA requires proxies similar to REA and omits about two-thirds of total sampling time through the use of dead bands. Dead bands can increase flux uncertainty by omitting parts of the signal due to incomplete sampling of the time series.

1.1.3 Disjunct eddy accumulation (DEA) and disjunct eddy covariance (DEC)

Disjunct eddy sampling (Rinne et al., 2000; Turnipseed et al., 2009) is based on considerations by Lenschow et al. (1994) for representing turbulent time series by temporal subsamples. Disjunct eddy covariance (DEC) takes very short grab samples (ca. 0.1 s), followed by a pause (e.g., 5 to 60 s) for gas analysis with relatively slow instruments. Similarly, disjunct eddy accumulation can be used to obtain short grab samples at a mass flow rate proportional to the magnitude of vertical wind velocity when continuous dynamic mass flow control can not be performed.

Figure 3 Eddy covariance and disjunct eddy covariance sampling scheme. Vertical wind, w (a), scalar density, CO₂ (b). Black solid lines indicate the continuous true atmospheric signal. Sampling resolution of EC and DEC are 10 Hz and 10 s, respectively. Note that active sampling time for DEC is only 1 % of the sampling time for EC.

The disjunct sampling principle is illustrated in Fig. 1 for the DEA method and in Fig. 3 for the DEC method. Comparing the few disjunct samplings at a resolution of 10 s of the DEA method (flow rate indicated by black vertical lines in Fig. 1c at times of 0, 10, and 20 s) relative to the continuous flow rate of the TEA method (red and orange vertical lines in Fig. 1c) illustrates the small fraction of the total time series being actually sampled by DEA.

Disjunct sampling allows more time for the analysis of the chemical species than continuous sampling. However, the uncertainty of disjunctly sampled scalar and wind time series, and as a result the flux uncertainty, is larger compared to continuous sampling. Turnipseed et al. (2009) found an additional uncertainty of ±30 % of the flux due to disjunct sampling and estimated the overall uncertainty of their DEA flux measurements as ±40 %.

1.1.4 Challenges of eddy accumulation

There are a number of challenges associated with eddy accumulation flux measurements (see also Hicks and McMillen, 1984). The first two listed below are specific to the TEA and DEA methods. The others are common to all eddy accumulation methods.

1. Mass flow control. The air sampling, i.e., the separation of updrafts and downdrafts as well as the response of the vertical wind velocity proportional mass flow control, needs to be sufficiently fast (10 to 20 Hz) and dynamic to resolve the relevant turbulent fluctuations. Further, the mass flow control needs to be accurate, even under dynamically changing flow rate conditions despite the compressibility of air. Finally, the dynamic range of the mass flow control, i.e., the ratio of the largest to the smallest accurately controllable mass flow, needs to be on the order of 100 or higher to limit flux errors (Hicks and McMillen, 1984). No commercially or otherwise readily available technology for fast, dynamic, and accurate control of mass flow rates exists or has been demonstrated to perform well in TEA.

2. Density fluctuation effects. Density fluctuations due to heat and water vapor transfer affect the flux of the scalar of interest. Corresponding corrections specific to the TEA method have been unavailable.

3. Spectra and co-spectra. No turbulence spectrum of the scalar nor the co-spectrum of the scalar and the vertical wind velocity can be obtained from the accumulated samples as they are mixed, and time-resolved analysis is therefore not possible. Spectral information on the wind is of course available as in any other method using a fast-response ultrasonic anemometer.

4. Coordinate rotation (see further details in Sect. 2.7). Sampling decisions need to be performed in real time, and they are definitive; i.e., they cannot be modified in post-processing. This is an important difference to the EC method. The separation of updrafts and downdrafts depends on the definition of vertical wind velocity, w. The mean w over the averaging period needs to be zero (see Sect. 2, Eq. 6). One way to minimize mean w is to align the coordinate system of the wind measurements to the mean streamlines through coordinate rotation (Wilczak et al., 2001).

   In eddy accumulation, the coordinate rotation needs to be performed in preprocessing to be available in real time. Coordinate rotation and any other operation attempting to nullify mean w over the flux averaging interval would require knowledge of w over the entire interval, including future observations. However, only past and present data are available in real time to approximate coordinate rotation and perform sampling decisions based on the sign and magnitude of w. Remaining nonzero mean w causes flux bias.

5. Decorrelation through sensor separation. Spatial separation of the gas sampling inlets and the wind sensing volume causes a time lag between the wind measurement and the time for obtaining the corresponding air sample. Not accounting for time lags leads to decorrelation of wind and scalar and therefore flux loss. Contrary to the EC method, where time lags can be detected through covariance maximization and corrected for, in eddy accumulation such post-processing is not possible because high-frequency scalar time series are not obtained. Therefore the wind measurement and the air inlet need to be co-located as close as possible.

6. Analyzer sensitivity. Trace gas concentration differences between reservoirs might be too small to be resolved by given gas analyzers.

7. Reliability. Eddy accumulation systems are mechanically and electronically complex machines. Particularly moving parts pose the risk of failure. Careful design is required for robust implementations and unattended long-term deployments.

We address the above challenges in the following ways:

1. A new type of digital and highly dynamic mass flow controller was deployed. The technology previously developed by the first author has a fast and dynamic response sufficient to resolve relevant turbulent scales at 10 Hz and above. The design accounts for the compressibility of air in dynamic sampling. See Sect. 3.2 for further details.

2. TEA-specific adaptation of flux corrections accounting for density fluctuations as proposed for the EC method by Webb et al. (1980) is subject to ongoing research by the authors.

3. Here we propose that, in absence of co-spectral information which is required to correct for flux attenuation due to sonic anemometer path length averaging and due to the separation of the air inlet and the sonic sensing volume, further research should investigate whether a proxy scalar such as air temperature or sonic temperature might be used to estimate the attenuation. While this would imply scalar similarity of the proxy and the scalar of interest, the impact of potential non-similarity would be limited to the spectral attenuation flux correction term, which is typically small relative to the flux itself, and the uncertainty from the scalar similarity constraint would be even smaller.

4. Nonzero mean vertical wind velocities were minimized through real-time coordinate rotation with continuous near-real-time updates of the rotation coefficients as well as a procedure to minimize remaining bias in the mean vertical wind velocity (see Sect. 2). Further, a correction accounting for volume mismatch between updraft and downdraft reservoirs due to nonzero mean vertical wind velocities (Turnipseed et al., 2009) has been applied (see Sect. 2).

5. Spatial separation of the wind sensing volume and the point of air sampling can be minimized through integration of the air inlet inside the wind sensing volume of the sonic anemometer. A certain degree of time lag between the wind signal and the scalar sampling will ultimately remain as long as the wind is sensed over a measurement volume rather than at a point and with discrete finite time resolution rather than truly instantly.

6. Performance analysis of a typical infrared gas analyzer model for the measurement of CO₂ gave satisfactory results in terms of the resolution but revealed limited stability. Subsequent work by the authors with current laser spectrometers gave superior results due to their improved stability. Details on the latter work will be reported separately.

7. Suggestions towards a robust design of an eddy accumulator are given in the conclusions of the current study.

1.2 Objectives

Out of all the methods discussed above, the true eddy accumulation method is the only alternative to the eddy covariance method for directly measuring the physical flux. Every effort should be made towards mastering the dynamic mass flow control necessary for direct TEA as well as addressing the other challenges listed above.

It is the objective of this work to deliver the proof of concept of true eddy accumulation trace gas flux measurements based on dynamic and accurate mass flow control proportional to vertical wind velocity and based on fully digital and real-time signal processing.

2 Materials and methods

2.1 Experimental design

The experimental design comprises three elements: novel true eddy accumulation flux measurements of CO₂ and water vapor, conventional eddy covariance flux measurements of CO₂ and water vapor, and meteorological measurements. All measurements were performed side by side. This allowed for the performance of the new TEA system to be evaluated on its ability to measure turbulent fluxes of CO₂ by relating the observed fluxes to meteorological drivers and by comparing the TEA CO₂ flux measurements to conventional EC CO₂ flux measurements. This section provides details on the measurement site, the meteorological measurements, the TEA method, technical implementation, and flux computation as well as information on the reference flux measurements by EC.

2.2 Experimental site

Flux measurements were performed at a grassland experimental field site of the University of Göttingen, Germany, located at 51^∘33^′3^′′ N, 9^∘57^′2^′′ E. The flux tower was installed at an altitude of 230 m above sea level on a flat area of about 50 m by 50 m situated on a south–southeast-facing hill, with a slope angle of 5^∘ and length of 800 m. Vegetation height of the grass was 0.2 m during the experiment. Vegetation further comprised patches of bushes and trees with a minimum distance from the flux tower of 50 m (west of the tower).

2.3 Experimental period

The TEA flux measurements presented in this study were conducted from 4 to 10 April 2015. After a cold and wet month of March, this period was characterized by increasing physiological activity of the grasses due to increasing light availability, increasing temperatures during the day, less frequent frost events, and increasing CO₂ fluxes.

Prior to the flux experiment, the TEA instrument and method were further developed and tested in the field, with continuous operation starting on 5 March 2015. The TEA deployment continued after 10 April until 17 June 2015. However, the frequent charging and discharging of the air sampling bags led to material fatigue and progressive leakage. Therefore, no meaningful flux measurements are available after 10 April 2015. The period from 10 April to 17 June 2015 was used for testing different kinds of air bags and for further developing the TEA method. Altogether, the TEA air sampling system in its current form was in continuous operation from 5 March to 17 June 2015, corresponding to more than 5000 30 min TEA flux sampling intervals.

2.4 Instrumentation

The instrumentation used for meteorological measurements, TEA flux measurements, and EC flux measurements and the respective variables measured are listed in Table 1.

Table 1 Instrumentation used for turbulent flux and meteorological measurements. Manufacturer key: Gill Instruments (Lymington, UK), Li-COR Environmental Inc. (Lincoln, Nebraska, USA), Kipp & Zonen B.V. (Delft, The Netherlands), Bosch Sensortec (Stuttgart, Germany), Adolf Thies GmbH & Co. KG (Göttingen, Germany), Imko Mikromodultechnik GmbH (Ettlingen, Germany), Hukseflux Thermal Sensors B.V. (Delft, the Netherlands).

Meteorological variables (Table 1) were logged using a DL16 data logger (Adolf Thies GmbH & Co. KG Göttingen, Germany). All raw data needed for TEA and EC flux measurements, including the sonic anemometer data and data from the two infrared gas analyzers LI-6262 and LI-7500, were synchronized and logged on the central TEA controller. Using a mobile network link, raw data were continuously mirrored to a central server for archival and flux processing.

2.5 Meteorological measurements

The following set of meteorological variables was measured on site to support the computation and interpretation of turbulent fluxes: global radiation, photosynthetically active radiation (PAR), net radiation, air temperature at 2 m a.g.l., air pressure, relative humidity at 2 m a.g.l., precipitation, wind velocity at 2 m a.g.l., wind direction at 2 m a.g.l., soil temperature at 0.3 m b.g.l. (three probes), soil moisture at 0.3 m b.g.l.(three probes), and soil heat flux.

2.6 Experimental data

The above-described meteorological measurements, as well as the turbulent flux measurements obtained by the true eddy accumulation and the two eddy covariance instrumental setups, have been published via open access (Siebicke and Emad, 2019) to foster open science in general and the scientific advancement of true eddy accumulation in particular.

2.7 Coordinate systems and net ecosystem exchange

If the trace gas source and sink strength of the ecosystem is of interest, as is typically the case when investigating biological, physiological, or biogeochemical processes or deriving trace gas budgets, then the total flux in and out of the ecosystem needs to be determined. For the exchange of the ecosystem with the atmosphere, the concept of a virtual control volume is often used. Net ecosystem exchange (NEE), i.e., the net flux across the surfaces of this control volume, can be written as (e.g., Aubinet et al., 2003; Siebicke et al., 2012)

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{NEE}& ={\displaystyle \frac{\mathrm{1}}{V_{\mathrm{m}}}}{\underset{\mathrm{0}}{\int }}^h\underset{\text{I}}{\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}\right)}\mathrm{d}z+{\displaystyle \frac{\mathrm{1}}{V_{\mathrm{m}}}}\underset{\mathrm{II}}{{\left(\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}\right)}_h}\\ & +{\displaystyle \frac{\mathrm{1}}{V_{\mathrm{m}}}}{\underset{\mathrm{0}}{\int }}^h\left(\stackrel{\mathrm{‾}}{w}\underset{\mathrm{IIIa}}{\left(z\right){\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial z}}}+\stackrel{\mathrm{‾}}{c}\underset{\mathrm{IIIb}}{\left(z\right){\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{w}}{\partial z}}}\right)\mathrm{d}z\\ & +{\displaystyle \frac{\mathrm{1}}{V_{\mathrm{m}}}}{\underset{\mathrm{0}}{\int }}^h\underset{\mathrm{IV}}{\left(\stackrel{\mathrm{‾}}{u}\left(z\right){\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial x}}+\stackrel{\mathrm{‾}}{v}\left(z\right){\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial y}}\right)}\mathrm{d}z\\ & +{\displaystyle \frac{\mathrm{1}}{V_{\mathrm{m}}}}{\underset{\mathrm{0}}{\int }}^h\underset{\text{V}}{\left({\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{u^{\prime }{c}^{\prime }}\right)}{\partial x}}+{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{v^{\prime }{c}^{\prime }}\right)}{\partial y}}\right)}\mathrm{d}z,\end{array}\end{array}$$

with the molar volume of dry air V_m, CO₂ concentration c, time t, horizontal distances x and y, vertical distance above ground z, height of the control volume h, horizontal wind velocity u along the x direction, horizontal wind velocity v along the y direction, and vertical wind velocity w along the z direction. Overbars denote temporal means, and primes denote the temporal fluctuations relative to the temporal mean.

The terms on the right hand side of Eq. (1) are the change of storage (I), the vertical turbulent flux (II), vertical advection (IIIa), vertical mass flow from the surface, e.g., due to evaporation (IIIb) according to Webb et al. (1980), horizontal advection (IV), and flux divergence (V). The form of NEE presented in Eq. (1) excludes the horizontal variation of the vertical turbulent flux and the horizontal variation of vertical advection. Most flux measurements typically only determine the vertical turbulent flux density (term II) and sometimes the storage flux density (term I), neglecting the remaining terms due to a lack of spatially distributed information.

The choice of the reference coordinate system (Finnigan, 2004) is important for the attribution of the total flux to its components (Eq. 1) and therefore for the interpretation of turbulent flux density measurements relative in their ability to approximate the net ecosystem exchange. If NEE is to be assessed, and available flux observations are restricted to the turbulent vertical flux density at a single location above the ecosystem, a reference coordinate system is needed which minimizes the remaining flux terms. Sun et al. (2007, Table 1) summarize coordinate systems. In streamline coordinates (Finnigan, 2004; Sun et al., 2007), the long-term flow is tangential to long-term streamlines. This means that, in streamline coordinates, the velocity normal to the streamlines becomes zero, implying that the long-term vertical advection vanishes. There are various methods for coordinate rotation (Finnigan et al., 2003), i.e., the transformation of the wind measurements from the coordinate reference frame of the sonic anemometer to the coordinate system of the flux measurements, also known as tilt correction (Tanner and Thurtell, 1969; Hyson et al., 1977; Kaimal and Finnigan, 1994; McMillen, 1988; Paw U et al., 2000; Wilczak et al., 2001).

Over flat terrain or planar terrain with a uniform slope, the mean streamlines close to the surface approximately follow the terrain surface. The planar fit method (Paw U et al., 2000; Wilczak et al., 2001) is often used to obtain long-term streamline coordinates. In contrast to the double rotation method (Kaimal and Finnigan, 1994), which nullifies the mean vertical wind velocity $\stackrel{\mathrm{‾}}{w}$ of the flux integration interval, planar fit rotated $\stackrel{\mathrm{‾}}{w}$, using the original formulation of Wilczak et al. (2001), is typically small but not zero. Even the long-term $\stackrel{\mathrm{‾}}{w}$ only becomes zero if mean streamlines are planar, and there is no instrumental bias in measurements of w. Nonzero $\stackrel{\mathrm{‾}}{w}$ would imply existence of a vertical advection term proportional to $\stackrel{\mathrm{‾}}{w}$ in the presence of vertical trace gas concentration gradients, which for CO₂ typically exist close to the surface.

A variant to the planar fit method proposed by Van Dijk et al. (2004) removes velocity bias relative to the flux integration interval, addressing instrument offsets. This procedure can lead to local misalignment of streamlines for nonplanar mean flow fields. Nullifying $\stackrel{\mathrm{‾}}{w}$ over the flux integration period would formally remove vertical advection terms from the flux budget equation (Eq. 1). However, this procedure would also still ignore the effect of misalignment of the reference coordinate system and the streamlines over the flux integration interval on the vertical turbulent flux. The mismatched length and timing of the planar fit period relative to the shorter individual flux integration intervals acts as a high pass filter and results in loss of low-frequency flux contributions and in unresolved distortion of co-spectra of the shorter flux intervals (Finnigan et al., 2003).

On a related matter, Siebicke et al. (2012) performed an explicit treatment of the length and timing of the reference period for planar fit rotation. They demonstrated changes of up to 50 % of advective CO₂ fluxes over forest depending on the window size of a new serial planar fit approach.

Over complex nonplanar terrain, the mean streamlines are not tangential to a plane. Even over planar surfaces, streamlines further away from the surface may not be tangential to the terrain surface due to vertical velocity divergence (Sun et al., 2007). For curved streamlines, other terms of the mass flow equation (Finnigan, 2004), in addition to vertical turbulent flux and vertical advection, become important. For curved streamlines, horizontal trace gas advection is not proportional to the gradient of trace gas concentrations along the streamlines.

Several authors have suggested variants to the planar fit method to account for steep slopes, where buoyancy forces are no longer normal to the terrain surface (Oldroyd et al., 2016), for obstructed flow (Griessbaum and Schmidt, 2009), and for complex topography, where $\stackrel{\mathrm{‾}}{w}$ becomes a function of wind azimuth angle. Consequently, several studies apply the planar fit rotation separately for different wind direction sectors (Foken and Napo, 2008; Yuan et al., 2011, and others). However, this introduces discontinuities in $\stackrel{\mathrm{‾}}{w}$ at wind direction sector limits and in the definition of the reference coordinate system.

Here we propose that a more general approach avoiding directional discontinuities would be to fit a surface rather than a plane, where the curvature of the surface adapts to long-term streamlines as a function of one or more parameters, i.e., wind direction and optionally other variables, such as horizontal wind velocity. A related approach has recently been proposed by Ross and Grant (2015), who also suggest tilt correction as a continuous function of wind direction instead of the relatively common discrete sectoral approach used currently. Siebicke et al. (2012) showed that the effect of atmospheric stratification, friction velocity, stationarity, and integral turbulence characteristics (Foken et al., 2004) on sectoral planar fit rotation was small relative to the wind direction effect over a forest.

In the presence of flow distortion due to obstacles, terrain features, towers and sensor mounts, or non-omnidirectional sonic anemometer designs (Li et al., 2013), distorted sectors need to be excluded from the definition of the coordinate system and subsequent flux derivations, unless distortions are characterized and corrected for (Van Dijk et al., 2004; Griessbaum and Schmidt, 2009; see also Sect. 2.8, current study).

All the above considerations on coordinate rotation apply to both the EC and EA methods, respectively, including their derivatives. However, there is one conceptual difference: in EC, high-frequency observations of both the scalar and the wind are measured and stored for post-processing. Therefore it is feasible to fit the coordinate rotation plane to all observations from one or several 30 min flux integration intervals, including observations which were obtained before and after any wind measurement which is to be rotated. This includes the possibility to center the planar fit period around the time of interest.

Conversely, in EA, in absence of high-frequency scalar observations, any decision on the reference coordinate system becomes final on obtaining and mixing individual high-frequency air samples, precluding post-processing and any reconsideration of the coordinate system. In EA, the reference period defining the coordinate system necessarily cannot coincide and never fully overlap with the flux integration interval for any sample but the last one in the flux integration interval, if any. Due to this conceptual difference, the contribution of other flux terms, in particular vertical advection, may not be identical for EC and EA if not using the same reference period to align the coordinate system to the mean streamlines.

Nonvanishing mean vertical wind velocities $\stackrel{\mathrm{‾}}{w}$ over the flux integration interval can nominally be removed in EC through subtraction of $\stackrel{\mathrm{‾}}{w}$ from instant vertical velocities (Van Dijk et al., 2004) and in EA through application of the volume mismatch correction (Turnipseed et al., 2009); see Eq. (11) above. However, distortion of co-spectra (Finnigan et al., 2003) remains uncorrected.

We distinguish in evaluating the implications of discussed deviations of the real flow from any chosen ideal reference conditions: (i) the case of deploying turbulent vertical flux measurements to estimate net ecosystem exchange and (ii) the case of comparing the EA method and instruments side by side to the chosen reference method EC for assessing whether the EA method's physical air sampling principle produces comparable results to the mathematical computation of covariances for the EC method. This study is concerned with the latter case only. Most of the above issues of coordinate frames and the spatiotemporal variability of the flow field afflict both methods alike. Only the unavoidable differences in the application of the coordinate rotations between the two methods, i.e., the nonmatching rotation periods, need to be of concern when evaluating the relative performance of the two turbulent flux observation methods. To eliminate this remaining difference, identical rotation procedures and planar fit reference periods need to be applied to both EA and EC, accepting the EA version as the reference.

2.8 Flow distortion and angle of attack correction

The physical structure of sonic anemometer probes distorts the air flow they intend to measure (Wyngaard, 1981), introducing systematic errors in flux measurements. Measurement errors, due to probe-induced flow distortion and self-sheltering of ultrasonic transducers (Gash and Dolman, 2003), depend on the “angle of attack” (Kaimal and Finnigan, 1994), i.e., the angle between horizontal and the instantaneous wind vector. A wind tunnel calibration for anemometer models R2 and R3 (Gill Instruments Ltd., UK) was provided by van der Molen et al. (2004) and updated by Nakai et al. (2006). The representativity of the wind tunnel calibrations for turbulent conditions in the field has been questioned (Högström and Smedman, 2004) and is still under debate (Huq et al., 2017). Nakai and Shimoyama (2012) proposed an improved correction based on field measurements under turbulent conditions for the R3 and Windmaster models (Gill Instruments Ltd., UK). There is still no consensus on whether this correction should be applied, and care must be taken as the correction applies to certain instrument models (Gill Windmaster) and serial numbers only.

In the current study, which uses two R3-type anemometers (Gill Instruments Ltd., UK), we do not apply any angle-of-attack correction because (i) the applicability of the wind tunnel calibration (Nakai et al., 2006) may or may not be applicable; (ii) there is contrasting information on the applicability of the calibration under turbulent conditions (Nakai and Shimoyama, 2012), specifically to the R3 model (recommended for R3 by original authors but not according to later information from Gill Instruments, UK, and LI-COR Env., USA); (iii) no angle-of-attack correction was available in the current TEA system software at the time of the field experiment nor can the TEA flux measurements be post-processed to include the correction. For the above reasons, no angle-of-attack correction was applied to the presented results of the current study, neither to TEA nor to EC fluxes. However, we did assess the impact of the angle-of-attack correction on EC CO₂ fluxes for the two sonic anemometers.

2.9 True eddy accumulation (TEA) flux measurements

2.9.1 TEA method

The true eddy accumulation method determines the flux, $\stackrel{\mathrm{‾}}{wc}$, of a scalar (such as a trace gas) as the sum of the covariance, $\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}$ of the scalar, c, and the vertical wind velocity, w, and the product of the time averages of scalar and vertical wind velocity, $\stackrel{\mathrm{‾}}{w}\stackrel{\mathrm{‾}}{c}$, as

$$\begin{array}{}\text{(2)}& \stackrel{\mathrm{‾}}{wc}=\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}+\stackrel{\mathrm{‾}}{w}\stackrel{\mathrm{‾}}{c},\end{array}$$

where overbars denote time averages over the averaging period, T_avg, and primes denote fluctuations from the mean. Eq. (2) is analog to the eddy covariance (EC) method.

However, in contrast to EC, which requires high-frequency observations of the scalar and the vertical wind velocity and mathematical derivation of the covariance through post-processing, in the case of TEA, the separate sampling of the wind and scalar time series is replaced by physically collecting separate air samples of updrafts and downdrafts proportionally to the magnitude of the vertical wind velocity. The TEA flux over a given averaging period, T_avg, can thus be obtained as (Desjardins, 1977; Hicks and McMillen, 1984)

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{‾}}{wc}={\displaystyle \frac{\mathrm{1}}{T_{\mathrm{avg}}}}\underset{\mathrm{0}}{\overset{T_{\mathrm{avg}}}{\int }}({\mathit{\delta }}^{+}cw+{\mathit{\delta }}^{-}cw)\mathrm{d}t,\end{array}$$

where ${\mathit{\delta }}^{+}=\mathrm{1}$ when $w>\stackrel{\mathrm{‾}}{w}$ and 0 when $w<\stackrel{\mathrm{‾}}{w}$, and ${\mathit{\delta }}^{-}=\mathrm{1}$ when $w<\stackrel{\mathrm{‾}}{w}$ and 0 when $w>\stackrel{\mathrm{‾}}{w}$. The amount of air, cw, sampled per unit time, dt, contains the molar fraction of the scalar of interest, c.

Assuming ideal conditions such that the mean vertical wind velocity over the integration period, $\stackrel{\mathrm{‾}}{w}$, was zero, the mean term, $\stackrel{\mathrm{‾}}{w}\stackrel{\mathrm{‾}}{c}$, becomes zero and the scalar flux, F_c, becomes

$$\begin{array}{}\text{(4)}& F_c=\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}\end{array}$$

in kinematic units of meters per second (m s⁻¹). Multiplying by moist air density ρ, we obtain the constituent mass flux, F_c, per unit area and unit time in units of kilograms per square meter per second ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) as

$$\begin{array}{}\text{(5)}& F_c=\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}\stackrel{\mathrm{‾}}{\mathit{\rho }}.\end{array}$$

Given $\stackrel{\mathrm{‾}}{w}=\mathrm{0}$, the volume of the air samples accumulated in the updraft reservoir is identical to the volume of the downdraft reservoir given

$$\begin{array}{}\text{(6)}& w^{+}+w^{-}=\stackrel{\mathrm{‾}}{w}=\mathrm{0},\end{array}$$

where w⁺ is vertical wind larger than $\stackrel{\mathrm{‾}}{w}$, and w⁻ is vertical wind smaller than $\stackrel{\mathrm{‾}}{w}$.

A practical implementation of TEA then determines the scalar flux, F_c, as half of the difference between the mole fraction of the scalar in the updraft reservoir, $\stackrel{\mathrm{‾}}{c^{+}}$, and the mole fraction of the scalar in the downdraft reservoir, $\stackrel{\mathrm{‾}}{c^{-}}$, multiplied by the mean of the absolute value of vertical wind velocity, assuming $\stackrel{\mathrm{‾}}{w}=\mathrm{0}$.

$$\begin{array}{}\text{(7)}& F_c={\displaystyle \frac{\stackrel{\mathrm{‾}}{\left|w\right|}}{\mathrm{2}}}(\stackrel{\mathrm{‾}}{c^{+}}-\stackrel{\mathrm{‾}}{c^{-}})\stackrel{\mathrm{‾}}{\mathit{\rho }}\end{array}$$

in units of kilograms per square meter per second ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) or

$$\begin{array}{}\text{(8)}& F_c={\displaystyle \frac{\stackrel{\mathrm{‾}}{\left|w\right|}}{\mathrm{2}}}(\stackrel{\mathrm{‾}}{c^{+}}-\stackrel{\mathrm{‾}}{c^{-}}){\displaystyle \frac{\mathrm{1}}{V_m}}\end{array}$$

in units of moles per square meter per second ($\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), with the molar volume of air, V_m (m³ mol⁻¹), according to the ideal gas law,

$$\begin{array}{}\text{(9)}& V_m={\displaystyle \frac{RT}{P}},\end{array}$$

with temperature, T, pressure, P, and the ideal gas constant, R. Note that Eqs. (7) and (8) write $\stackrel{\mathrm{‾}}{\left|w\right|}$, not $\left|\stackrel{\mathrm{‾}}{w}\right|$; i.e., the order of math operations is important.

For comparison of TEA and REA methods, we include here the formulation of the flux using the REA method (Businger and Oncley, 1990), expressed as mole fraction measurements of the constituent, c,

$$\begin{array}{}\text{(10)}& F_c=\mathit{\beta }{\mathit{\sigma }}_w\left(\stackrel{\mathrm{‾}}{c^{+}}-\stackrel{\mathrm{‾}}{c^{-}}\right)\stackrel{\mathrm{‾}}{\mathit{\rho }},\end{array}$$

with the standard deviation of vertical wind velocity, σ_w, and the REA flux-variance similarity factor, β.

In any practical flux measurement application the observed mean vertical wind velocity over the integration period is likely unequal to zero, i.e., $\stackrel{\mathrm{‾}}{w}\ne \mathrm{0}$. In eddy covariance the assumption of $\stackrel{\mathrm{‾}}{w}=\mathrm{0}$ is satisfied in post-processing once all observations of the entire integration period are available. This is commonly achieved by rotating the coordinate frame of the wind measurements to minimize or even nullify $\stackrel{\mathrm{‾}}{w}$ followed by subtraction of $\stackrel{\mathrm{‾}}{w}$ from individual vertical wind velocity measurements, w, so $\stackrel{\mathrm{‾}}{w}$ becomes zero.

Conversely, for TEA, knowledge of the mean vertical wind velocity over the flux averaging period, $\stackrel{\mathrm{‾}}{w}$, is required at any time throughout the averaging period, in order to be able to classify vertical winds as updrafts or downdrafts and accumulate air samples in the corresponding reservoirs. At any time during the averaging interval only the past and present vertical wind measurements are known. Therefore, any attempt to obtain $\stackrel{\mathrm{‾}}{w}=\mathrm{0}$ needs to rely in part on an estimate of the mean vertical wind velocity of the entire averaging period without knowledge of future observations from the present through to the end of the averaging period. In practice, this situation can lead to $\stackrel{\mathrm{‾}}{w}\ne \mathrm{0}$, resulting in unequal sample volumes accumulated over the averaging period in the updraft and downdraft reservoirs.

Following Turnipseed et al. (2009) the flux needs to be corrected by a term accounting for the mismatch between the volume of accumulated updrafts, V⁺, and downdrafts, V⁻, respectively, for $V^{+}\ne V^{-}$. The flux due to a mismatch of volumes V⁺ and V⁻ is

$$\begin{array}{}\text{(11)}& F_{c,\mathrm{volume}\mathrm{\_}\mathrm{mismatch}}=\stackrel{\mathrm{‾}}{w}\left({\displaystyle \frac{(\stackrel{\mathrm{‾}}{c^{+}}+\stackrel{\mathrm{‾}}{c^{-}})}{\mathrm{2}}}-\stackrel{\mathrm{‾}}{{}^{v}c}\right){\displaystyle \frac{\mathrm{1}}{V_m}},\end{array}$$

where $\stackrel{\mathrm{‾}}{w}$ is the mean vertical velocity for the averaging period. The volume mismatch correction term is the difference between the unweighted mean density of the reservoirs, $(\stackrel{\mathrm{‾}}{c^{+}}+\stackrel{\mathrm{‾}}{c^{-}})/\mathrm{2}$, and the volume weighted mean density,

$$\begin{array}{}\text{(12)}& \stackrel{\mathrm{‾}}{{}^{v}c}={\displaystyle \frac{(\stackrel{\mathrm{‾}}{c^{+}}{V}^{+}+\stackrel{\mathrm{‾}}{c^{-}}{V}^{-})/\mathrm{2}}{(V^{+}+V^{-})/\mathrm{2}}},\end{array}$$

weighted by the updraft and downdraft volumes, V⁺ and V⁻, respectively.

Inserting Eq. (11) in Eq. (8) yields the volume-mismatch-corrected TEA flux,

$$\begin{array}{}\text{(13)}& F_c=\left({\displaystyle \frac{\stackrel{\mathrm{‾}}{\left|w\right|}}{\mathrm{2}}}\right(\stackrel{\mathrm{‾}}{c^{+}}-\stackrel{\mathrm{‾}}{c^{-}})+\stackrel{\mathrm{‾}}{w}\left({\displaystyle \frac{(\stackrel{\mathrm{‾}}{c^{+}}+\stackrel{\mathrm{‾}}{c^{-}})}{\mathrm{2}}}-\stackrel{\mathrm{‾}}{{}^{v}c}\right)){\displaystyle \frac{\mathrm{1}}{V_m}}.\end{array}$$

In practical applications, the instant sampling volume per unit time, V_i, is related to instant vertical wind velocity, w_i, through a proportionality factor, k, as

$$\begin{array}{}\text{(14)}& V_i=k\left|w_i\right|.\end{array}$$

2.9.2 Correction of trace gas mole fractions for the effects of water vapor

When measuring trace gases such as CO₂ in moist air with infrared gas analyzers, two corrections are required to remove the effect of water vapor on the measurement of the mole fraction of the trace gas. The first correction accounts for pressure broadening due to the presence of water vapor. This is known as the “pressure broadening correction” (Licor, 1996, Li-COR LI-6262 manual, Sect. 3.5, 25 pp., Eq. 3–30).

The second correction accounts for the dilution of the trace gas by water vapor in the sample. This correction is known as “dilution correction” and is required to convert the wet mole fraction, c_wet, to the dry mole fraction, c_dry. The dry mole fraction is needed for calculating the trace gas flux, F_c, based on Eq. (5). The dry mole fraction is obtained from the wet mole fraction following the instructions for the infrared gas analyzer (Licor, 1996):

$$\begin{array}{}\text{(15)}& c_{\mathrm{s}}^{\mathrm{wr}}=c_{\mathrm{s}}^{\mathrm{ws}}\left({\displaystyle \frac{\mathrm{1}-X_{\mathrm{w},\mathrm{ref}}/\mathrm{1000}}{\mathrm{1}-X_{\mathrm{w}}/\mathrm{1000}}}\right),\end{array}$$

where X_w is the mole fraction of water vapor in the sample cell, X_w,ref is the water vapor mole fraction in the reference cell, $c_{\mathrm{s}}^{\mathrm{ws}}$ is the actual mole fraction of the trace gas in the sample cell diluted by X_w, and $c_{\mathrm{s}}^{\mathrm{wr}}$ is the equivalent sample cell CO₂ mole fraction if it were diluted by X_w,ref.

2.9.3 TEA instrumentation and technical implementation

The TEA instrumentation used in this study was developed with particular attention to accurate and dynamic sampling of air and to real-time processing of wind data. The system was further designed to minimize time lags and jitter in wind data processing and in air sampling and to minimize dead volumes in the gas sampling system.

Vertical and horizontal wind velocities and the sonic temperature were measured using an ultrasonic anemometer R3 (Gill Instruments Ltd.), the same type which was also used for the side-by-side eddy covariance reference flux measurements. Wind velocity data were logged at a 10 Hz frequency.

Instant observations of vertical wind velocity were subject to real-time coordinate rotation to align the coordinate system of the sonic anemometer with the mean streamlines prior to controlling the sampling of air into updraft or downdraft reservoirs using the planar fit method (Wilczak et al., 2001). To the best of our knowledge, this study is the first eddy accumulation study to deploy a real-time coordinate rotation based on the planar fit rotation (Wilczak et al., 2001). A moving window of 1 d was used to estimate planar fit rotation coefficients, with an update frequency of the rotation coefficients of 30 min. The coefficients were then applied to rotate the instant raw wind measurements, w_i, 10 times per second.

To minimize $\stackrel{\mathrm{‾}}{w}$ over the flux averaging period the following procedure was applied: the mean vertical wind velocity of the current accumulation interval was approximated by the mean rotated vertical wind velocity over the most recent samples over a period with length equal to the length of the accumulation period, in this case 30 min. This estimate of mean vertical wind velocity was updated every 2 min and subtracted from every instantaneous vertical wind velocity measurement after coordinate rotation, i.e., 10 times per second.

The decision on updrafts and downdrafts was based on the sign of rotated w_i. Sample volumes were computed following Eq. (14). k was determined such that instant flow rates would not exceed the maximum possible flow rate of 3 L min⁻¹, with a probability of 99 % based on absolute wind data in the period from 30 min ago to present. The proportionality factor, k, which was based on the 99 % quantile multiplied by a factor of 2, was updated every 30 min.

Figure 4 Schematic functional flowchart of true eddy accumulation system. Sampling of air with atmospheric constituents (scalars) shown on the left, sampling of wind vector shown on the right. Updrafts and downdrafts were sampled and stored via separate lines. One single analyzer was used with the sample supply alternating between the updraft and downdraft reservoirs every 150 s (see Fig. 6). During a particular 30 min period, bag set A was being filled with updraft and downdraft samples in the updraft and downdraft bags, respectively. At the same time bag set B was being analyzed and discharged, with the analysis alternating between updrafts and downdrafts. Every 30 min, the operation of bag sets A and B, either filling or discharging, was swapped. The vertical wind velocity data (right) control the air sampling via instantaneous wind measurements and wind statistics.

Figure 5 Hydraulic and functional schematic of the true eddy accumulation system with components layout, properties, and operating conditions. System components shown are the piping layout, mass flow controllers (MFCs), pumps, valves, filters, reservoirs, and the gas analyzer. Piping lengths and volume of air reservoirs are shown alongside the components. The operating pressure is shown for three distinct regions of the system below the system layout (colored bars) together with the average transit times of air through specific sections of the system. The air sampling bags are marked with “up” for collecting updrafts and “down” for collecting downdrafts. Bag sets A and B alternate their function (charge, discharge) every 30 min.

Air inlets were co-located with the wind measurement and positioned 18 cm below and 2 cm beside the center of the sonic anemometer. Two separate air sampling lines were used, one for obtaining the updraft samples and one for the downdraft samples, respectively. In contrast to many previous eddy accumulation studies, which used a single air inlet and a 3∕2-way valve to direct the samples towards the appropriate updraft or downdraft reservoir, the current design with two separate sampling lines avoids any undesired mixing of updrafts with downdrafts in the system. A detailed technical description and layout of the true eddy accumulation system with the two separate sample lines is presented in Fig. 5, including piping layout and system components, operating pressure conditions, pressure drops, dead volumes, and transit times of air samples through the system.

The intake of air was controlled by fast-response mass flow controllers with a dynamic response resolving turbulent eddies at 10 Hz. The mass flow controllers used for this study were calibrated using conventional thermal mass flow controllers (Vögtlin red-y smart series). The accuracy of the new dynamic mass flow controllers was equal to or better than 0.3 %, which corresponds to the accuracy of the conventional mass flow controller model used for calibration.

The new type of mass flow controller was tested for potential leaks, both relative to the ambient air and also for potential flow during times when the respective sampling line was inactive. The tests showed that the mass flow controller (MFC) was virtually leak-free. The combined leak rate of the MFCs, the pumps, the filters, tubing, and fittings was 0.0086±0.0003, expressed as the leak rate relative to the average inlet flow rate. In terms of flux uncertainty this would mean that a theoretical maximum of 0.86±0.03 % and, considering the other system components, likely far less of the flux uncertainty would be related to potential leakage or undesired nonzero flow of the mass flow controllers.

Air sampling lines were made of Teflon with a 6 mm outer and 4 mm inner diameter and a length of 5 m between intake and accumulation reservoirs. The air was filtered before entering the pumps and the bag reservoirs using PTFE membrane Gelman Acro Disc filters with a 50 mm diameter and a 2 µm pore size. Another filter was placed directly upstream of the gas analyzer.

At any time, one of the air sampling lines was active, with the selection of the line depending on the sign of the vertical wind velocity. The wind was measured at a frequency of 10 Hz using the sonic anemometer. With every new reading of vertical wind velocity, i.e., every 100 ms, the selection of the active inlet (updraft or downdraft) was updated depending on the sign of the vertical wind velocity and an air sample obtained with mass proportional to vertical wind velocity.

In contrast to common practice in many relaxed eddy accumulation studies, which typically define a minimum vertical wind velocity for air sampling, in the current design, air samples were obtained for all magnitudes of vertical wind velocity (except for the 0.5 % most positive and most negative values of w, respectively, as per the definition of k above).

Air samples were collected using two variable speed separate brushless DC membrane pumps (KNF Neuberger GmbH, Germany) with a maximum flow rate of 3 L min⁻¹, each, feeding air into the bag reservoirs at flow rates between 0 and 3 L min⁻¹.

Air was collected in lab grade, chemically inert Alumini^® air sample bags (Westphalen AG, Germany) with a volume of 28 L. The composite wall of the bags was made of (from outside to inside) polyethylene terephthalate (PET), polyethylene (PE), aluminum (ALU), oriented polyamide (OPA), and polyethylene (PE).

The layout of the TEA system is shown in Fig. 4. The system was designed for continuous operation, with continuous sampling of air and continuous on-site gas analysis. Air was collected in bags over periods of 30 min. Subsequently the air was analyzed over the following 30 min periods. A total of four bags was used at any time. A set of two bags (marked with “A” in Fig. 4) was charged with samples over 30 min, with one bag accumulating samples corresponding to updrafts (marked “updraft”) and one bag accumulating samples corresponding to downdrafts (marked “downdraft”). A second set of two bags (marked with “B” in Fig. 4), which contained air samples from the previous 30 min interval, was discharged and analyzed in parallel to the filling of the first set of bags. Every 30 min the two bag sets A and B would swap their function from being charged to being discharged. Note that data analysis revealed a leak in one of the bag sets, so that the data of every other half hour had to be discarded.

The accumulated air was analyzed for molar density of CO₂ (µmol mol−1) and of water vapor (mmol mol−1) using a dual cell infrared gas analyzer, model LI-6262 (Li-COR Env. Inc, USA). Air samples were discharged by a membrane pump (KNF Neuberger, Germany) from the sample bags through the sample cell of the gas analyzer at a flow rate of 0.6 L min⁻¹. The sampling frequency of the gas analyzer was 1 Hz. The reference cell was purged by dry and CO₂-free zero gas obtained by circulating air through a scrubber filled with soda lime and Drierite desiccant. 3∕2-way solenoid valves were used to select the appropriate gas bag for gas analysis.

During any 30 min period, the gas analysis alternated between sampling the updraft and downdraft reservoir, respectively (Fig. 6). Each bag was sampled for 150 s at a time. The first eight measurements periods of 150 s each were used for further analysis, resulting in four replicate measurements of the gas densities of the updraft reservoir and likewise four replicates of the downdraft reservoirs. Gas density measurements were tagged by the TEA controller with the respective active channel, either updraft or downdraft reservoir, for flux processing. The alternating sampling sequence lasted 1200 s. Over the remaining 600 s of the 30 min period, the remaining air was discharged from the bags to the atmosphere at a flow rate of 1.5 L min⁻¹ until depletion.

Figure 6 CO₂ dry molar fraction measured as a sequence alternating between updraft and downdraft reservoirs. First the updraft reservoir was measured for 150 s, then the downdraft reservoir for 150 s, discarding the initial 30 s of measurements from each block. This sequence was repeated four times. During the remaining 10 min of the 30 min period any remaining air in the reservoirs was purged and not used for analysis (gray shaded area). The means of the despiked updraft and downdraft dry molar fractions are indicated by a red and blue solid line, respectively. The dashed lines indicate the mean ±1 standard deviation. The start date of the time series is 10 April 2015 13:00:00 UTC.

The fully automatic TEA system was controlled by an ARM-based single-board Linux computer of the type “Raspberry Pi” (Raspberry Pi Foundation, UK), the “TEA controller”. All sensor measurement data, including the wind and gas density measurements, were synchronized, logged, and processed on the same TEA controller.

The following raw data were logged for subsequent turbulent flux calculations: horizontal and vertical wind velocity components, $u,v,w$; sonic temperature, T_s; CO₂ and H₂O molar densities; analyzer cell temperature and cell pressure; ambient air temperature, T; and air pressure, P. Further data on the state of the TEA sampling system and the analysis were logged for attribution of the gas analyzer measurements to updraft and downdraft samples, for the selection of the bag sets, and for system monitoring and quality control.

The energy-efficient TEA system of the current study consumed 15 W of electrical power (excluding the gas analyzer), about 10 W of which was used by the three pumps. The pumping power required for the current TEA system was 2 orders of magnitude smaller than that of fast-flow closed-path eddy covariance systems using infrared gas analyzers or laser spectrometers of ca. 1 kW. The additional power consumption of typical current laser spectrometers would be on the order of 0.25 kW. The difference in pumping power between TEA and EC scales with the flow rate and sample cell vacuum of the EC application and is therefore even more important for most closed-path laser spectrometers than for infrared gas analyzers.

2.9.4 TEA flux computations

Turbulent fluxes of CO₂ were calculated from raw data of horizontal, u,v, and vertical, w, wind velocity components, sonic temperature, T_s (10 Hz), CO₂ dry mole fraction (converted from wet mole fraction) and H₂O dry mole fraction (converted from wet mole fraction) (1 Hz), ambient air temperature, T (10 min resolution), and pressure, P (10 min resolution).

Raw gas density measurements were processed prior to flux computations in order to filter the data for noise and aggregate individual readings to single representations of the gas density, one for the updraft and one for the downdraft reservoir, during any one 30 min period. Blocks of 150 s (see Fig. 6) of measurements at 1 Hz were filtered and aggregated to 30 min values (see Sect. 3, Fig. 11 for results). The following statistically robust procedure was used:

• Raw voltage signals of CO₂ and H₂O were converted to physical units of micromoles per mole (µmol mol⁻¹), and the band broadening effect of pressure on CO₂ observations was corrected for according to Licor (1996).

• The CO₂ wet mole fraction was converted to the dry mole fraction.

• Gas density data were checked for plausibility based on preset minimum and maximum values.

• The initial 30 s (dead-band filter) after switching channels was omitted to allow for purging of the shared gas handling components, i.e., valves, sample line, and gas analyzer sample cell.

• Raw data were de-spiked (spike filter) using the function “despike” from the R package “oce” (Kelley and Richards, 2017), replacing discarded values with the median of remaining values. The method identifies spikes with respect to a “reference” time series and replaces these spikes with the reference value, i.e., here the median.

• The time series were smoothed (smoothing filter) using the function “loess” (Cleveland et al., 1992) from the R package “stats” (R Core Team, 2017). The loess function fits a polynomial surface determined by one or more numerical predictors, using local fitting.

• Stable readings were selected (stationarity filter) by limiting maximum permissible gradients between individual samples for channel-specific data blocks of 150 s (120 s remaining after dead-band filter), with a maximum permissible change of 0.002 $\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$.

• Sufficient availability of data after filtering was checked for (availability filter): a data block was discarded if fewer than 30 (out of 120) values remained available after the stationarity filter.

• The remaining filtered data were aggregated per channel-specific data block of 150 s using the median function.

• The four replicates were aggregated per channel into a single value per 30 min as the weighted mean of the four samples, weighted by the number of accepted raw measurements in each of the four replicates, separately for updrafts and downdrafts.

• To quantify precision of the CO₂ molar fraction measurements using the LI-6262 infrared gas analyzer, the minimum and maximum molar fractions over the four replicated samples were estimated, separately for updrafts and downdrafts and propagated into minimum and maximum flux estimates, respectively.

Fluxes were then calculated through the following:

• plausibility check of wind data based on preset minimum and maximum values;

• de-spiking of raw wind data (Vickers and Mahrt, 1997);

• computation of the mole fraction difference between updraft and downdraft reservoirs per 30 min period;

• computation of uncorrected turbulent fluxes for 30 min intervals according to Eq. (8);

• computation of turbulent fluxes for 30 min intervals, corrected for volume mismatch between updraft and downdraft reservoirs, according to Eq. (13).

2.10 Eddy covariance (EC) reference flux measurements

2.10.1 EC instrumentation

A conventional EC system was set up for flux measurements of CO₂, sensible heat, and latent heat. The EC setup served as a reference for the TEA flux measurements. Instruments used were a three-dimensional sonic anemometer of type R3 (Gill Instruments Ltd, UK) and an infrared gas analyzer, type LI-7500 (Li-COR Env. Inc., USA). Wind and mole fraction data were recorded at a 20 Hz frequency at a height of 2.5 m a.g.l. The EC sensors were mounted side by side to a separate sonic anemometer and air inlet used for TEA. The two sonic anemometers were separated by a distance of 1 m. For quality assurance, in addition to the above primary eddy covariance setup, we used the sonic anemometer of the eddy covariance in combination with the open-path fast response gas analyzer (IRGA) for an alternative eddy covariance flux estimate. The horizontal separation between the EC sonic and the IRGA was 0.35 m, and the separation between the TEA sonic and the IRGA was 0.7 m.

2.10.2 EC flux computations

Eddy covariance raw data were post-processed to obtain fluxes at a resolution of 30 min using the EddyPro^® software (LI-COR Env. Inc., USA), version 5.0.0. The flux processing comprised the following steps:

• statistical tests for raw data screening according to Vickers and Mahrt (1997), including spike count and removal, amplitude resolution, dropouts, absolute limits, skewness, and kurtosis;

• de-trending of raw time series by block averaging;

• compensation of time lag between sonic anemometer and gas analyzer measurements by covariance maximization;

• axis rotation for tilt correction using the planar fit method (Wilczak et al., 2001) with removal of the velocity bias (Van Dijk et al., 2004), in running window mode (1 d window, updated every 30 min) (TEA sonic) and fixed period (7 d) mode (EC sonic);

• flux quality check according to Foken et al. (2004), selecting classes 0 and 1 for further analysis on a scale of 0, 1, and 2.

3 Results and discussion

This section is organized as follows: (1) meteorological conditions during the experiment are presented, followed by (2) a characterization of mass flow control performance, a prerequisite for the (3) determination of concentration differences between accumulated updrafts and downdrafts, which, in combination with (4) vertical wind measurements, finally result in (5) trace gas fluxes. To inform the discussion on uncertainties of the eddy accumulation method, (6) coordinate rotation results, (7) uncertainties of vertical wind distributions, and (8) instrumental errors of the sonic anemometers and infrared gas analyzers used are presented.

Figure 7 Meteorological conditions and turbulent energy fluxes during the experimental period from 4 to 10 April 2015: photosynthetically active radiation (PAR), air temperature, relative humidity, soil temperature, wind velocity, wind direction, sensible heat flux (H), latent heat flux (LE), and friction velocity. The wind and eddy covariance flux observations shown in panels (e)–(h) were obtained from the R3 sonic of the eddy accumulation system in combination with the LI-7500 gas analyzer.

3.1 Meteorological conditions

Meteorological conditions (Fig. 7) during the experimental period from 4 to 10 April 2015 were characterized by fair weather conditions with photosynthetically active radiation peaking at around 1500 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ at noon (a). Air temperature was initially below 10 ^∘C on the first day with frost during the night but then rapidly increased to more than 20 ^∘C on the last day (b). Similarly, there was a positive soil temperature trend (d). No precipitation was observed during this period. Wind direction (f) was dominated by easterly winds except for 7 and 8 April with mostly westerly and stronger winds (e) and higher relative humidity (c). 7 April stands out as the day with highest wind speed (e), high friction velocity (h), low radiation (a), and low sensible heat flux (g). The sensible and latent heat fluxes (g) otherwise largely tracked radiation levels with the highest latent heat fluxes observed on 9 and 10 April and with a decreasing Bowen ratio.

Figure 8 Linearity of two digital mass flow controllers verified by a conventional thermal mass flow controller for a series of 100 constant flow rate levels in the range of 0 to 3000 sml min⁻¹. The thermal mass flow controller reading is shown on the y axis versus the set point of the digital mass flow controllers on the x axis (black dots). The set point range on the x axis from 0 to 1.0 is linearly related to a flow rate range of 0 to 3000 sml min⁻¹. Observed linearity errors of the digital mass flow controllers shown here are below 0.3 % of the reading. Linear model fit of the digital mass flow controller of the “updraft channel” of the TEA system (dashed red line), and of the “downdraft channel” (dashed blue line), respectively.

Figure 9 Mass flow rate readings of a conventional thermal mass flow controller (Vögtlin red-y smart series), in blue, measuring a dynamically changing reference flow (black). The magnitude of the reference flow rate, which is proportional to the magnitude of recorded vertical wind velocity data, varied at a frequency of 10 Hz. The reference flow (black) was generated using the new digital mass flow controller. The uncertainty range of the digital mass flow controller (red) is small. Conversely, the mass flow reading of the conventional thermal mass flow controller deviates significantly from the reference as the conventional mass flow controller is unable to follow the highly dynamic reference signal.

3.2 Mass flow controller performance

Laboratory tests of the new digital mass flow controller used in this study showed that the accuracy and linearity of the digital mass flow controller for stationary flow were excellent: the maximum deviation of the new design from the conventional thermal mass flow controller used as reference was 0.3 % over the full operating range (Fig. 8). The minimum and maximum flow rates of the digital mass flow controller of 0.025 and 3 sml min⁻¹ correspond to a dynamic range of 120, exceeding minimum performance requirements formulated by Hicks and McMillen (1984) for eddy accumulation, i.e., a dynamic range of 100 or higher. The tests also showed that the two controllers used for the updraft and downdraft channels of the TEA system, respectively, performed to the same level (red and blue line in Fig. 8).

We further compared the performance of the two mass flow controller designs under dynamic flow rate conditions (Fig. 9). The new digital mass flow controller was used to generate a highly dynamic reference signal for comparison with the conventional thermal mass flow controller. The uncertainty of the new digital mass flow controller (red uncertainty ranges in Fig. 9) was small relative to the dynamic signal itself. Conversely, the conventional thermal mass flow controller (blue line in Fig. 9) was unable to follow the dynamic reference signal, showing strong overshoot and overestimation of the reference flow rate.

Figure 10 Regression of (a) the thermal mass flow controller flow rate versus the reference flow rate and (b) the digital mass flow controller flow rate with errors versus the reference flow rate. The reference flow rate was varied at a frequency of 10 Hz proportional to measured vertical wind velocity data. The black solid line denotes the 1:1 line, and the dashed red line denotes the linear model fit. Regression statistics: linear model equation; coefficient of determination, R²; and root mean square error, RMSE.

Statistics of the mass flow controller comparison for the same data, as in the previous figure, are shown in the regression plots in Fig. 10. The new digital mass flow controller showed a perfect slope of 1.00, a high coefficient of determination of R²=0.999, and a root mean square error (RMSE) of 8.8 sml min⁻¹ (Fig. 10b), which is 50 times smaller than the RMSE of the conventional mass flow controller of 441.8 sml min⁻¹ (Fig. 10a). The conventional mass flow controller further overestimated the reference by 300 % and matched less than half of the variance of the reference signal (R²=0.43).

In conclusion, the dynamic response, precision, and accuracy of the digital mass flow controller are suitable for eddy accumulation sampling, while the conventional thermal mass flow controller is not.

Figure 11 Dry molar fraction of CO₂ in the updraft and downdraft reservoirs of the true eddy accumulation device, in red and blue, respectively. The difference in molar fraction is a result of the vertical CO₂ flux.

3.3 CO₂ molar fraction and differences between accumulated updrafts and downdrafts

Time series of the CO₂ mole fraction obtained by the TEA system separately for the accumulated updrafts and downdrafts are shown in Fig. 11. Both the accumulated updrafts and downdrafts followed the common diurnal pattern of the CO₂ mole fraction with minimal CO₂ densities during the day when photosynthetic activity of the vegetation is at its maximum, a gradual buildup of CO₂ from the late afternoon through the night, and finally a rapid decrease of CO₂ in the morning when the daytime turbulence removes nightly accumulation of trace gases and photosynthesis then further draws down the ambient CO₂ mole fraction. As expected, despite the generally similar course of the CO₂ mole fraction of the updraft and downdraft reservoirs, there was a small but systematic difference between the two, with the CO₂ mole fraction of the updrafts (red line in Fig. 11) being lower than the downdrafts (blue line in Fig. 11). This difference was caused by the relative CO₂ depletion of updraft air due to photosynthesis during the day. At night, the inverse pattern was observed whereby updraft air was systematically enriched in CO₂ through respiration from soil and vegetation.

Figure 12 Difference in dry mole fractions of CO₂ between the updraft and downdraft reservoirs of the true eddy accumulation device. The 30 min raw time series shown is the result of subtracting the mole fraction in the downdraft reservoir from the mole fraction in the updraft reservoir (Fig. 11). A positive CO₂ mole fraction difference indicates a CO₂ flux away from the surface (respiration), and a negative CO₂ mole fraction difference indicates a CO₂ flux towards the surface (assimilation).

The difference in the CO₂ mole fraction between updraft and downdraft reservoirs is shown in Fig. 12. This difference was positive during the night and negative during the day. Windy conditions during the day cause a smaller magnitude of CO₂ difference as seen on 7 April (see Fig. 7 for wind and Fig. 12 for CO₂). Likewise, calm conditions enhance the CO₂ difference between updraft and downdraft reservoirs (see 9 and 10 April 2015, Fig. 12).

Figure 13 Vertical wind velocity of updrafts and downdrafts, in red and blue, respectively, averaged to 30 min resolution and shown as absolute values. Vertical velocity is subject to a 1 d running window real-time planar fit coordinate rotation as obtained from the TEA system. Mean vertical wind velocity at 30 min resolution after running window real-time planar fit coordinate rotation, in blue, and after a post-processing planar fit coordinate rotation with the fit period corresponding to the full period shown, in black. Note that updrafts (red) and downdrafts (blue) shown here are different from $\stackrel{\mathrm{‾}}{\left|w\right|}$ used for flux calculation according to Eq. (7).

3.4 Mean absolute vertical wind velocity

Vertical wind velocity measurements from the TEA system are shown in Fig. 13, separately for updrafts and downdrafts (red and blue lines, respectively). Both updrafts and downdrafts show similar magnitude, which is to be expected for a mean vertical wind velocity close to zero (black and cyan line in Fig. 13). On 9 and 10 April, absolute vertical wind velocity w during the day was lower than for other days (Fig. 13). Lower absolute w, indicating less vertical mixing, corresponded with more pronounced differences in the CO₂ molar fraction between updrafts and downdrafts, i.e., a more negative difference (Fig. 12) on the same 2 d. Under conditions of low winds and low turbulence, but intense radiation, the air close to the surface and the vegetation, which would be sensed as updrafts, was depleted in CO₂ through photosynthesis, relative to the air above.

Mean vertical wind velocity, $\stackrel{\mathrm{‾}}{w}$, which ideally is zero over the 30 min flux integration intervals, rarely exceeded ±4 cm s⁻¹ (dashed black lines), a threshold which, according to a simulation by Hicks and McMillen (1984), should not be exceeded to avoid significant flux errors. On two occasions, $\stackrel{\mathrm{‾}}{w}$ from the running window planar fit showed larger deviations from zero on 4 and 9 April 2015. Overall, the amplitude of $\stackrel{\mathrm{‾}}{w}$ was smaller for the planar fit rotation using a single 7 d window compared to the running planar fit using 1 d windows.

Figure 14 Comparison of CO₂ fluxes measured by true eddy accumulation (TEA) and eddy covariance (EC). Flux integration interval for both methods is 30 min. TEA fluxes were obtained every 60 min and EC fluxes every 30 min. The EC fluxes are quality-flag-filtered, accepting flags ≤1 on a scale of 0, 1, and 2 (Foken et al., 2004). The two alternative EC flux calculations shown, in blue and red, respectively, are from two separate Gill-R3 sonic anemometers using CO₂ density data from the same LI-7500 gas analyzer. The symbols indicate if EC flux estimates from the two sonic anemometers were within 50 % from each other (full circles) or not (crosses).

Figure 15 Mean diurnal cycle (a) and regression (b) of CO₂ fluxes measured by true eddy accumulation (TEA) and eddy covariance (EC) for the period 4–7 April 2015. Flux integration interval for both methods is 30 min. TEA fluxes were obtained every 60 min and EC fluxes every 30 min. The two alternative EC flux calculations shown, in blue and red, respectively, are from two separate Gill-R3 sonic anemometers using CO₂ density data from a single LI-7500 gas analyzer. The EC fluxes are quality-flag-filtered, accepting flags ≤1 on a scale of 0, 1, and 2 (Foken et al., 2004) and filtered for consistency of EC fluxes from the two anemometers, accepting EC fluxes with a maximum relative flux difference of 50 %. Error bars indicate ±1 standard deviation. (b) Scatterplot and linear model fit for TEA fluxes versus EC fluxes. The three EC flux versions shown are fluxes from the TEA sonic (in red), from the EC sonic (in blue), and the mean of the fluxes from the TEA and the EC sonic (in brown). Also shown are EC fluxes from EC sonic versus EC fluxes from TEA sonic (in black). Regression statistics shown are the coefficient of determination, R², and the linear model equation from a standard major axis regression (Legendre and Legendre, 1998; Sokal and Rohlf, 1995). The same color code (red, blue, brown, black) applies to all graphical elements, i.e., points, regression lines (solid lines), and regression coefficients. The dashed black line is the 1:1 line.

3.5 CO₂ fluxes

The turbulent exchange of CO₂ between the vegetation and the atmosphere as observed by TEA is displayed in Fig. 14. Fluxes clearly show CO₂ uptake during the day (photosynthesis, with values up to ca. −20 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and CO₂ release during the night (respiration, with values up to ca. +5 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$); see also the mean diurnal cycle in Fig. 15a. The temporal variability of CO₂ fluxes measured by TEA (Fig. 14) reflects variation in photosynthetically active radiation (PAR; Fig. 7): 8, 9, and 10 April 2015 with high levels of radiation and the least amount of cloud also show the highest fluxes. 6 April, which experienced more clouds during the day and therefore less abundant photosynthetically active radiation, also showed relatively small CO₂ fluxes during the day. Similarly, on 5 and 7 April, which were affected by clouds and reduced radiation in the early afternoon, the CO₂ fluxes during the afternoon were reduced compared to cloud-free days, e.g., 9 April.

Nighttime CO₂ fluxes measured by TEA showed a trend of increasing fluxes over the experimental period (Fig. 14). This trend corresponded to the observed trends of increasing air temperature and soil temperature (Fig. 7) over the same period. The observation of a positive correlation of positive CO₂ fluxes (respiration) with soil and air temperature is in line with the widely accepted mechanistic understanding that soil respiration is a function of soil temperature with a positive correlation over the temperature range presented. Overall, from an ecophysiological point of view, the observed CO₂ fluxes corresponded well with their meteorological drivers.

From a methodological point of view, a key step in assessing the TEA method's performance is the comparison of the TEA method with the established EC method. Fig. 14 shows CO₂ fluxes measured by both TEA (black line) and EC (red and blue points). The eddy covariance fluxes shown in Fig. 14 are quality-filtered, accepting flags ≤1 (Foken et al., 2004; on a scale of 0, 1, and 2). This filtering reduced the number of available 30 min eddy covariance flux estimates to 90 % and 93 %, and filtering for flag =0 reduced the fraction of available data to 45 % and 46 % for the sonic anemometers of the TEA and the EC system, respectively (both used for eddy covariance fluxes).

Generally, good agreement was observed between the TEA and EC methods. However, the intercomparison of the methods was complicated by the presence of high noise levels in the 30 min EC flux estimates. Analysis of R3 sonic anemometer raw data revealed that the EC sonic anemometer (blue line Fig. 14) was affected by correlated noise in the high-frequency wind and sonic temperature measurements, resulting in erroneously high sensible heat flux and momentum flux estimates. The erroneously high variance of the horizontal and vertical wind components was particularly pronounced during the night and decreased over the experimental period. Conversely, erroneously high sonic temperature variance increased over the experimental period. CO₂ fluxes were affected to a lesser degree because the noise in the R3 sonic anemometer data was not necessarily correlated with the LI-7500 infrared gas analyzer measurements.

The sonic anemometer of the TEA system seemed less affected by noise. Therefore we used it as an alternative input to the eddy covariance flux computations. However, the horizontal separation between the latter sonic anemometer and the open-path gas analyzer was large (0.7 m) relative to the low measurement height. Eddy covariance CO₂ fluxes from the TEA sonic anemometer were on average 18 % smaller compared to the EC sonic anemometer. However, after eliminating inconsistent fluxes whereby the fluxes from the two sonic anemometers disagreed by more than 50 % or alternatively 25 %, eddy covariance CO₂ fluxes from the TEA sonic anemometer were on average 6 % smaller and 2 % larger, respectively, compared to the EC sonic anemometer (Fig. 15).

A regression of CO₂ fluxes from the TEA method versus the EC method is shown in Fig. 15b, comparing TEA fluxes to the two alternative EC flux estimates after quality filtering eddy covariance fluxes with flags ≤1 (Foken et al., 2004). Before applying the above consistency filter, the coefficients of determination were R²=76 % and 67 %, and the linear model slopes were 1.04 and 0.87 when quality filtering eddy covariance fluxes for flags ≤1 (Foken et al., 2004) and improved to R²=71 % and 76 % and slopes of 0.99 and 0.92 for flags =0, for (i) TEA fluxes versus EC fluxes (using the TEA sonic anemometer for EC) and (ii) TEA versus EC fluxes (using the EC sonic anemometer for EC), respectively. The correlation further improved when rejecting conflicting EC fluxes: limiting the difference of the two EC flux estimates from the two independent sonic anemometers to 50 % of the mean of the two fluxes (flags ≤1), the coefficients of determination of TEA versus EC fluxes increased further to 80 %, 80 %, and 82 %, with linear model slopes of 1.01, 0.96, and 0.99, for (i) TEA fluxes versus EC fluxes (using the TEA sonic anemometers for EC), (ii) TEA versus EC fluxes (using the EC sonic anemometer for EC), and (iii) TEA versus EC fluxes (using the average of the TEA and EC sonic anemometers for EC), respectively (case shown in Fig. 15b).

While a certain fraction of the observed deviations between TEA and EC flux estimates can be attributed to the methodological differences, the two independent eddy covariance flux estimates also showed deficiencies in agreement, with coefficients of determination of the EC to EC regression of R²=84 %, 86 %, and 85 % and slopes of 0.81, 0.82, 0.85, leaving 16 %, 14 %, and 15 % of the flux variance unexplained for quality flag filter thresholds of 2, 1, and 0. The percentage of additional unexplained variance in the TEA fluxes with R²=74 % (regression of TEA fluxes versus mean EC fluxes from the two sonic anemometers, flag ≤1) relative to EC fluxes with R²=86 % was 12 % only. There was a positive intercept for the above regression cases between TEA and EC fluxes which ranged from 1.8 to 2.4 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$.

While the above-reported coefficients of determination for the EC versus EC flux regressions using various filter options were somewhat higher than for TEA versus EC fluxes, the differences between the TEA and EC methods were still of the same order of magnitude as the differences between the two EC flux estimates, which use the same infrared gas analyzer. It might seem from Fig. 15a and the above-reported statistics that the two EC flux estimates match more closely than TEA and EC fluxes. However, we note that the EC–EC agreement is to some degree artificial; i.e., only a part of this apparent agreement relates to differences between TEA and EC methods, and another part is simply due to the “consistency filter” applied to the two EC flux estimates. This filter, which by definition discards EC flux values which are not similar, is not normally required but was introduced to deal with the above-mentioned deficiencies of the two EC setups.

We would expect from a side-by-side comparison of eddy covariance flux measurements, using identical models of research class sonic anemometers and sharing the same gas analyzer, that the R² would exceed 90 %. We interpret the compromised match of the two current EC estimates largely as a result of compromised wind and sonic temperature measurements by the two R3 anemometers. When excluding the first 4 d which were relatively more affected by erroneous w, filtering for quality flags ≤1 and again filtering EC fluxes to not exceed a relative difference between the two EC flux estimates of 50 %, then the match of TEA versus EC CO₂ fluxes improved further, yielding R² values of 84 %, 86 %, and 86 %, slopes of 1.04, 0.91, and 0.98, and intercepts of 2.0, 1.9, and 2.0 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for TEA fluxes versus EC fluxes using the TEA sonic anemometers, the EC sonic anemometers, and the mean of the TEA and EC sonic anemometers, respectively.

In relation to previous works on true eddy accumulation trace gas flux measurements, we note that, despite compromised data quality of one of the R3 sonic anemometers of the current study, the match of true eddy accumulation and eddy covariance CO₂ fluxes exceeded the match in any previously published true eddy accumulation experiments we are aware of. The closest results are those from Komori et al. (2004), who obtained a coefficient of determination for TEA versus EC CO₂ fluxes of R²=0.64, a slope of 0.95, and a relatively high intercept of 8.6 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for 17 flux integration intervals.

3.6 Uncertainty of vertical wind measurements

Vertical wind measurements contribute to flux uncertainty in two ways: (i) instrumental errors of the sonic anemometer and (ii) nonideal wind field and nonzero mean vertical wind velocity over the flux integration period.

3.6.1 Sonic anemometer measurement errors

Systematic and random errors in sonic anemometer measurements contribute to scalar flux uncertainty both in the EC and TEA methods. For a detailed analysis of sonic anemometer measurement errors we refer to instrument comparisons and quantifications of measurement errors for common sonic anemometer types (Loescher et al., 2005; Mauder and Zeeman, 2018; Foken et al., 2019, and others). These studies suggest that measurement errors of sonic anemometers, including differences between different types of sonic anemometers and differences between different units of the same type of sonic anemometer, may account for anywhere from several percent up to about 25 % of the error in scalar flux measurements.

Regarding the R3-type sonic anemometer used in the current study, Loescher et al. (2005, Fig. 5) found in wind tunnel tests that the R3 sonic anemometer, like other post-mounted designs, suffered from flow distortion, systematically overestimating vertical wind velocity. The R3 overestimated vertical wind velocity for vertical velocities below 0.15 m s⁻¹ by up to ca. 0.05 m s⁻¹ for vertical velocities close to zero and underestimated vertical velocity by up to the same amount for vertical velocities up to 0.3 m s⁻¹ (gain error). In addition, when the stanchions supporting the upper transducers were in the flow path, the vertical wind velocity response was nonlinear. Nonlinearity and gain errors can result in misalignment of the coordinate system with the mean streamlines (see Sect. 3.6.2) and apparent asymmetry of vertical wind distributions (see Sect. 3.6.3).

Figure 16 30 min means of vertical wind velocity, w, before (a, c) and after (b, d) coordinate rotation as a function of wind direction (a, b) and as a function of along-slope wind, U (c, d). The vertical wind velocity w in (a) and (b) is normalized by horizontal wind velocity. Panels (b) and (d) both differentiate between two approaches to the planar-fit coordinate rotation: the real-time planar fit applied in the TEA measurements, deploying a moving window of 1 d (“running PF”, in red and orange), and a single planar fit, covering the full measurement period (7 d long) obtained in post-processing (“fix PF”, in black and gray). The labels “rot” and “unrot” denote the rotated and unrotated vertical wind velocities, respectively.

3.6.2 Coordinate rotation

The following analyses present nonzero mean vertical wind velocities, which cause scalar flux uncertainty. Figure 16 shows vertical wind before and after coordinate rotation. Vertical wind velocity scaled by horizontal wind velocity ideally follows a (co)sine function when anemometer coordinates are tilted relative to streamline coordinates (Fig. 16a). Planar fit rotation reduced the range of $\stackrel{\mathrm{‾}}{w}$ from about ±0.1 m s⁻¹ to about ±0.05 m s⁻¹ (Fig. 16b). The range of $\stackrel{\mathrm{‾}}{w}$ was slightly smaller for the 7 d rotation window compared to the running 1 d window. Figure 16c indicates that before coordinate rotation the sonic coordinate frame was tilted relative to the streamlines, which followed the terrain slope. The slope of the relation of $\stackrel{\mathrm{‾}}{w}$ over along-slope wind U (c) vanished after planar fit rotation using the 7 d window (d). The running planar fit with a 1 d window did not fully remove the dependence of $\stackrel{\mathrm{‾}}{w}$ from along-slope wind velocity, biasing $\stackrel{\mathrm{‾}}{w}$ by up to ±0.02 m s⁻¹ over the range of along-slope velocities shown (d), corresponding to about 15 % of bias before tilt correction.

Figure 17 30 min mean vertical wind velocity w as a function of horizontal wind velocities u and v before and after planar fit coordinate rotation for two different planar fit procedures and two periods. Wind velocities were obtained by the Gill-R3 sonic anemometer which was used for true eddy accumulation. (a) Unrotated vertical wind velocity. (b) Rotated vertical wind velocity from real-time 1 d moving window planar fit rotation performed by TEA system. (c) Rotated vertical wind velocity with single planar fit rotation period performed in post-processing. (d) Same as (c) but for a longer period, i.e., from 1 April to 31 May 2015. Panels (a) to (c) show data from the experimental period of the current study from 4 to 11 April 2015. Black dots indicate the location of individual 30 min mean vertical wind velocity readings in the u–v velocity space. Red arrows indicate the direction to north. Dashed lines, in white, indicate the azimuth of the vertical stanchions of the sonic head structure relative to the center of the sonic coordinate system.

Figure 17 presents residuals of mean vertical wind velocity $\stackrel{\mathrm{‾}}{w}$ before and after coordinate rotation in u and v horizontal velocity space. This analysis relates $\stackrel{\mathrm{‾}}{w}$ to wind direction, horizontal wind velocity, and obstacles causing flow distortion. $\stackrel{\mathrm{‾}}{w}$ (in Fig. 17a) before tilt correction shows the effect of the terrain slope on wind measurements in sonic coordinates. The slope effect was fully removed through planar fit rotation (Fig. 17b). The 1 d running window planar fit led to a slight overcorrection of the tilt (Fig. 17b), as already noted regarding Fig. 16d, which was not the case for the planar fit rotation using a 7 d window. Residuals of $\stackrel{\mathrm{‾}}{w}$ were small with the 1 d running planar fit, resulting in relatively larger $\stackrel{\mathrm{‾}}{w}$ residuals up to about ±0.05 m s⁻¹ (Fig. 17b) compared to up to about ±0.02 m s⁻¹ for the 7 d planar fit (Fig. 17c).

Some dependence of $\stackrel{\mathrm{‾}}{w}$ on horizontal wind velocity and direction was observed: for horizontal wind velocities of more than about 1 m s⁻¹, residuals of $\stackrel{\mathrm{‾}}{w}$ were mostly positive, particularly for southwesterly and northeasterly winds (Fig. 17b and c). This was confirmed by the planar fit rotation of a 2-month-long data set (Fig. 17d). Possible interpretations include flow distortion from trees and bushes in the southwest and northeast and velocity-dependent flow distortion of the sonic anemometer or nearby structures. No obvious influence of the stanchions of the anemometer on $\stackrel{\mathrm{‾}}{w}$ was identified (Fig. 17a–d).

3.6.3 Updraft–downdraft asymmetries

The TEA method requires $\stackrel{\mathrm{‾}}{w}=\mathrm{0}$, which would result from symmetry in updraft and downdraft statistics. However, we observed asymmetry in the mean, the count, and the sum of updraft and downdraft samples. Quantification of the flux uncertainty due to asymmetric distributions of w would require co-spectral information of w and CO₂ densities, which is generally absent for TEA measurements. Instead, we present a quantification of observed asymmetries of w, informing about the magnitude of the asymmetries and their variability over time.

Figure 18 Ratio of statistical measures of vertical wind velocity of updrafts and downdrafts per 30 min flux integration interval. The statistical measures are the ratio of the count, the mean, and the sum of vertical wind velocity records during updrafts and downdrafts, respectively (solid lines). The dashed lines indicate the temporal mean of the above statistics over the period displayed.

Figure 18 shows that (i) on average, the count of updrafts was larger than the count of downdrafts; (ii) however, on average, the mean of updrafts was smaller than the mean of downdrafts; and (iii) on average, the sum of updrafts was slightly smaller than the sum of downdrafts. It is noteworthy that the mean of w can still be zero while the 7 d mean of the 30 min mean of updrafts and the 30 min mean of downdrafts are nonzero; or in other words, the ratio of the 30 min mean of the updrafts to the 30 min mean of the downdrafts is different from unity as observed here. This can be understood by considering the different weights of updrafts and downdrafts in the mean which reflect the asymmetry in the counts of updrafts and downdrafts, respectively.

Figure 19 Probability density distributions of the relative uncertainty (RU) of the mean absolute vertical wind velocity, $\stackrel{\mathrm{‾}}{\left|w\right|}$, due to nonzero mean vertical wind velocity over the 30 min flux integration periods (blue solid line). $\stackrel{\mathrm{‾}}{\left|w\right|}$ is required for the determination of trace gas fluxes according to Eq. (2). The relative uncertainty of the mean vertical wind velocity was calculated as $\mathrm{RU}=\left(\right(\stackrel{\mathrm{‾}}{w_{\mathrm{up}}}-\stackrel{\mathrm{‾}}{\left|w_{\mathrm{down}}\right|})/\mathrm{2})/\left(\right(\stackrel{\mathrm{‾}}{w_{\mathrm{up}}}+\stackrel{\mathrm{‾}}{\left|w_{\mathrm{down}}\right|})/\mathrm{2})$, where the overline denotes the temporal mean of vertical wind velocity over the 30 min flux integration interval, and subscripts “up” and “down” refer to updrafts and downdrafts, respectively. Results based on w from the real-time moving window planar fit coordinate rotation of the TEA system and from the fixed window post-processing planar fit rotation (FPF) are shown in blue and red, respectively. “Mean RU of means”, “Mode RU of means”, and “Quant. RU of means” indicate the mean, the mode, and the quantiles, respectively, for a probability of 10 % and 90 %. While mean vertical wind velocity $\stackrel{\mathrm{‾}}{\left|w\right|}$ is needed for flux derivation, the sum of vertical wind velocity over the flux integration interval, $\sum \left|w\right|$, relates to the accumulated air sample volumes. The relative uncertainty of the sums, “RU of sums”, corresponds to the contribution of vertical wind velocity to the volume mismatch correction defined in Sect. 2. Relative uncertainty of the sum of vertical wind velocity per 30 min flux integration interval was calculated as $\mathrm{RU}=\left(\right(\sum w_{\mathrm{up}}-\sum \left|w_{\mathrm{down}}\right|)/\mathrm{2})/\left(\right(\sum w_{\mathrm{up}}+\sum \left|w_{\mathrm{down}}\right|)/\mathrm{2})$. “Mean RU of sums”, “Mode RU of sums”, and “Quant. RU of sums” indicate the mean, the mode, and the quantiles of the distribution, respectively, for 10 % and 90 % probability.

Figure 19 shows probability density distributions of the relative uncertainty of the mean absolute vertical wind as defined in the figure caption. The distribution of the sums of updrafts and downdrafts is centered around zero, with only a small negative bias of −0.01 for the mode of the distribution and a larger bias of −0.04 for the mean of the distribution. Conversely, the distribution of the relative uncertainty of the means peaked for both the mode and the mean of the distribution at a more negative bias of −0.11, as observed for the running window planar fit (Fig. 19, blue line). Similar results were obtained for the 7 d stationary planar fit although with a smaller negative bias of −0.08 (Fig. 19, red line). Less than 10 % percent of relative uncertainty values were more negative than −0.3, and less than 10 % of relative uncertainties were larger than +0.14. In summary, updraft–downdraft asymmetries were on the order of 10 % of the mean absolute vertical wind velocities used in flux calculations.

3.7 Uncertainty of trace gas concentration measurements

Trace gas flux errors are a function of accuracy and precision of the trace gas analysis. Regarding accuracy, bias between the two infrared gas analyzers used for the TEA and EC methods is more relevant than absolute accuracy for comparing the TEA method to the EC reference method. By comparing time-averaged time series of CO₂ concentrations of the LI-6262 and LI-7500 infrared gas analyzers, a time-variable bias was found, which accounted for up to 5 % of the scalar flux.

Regarding scalar flux errors from the TEA method, systematic errors leading to bias between measurements of updrafts and downdrafts and precision are important. Systematic errors biasing the concentration difference between updrafts and downdrafts are difficult to quantify. The following results quantify the precision of the gas analysis, based on the analysis of four replicated measurements of 120 s each of the accumulated updraft and downdraft concentrations, per 30 min flux integration interval. The results comprise precision of the gas analyzer at the measurement frequency of 1 Hz as well as precision of the TEA gas sampling storage and delivery system, feeding samples to the gas analyzer. The latter includes drift of the gas analyzer signal and of the trace gas concentration over the time required to determine a concentration difference between updrafts and downdrafts, i.e., 2 times 150 s.

Regarding the observed CO₂ concentration signal, 90 % of CO₂ dry molar fraction differences between updrafts and downdrafts at 30 min integration were between 0.14 and 2.08 µmol mol⁻¹ (Fig. 20). Regarding the observed precision of the total gas analysis system under field conditions, for 90 % of flux integration intervals, the CO₂ dry molar fraction measurements over four replicated measurements of the updraft reservoir varied in the range of 0.033 to 0.48 µmol mol⁻¹. The precision of downdraft measurements was 50 % lower, with 90 % of the downdraft measurements showing a range of the four replicates of 0.05 and 0.73 µmol mol⁻¹.

Figure 20 Density distributions and selected probabilities of the trace gas concentration difference signal and of the measurement uncertainties, and gas analyzer precision. The signal is the difference in molar fraction between accumulated updrafts and downdrafts (black line). The uncertainty of the measurements is expressed as the range of four replicated samples of the updraft and downdraft concentration, in red and blue, respectively. Vertical dashed lines indicate various probabilities of the above distributions. Vertical solid lines indicate the nominal precision of the LI-6262 (LI-COR, USA) infrared gas analyzer used in this study, calculated for an integration time of 120 s (cyan). Precision of two types of laser spectrometers are also shown for reference: Allan deviation of FGGA (Los Gatos Research, USA), in magenta, and G2301 (Picarro, USA), in orange, both determined in the laboratory.

For 85 % of the flux integration intervals, the signal, i.e., the dry molar fraction difference between updrafts and downdrafts, was larger than the nominal and extrapolated 120 s precision of the LI-6262 infrared gas analyzer used in this study, as well as the precision, i.e., Allan deviation, of two laser spectrometers we tested in the laboratory. The latter two instruments were not used in the current study, but characteristics are provided to put the instrument used in this study in the context of current state-of-the-art greenhouse gas monitors. Note that the indicated precision of the LI-6262 of 0.29 µmol mol⁻¹ is an extrapolation of nominal precision and drift values to 120 s, for which nominal precision was given as peak-to-peak noise, rather than Allan deviation, which was used to characterize the laser spectrometers and is by definition smaller than or equal to the peak-to-peak noise.

Also note that the Allan deviation at 120 s integration of the G2301 (Picarro) instrument of 0.0125 µmol mol⁻¹ appeared to be only 7.4 % of the Allan deviation of the FGGA (Los Gatos Research) laser spectrometer of 0.17 µmol mol⁻¹. However, in addition to differences in the design of the two spectrometers and any potential differences in test conditions, it appears from the analysis of Allan deviation that the G2301 (Picarro) instrument may be subject to some degree of internal smoothing of the gas concentration readings. We cannot say with certainty to which degree such potential smoothing might have affected the Allan deviation at 120 s integration time because the manufacturer of the instrument was unable to provide further information on the suspected filter beyond acknowledging its existence.

The observed signal-to-noise ratio (SNR) of the total trace gas analysis system under field conditions ranged between −9 and +21 (and one value at −66.7) and improved with the magnitude of the signal itself (Fig. 21a). The slope and intercept of an ordinary least-squares linear model fit to this relation were 3.2 µmol⁻¹ mol and 0.48, respectively. The fact that this quasilinear relationship and the slope significantly differed from zero means that larger fluxes have relatively smaller errors of the type considered, a feature reducing absolute uncertainty of trace gas flux budgets.

Figure 21 (a) Signal-to-noise ratio (SNR) of the trace gas dry molar fraction difference measurements as a function of the difference, the latter being proportional to the trace gas flux. The solid line is a linear model fit using ordinary least-squares regression. (b) Density distribution of the signal-to-noise ratio. Dashed lines in (a) and (b) indicate a signal-to-noise ratio of ±1, corresponding to probabilities of 37.7 % and 49.5 %, respectively. Consequently, in 88 % of the cases the signal-to-noise ratio is higher than 1, leaving 12 % of the cases unresolved. The display in (a) excludes one extreme value at SNR $=-\mathrm{66.7}$.

For 88 % of 30 min flux integration intervals, the signal-to-noise ratio was larger than 1 (Fig. 21b). Over the period of the experiment, the sum of the same noise data as above accounted for up to 25 % of the sum of the trace gas concentration signal, i.e., of the difference in CO₂ dry molar fraction between updraft and downdraft reservoirs. The percentage of 25 % is a maximum estimate for this type of noise, as it was determined from the range of the four replicates of concentration difference measurements, which is sensitive to extremes.

Figure 22 Range of CO₂ trace gas flux estimates observed by the true eddy accumulation method, accounting for the range of possible flux estimates from four replicated measurements of the dry molar fraction difference between accumulated updrafts and downdrafts per 30 min flux integration interval.

3.8 Uncertainty of trace gas flux measurements

The uncertainty of trace gas fluxes due to the uncertainty of the gas analysis is shown in Fig. 22. Over the period of the experiment, the sum of the noise range, i.e., the absolute value of the difference between the largest and smallest flux estimate, accounted for 37 % of the sum of the signal, i.e., the absolute value of the mean flux. As stated above for CO₂ concentrations, the percentage of 37 % is a maximum estimate of this type of noise because it was determined from the range of the four replicates of concentration difference measurements, which is sensitive to extremes, and because this estimate is additive. In practice, it is highly unlikely that this uncertainty range leads to additive errors; instead, some of the errors would cancel, leading to much smaller actual uncertainties.

Regarding the time series of trace gas flux noise, i.e., the range of maximum and minimum flux estimates (Fig. 22), calm conditions with low wind speeds and low friction velocities, e.g., on 10 April, result in relatively large concentration differences and relatively small vertical wind terms contributing to the trace gas flux calculations and therefore result in relatively low uncertainty of the flux due to the uncertainty of the gas analysis. The opposite can be observed for windy conditions with high friction velocity, e.g., on 7 April, which result in relatively small concentration differences and a relatively high contribution of the vertical wind term to the flux calculations and therefore a relatively high uncertainty of the flux due to the uncertainty of the gas analysis.

The total uncertainty of the trace gas flux needs to account for the uncertainty of the mass flow control, the uncertainty of the concentration differences, and the uncertainty of the vertical wind signal. Ideally, such analysis would incorporate the effect of different approaches of coordinate rotation not just on the residuals of $\stackrel{\mathrm{‾}}{w}$ but also the effect on the fluxes themselves. This would require consideration of co-spectral information of wind and scalar using time-resolved high-frequency data and simulations of true eddy accumulation with different coordinate rotation approaches.

4 Conclusions

The following conclusions intend to summarize the performance of the true eddy accumulation method, put the results of the current experiment into context relative to existing published studies, summarize and quantify main sources of uncertainty, report on limitations and lessons learned during the current experiment, suggest future improvements regarding technical and methodological aspects, and finally identify applications in which true eddy accumulation can facilitate novel flux measurements in the future.

The current study has presented CO₂ fluxes measured by true eddy accumulation. The TEA system measured continuously and automatically fluxes at 30 min resolution over a duration of more than 7 d. The TEA measurements were able to capture fluxes representing the biological activity of the system. TEA flux measurements compared favorably with eddy covariance reference measurements, with R² values of up to 86 % and a regression slope of 0.98.

A novel implementation of dynamic mass flow control was key to the success. It was 50 times more accurate in terms of root mean square error than the conventional thermal mass flow controller reference during laboratory tests and proved to be robust and without failure during more than 3 years of operating time in the field. Further innovative features were the digital signal processing and the real-time sampling decisions, incorporating online coordinate rotation and correction of the mean vertical wind, and, finally, the elimination of dead volumes in the gas sampling system.

Compared to earlier studies published on true eddy accumulation flux measurements (Desjardins, 1977, on temperature fluxes; Speer et al., 1985, and Neumann et al., 1989, on water vapor fluxes; Rinne et al., 2000, on isoprene fluxes; and Komori et al., 2004, on CO₂ fluxes), the current study obtained the best fit of TEA fluxes to EC fluxes of any trace gas or scalar. The current study also presents the longest continuous CO₂ flux measurements by TEA.

A detailed analysis of uncertainties of the TEA method was presented in terms of the uncertainty of the mass flow controllers, the uncertainty of the trace gas handling and analysis system, and the uncertainty of the vertical wind velocity measurements and 30 min means. Uncertainties of the eddy covariance method and instruments were partially quantified through two replicated flux computations using two alternative sonic anemometers. Uncertainties of the EC fluxes explained a significant fraction of the mismatch between the TEA and EC methods. The signal-to-noise ratio of the TEA trace gas analysis system allowed the concentration difference signal to be detected in 88 % of 30 min flux intervals. Maximum uncertainty estimates of the TEA trace gas measurement precision accounted for up to 25 % of the concentration differences and up to 37 % of the fluxes. A comparison of the precision of three gas analyzers suggests that deployment of state-of-the-art laser spectrometers would significantly reduce TEA flux uncertainty due to uncertainties in the gas analysis, with preliminary analysis suggesting an improvement in precision by a factor of 10 or more for some instrument models. This would likely reduce the flux uncertainty due to the gas analysis to about 5 % or less. Residual mean vertical wind velocities were generally smaller than 0.05 m s⁻¹. Uncertainties of the mean of absolute vertical wind velocities, which are needed for flux calculations, in terms of undesired residuals of mean vertical wind velocities after coordinate rotation were frequently on the order of 5 %. The uncertainties of the mass flow control were relatively small compared to uncertainties of the gas analysis, uncertainties of residual mean vertical wind velocities, and uncertainties of the eddy covariance flux estimates.

The following two design limitations were discovered: firstly, the continuous and long-term operation with frequent charging and discharging of the air sampling bags with on the order of 1500 charge cycles per month over time leads to increasing levels of fatigue of the material and in turn after a few weeks to a significant amount of leakage and therefore contamination of the samples with ambient air. The second observation relates to the intermittent nature of the gas flow, the variable accumulation volumes, and the intermittent gas analysis in the current bag-based accumulation design. Intermittent operation causes instationarity of the following parameters: air pressure in the gas handling system, temperatures of air and system components, and interactions of air constituents with the internal surfaces of the device such as adsorption and desorption of gas molecules on internal surfaces. Instationary conditions can lead to signal drift, variation of moisture content, and subsequently to less accurate flux measurements.

To address the above-mentioned limitations, we suggest exploring the idea of a new system design for TEA using rigid air containers of constant volume and with continuous-flow operation replacing flexible air bags. In such a new design the charging and discharging of the air reservoirs would happen continuously and at the same time. This new design principle would overcome the issue of material fatigue and compromised accuracy due to instationarities in the operation. A key methodological advantage of the new continuous-flow design is furthering the opportunity to merge the principles of true eddy accumulation sampling with eddy covariance sampling simultaneously with the very same measurement device, the same air samples, and the same gas analyzer.

Using a precise state-of-the-art laser spectrometer we have since implemented such a continuous-flow system suggested above and demonstrated its superior performance compared to conventional discrete cyclic charging of air bags. True eddy accumulation CO₂ fluxes observed with the new continuous-flow system were tightly correlated with eddy covariance fluxes, with R² values of up to 96 %. More details on the latter study will be reported separately.

The impact of coordinate rotation on true eddy accumulation fluxes has been discussed. We have suggested a new type of coordinate rotation, which we refer to as “surface fit”. Similar to the planar fit method, it aligns the coordinate system with the mean streamlines, accounting for a multidimensional parameter set including wind direction, flow distortion, and optionally other independent variables in an integrative, continuous way.

We would like to highlight the need for research on flux corrections for TEA in a comprehensive way similar to the body of work which exists on EC flux corrections. Future work needs to investigate and establish flux corrections specifically for the TEA method, including the equivalent to the correction of trace gas fluxes due to density fluctuations caused by simultaneous transfer of heat and water vapor (Webb et al., 1980). The derivation of this and other corrections specifically for the TEA method is nontrivial and will be addressed in separate work.

The current implementation of TEA suggests that this method has the potential to facilitate flux measurements of trace gases and other atmospheric constituents for which no fast gas analyzers are available. TEA is an alternative when the precision and accuracy of currently available analyzers are insufficient for high-frequency EC applications. The low power consumption of the current TEA systems with low sample flow rates will enable new applications, including off-grid use in solar- and battery-powered stationary and mobile applications. The long sample integration times give TEA a further advantage over EC, allowing for a simpler analyzer design compared to high-frequency analyzers at the same precision or alternatively providing ultimate precision through long integration times when using a high-quality analyzer.

It is evident that CO₂ fluxes in particular can be readily observed with alternative methods. However, the nonreactive and nonpolar trace gas CO₂ is an ideal candidate to assess the performance of the TEA method. The current experiment is a successful proof of concept demonstrating that true eddy accumulation with dynamic and accurate air sampling proportional to vertical wind velocity can be achieved in practice today. The lessons learned during the present work provide concise avenues including the above-outlined machine design considerations and required flux corrections for further improving the true eddy accumulation method to enable accurate and reliable flux measurements of more trace gases and atmospheric constituents than ever before.