Eddy covariance with slow-response greenhouse gas analysers on tall towers: bridging atmospheric and ecosystem greenhouse gas networks

Greenhouse gas monitoring is important to ensure climate goals are being achieved. This study unveils the potential of using atmospheric tall towers in direct flux measurements, bridging the gap between atmospheric and ecosystem monitoring networks. The ICOS Cities (PAUL) project aims to monitor CO₂ emissions in urban areas, where concentrated emissions make them key targets for climate change mitigation. This study explores the synergy between ICOS atmospheric and ecosystem networks by utilizing slow-response analysers (∼ 3 s) on tall atmospheric towers for ecosystem studies using the eddy covariance method. A standard setup with an ultrasonic anemometer and an infrared (IR) fast-response CO₂ analyser was installed and compared with measurements from an existing cavity ring-down spectroscopy (CRDS) analyser measuring CO₂, CO, and CH₄. Deployed on the 100 m Saclay tower near Paris, covering a 43.9 km² 80 % footprint with heavy traffic roads, a nearby heating plant, and a forest, the setup addressed technical challenges and height-induced complexities. Corrections for flux attenuation by high-frequency losses were limited to < 20 % on average for all stabilities and around 11 % for unstable conditions. Elevated mean fluxes for CO₂ (10 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and CH₄ (200 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) were observed from the heating plant wind direction during December and January. Conversely, the forest direction exhibited the strongest sink among all wind directions, with −4 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ during July and August. Storage and vertical advection were estimated using the routine three-level profile measurements done in ICOS atmospheric towers. Storage term was of the same magnitude as turbulent flux, increasing at night and de-stocking during the first half of the day. Vertical advection averaged zero on a monthly basis. These results demonstrate the feasibility and versatility of utilizing atmospheric towers for urban emission monitoring, offering valuable insights for emission monitoring strategies worldwide.

1 Introduction

Global surface temperature is 1.6 °C warmer on land compared with the pre-industrial era (IPCC, 2021), and projections show more than 2 °C warming in 2100 (IPCC, 2022). Warming results from the increase in greenhouse gas (GHG) concentration in the atmosphere, mainly driven by anthropogenic emissions (IPCC, 2021), of which 86 % comes from fossil fuel CO₂ (Canadell et al., 2023). Agriculture, forestry, and other land use (AFOLU) are a significant source of GHG (12.0 ± 2.9 Gt CO₂-eq yr⁻¹) while concurrently possessing the potential to remove CO₂ from the atmosphere (Jia et al., 2022).

Urban areas concentrate human activities and represent a significant source of GHG emissions, consequently making it one of the targets for mitigating climate change. Many northern countries' cities have ambitious GHG emission reduction plans over the next 2 decades that consist of electrifying the energy grid, implementing car-free zones, and investing in insulation improvement. Consequently, there arises an imperative for robust monitoring of urban areas' emissions reduction. Several works have tried to decompose eddy covariance measurements in (sub-)urban setups with different degrees of uncertainty (Velasco et al., 2009; Bergeron and Strachan, 2011; Ueyama and Takano, 2022). Currently, in Europe, the ICOS Cities (PAUL) project aims to advance technologies for monitoring CO₂ concentrations in urban areas of three different-sized pilot cities (Munich, Paris, and Zurich).

Monitoring GHGs in the atmosphere, ocean, and ecosystem is the objective of world-distributed research infrastructures such as ICOS in Europe (Heiskanen et al., 2022). For this purpose, different methods are used on terrestrial sites. Ecosystem sites focus on local flux monitoring using high-frequency measurements. These sites measure the different terms of the surface fluxes often representing a specific biome determined by the tower's footprint. In contrast, atmospheric towers precisely measure the concentrations as an imprint of larger-scale fluxes. These sites have a footprint spanning several hundreds of square kilometres and may be used to identify anomalies in CO₂ surface fluxes based on concentration (Ramonet et al., 2020) or retrieve surface flux by inverse modelling (Ciais et al., 2011).

At local scale, eddy covariance (EC) is the reference method for GHG monitoring. The method is praised for directly and continuously measuring surface turbulent flux and has largely been applied since early measurements of different gases, including water vapour, CO₂, CH₄, and N₂O (Valentini et al., 1996; Moncrieff et al., 1996; Fowler et al., 1995). Standard measurements require fast-response instruments, which is a technical limitation in measuring certain compounds' concentrations. Long-term measurement sites are equipped with CO₂ and H₂O gas analysers and in some wet or agricultural sites with N₂O or CH₄ analysers (Nemitz et al., 2018). At larger scales, atmospheric concentration measurements are often used alongside mesoscale–continental-scale transport models to solve surface flux (Lauvaux et al., 2012). This top-down approach is often validated locally by direct EC measurements (Vuichard et al., 2016).

The differences between typical atmospheric and flux tower monitoring setups are the following:

1. Atmospheric towers are generally taller, exceeding 100 m, while flux towers range from 2 to 40 m. This height difference is due to atmospheric measurements being designed to capture seasonal and annual trends in regional atmospheric background concentrations, necessitating the minimization of local source impacts (El Yazidi et al., 2018), while flux towers focus on those local sources.

2. Atmospheric towers feature more precise but slower measurements, not cadenced, with sampling intervals spanning several seconds, in contrast to the eddy covariance method used by flux towers, which samples at constant rates between 5 and 20 Hz.

3. Ecosystem stations include monitoring of vegetation and soil, which is absent in atmospheric stations, while atmospheric stations measure additional gaseous compounds such as CH₄, N₂O, and CO (Hazan et al., 2016).

Being able to use slow-response analysers to calculate flux by eddy covariance has been identified as a useful strategy to expand the flux networks to other compounds (Wohlfahrt et al., 2009). Atmospheric towers have high-precision analysers, which, if we can use them to compute eddy covariance fluxes, would provide multi-species flux measurements that would expand the flux network. This would require a fast 3D anemometer and continuous data logging at these sites. Not all towers may be suitable though, as flux towers have more constraints related to the topography and surrounding obstacles.

Moreover, in atmospheric towers, ICOS focuses on measuring not only CO₂ but also CO and CH₄ concentrations routinely. Therefore, measuring fluxes on these towers potentially enables the measurement of CO and CH₄ fluxes in the surrounding areas of each tower. These gases can bring relevant information on the site greenhouse gas budget, directly as CH₄ is a greenhouse gas and indirectly by revealing processes otherwise mixed in the CO₂ flux. Carbon monoxide (CO) emissions come in the most part from direct fuel combustion in the Northern Hemisphere mid-latitudes (Zheng et al., 2019). These emissions are likely to increase during winter due to reduced combustion efficiency at lower temperatures (Helfter et al., 2016). Methane (CH₄) emissions are primarily sourced from natural wetlands, enteric fermentation, and fossil fuels (Saunois et al., 2020). Locally, lakes may be a significant contributor (e.g. 0.1–0.3 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ in Iwata et al., 2020), but their contribution depends on the amount of sediments (Delwiche et al., 2021). Soils act primarily as sinks for both CO (Inman et al., 1971; Conrad and Seiler, 1980; Conrad, 1996) and CH₄ (Canadell et al., 2023; Dutaur and Verchot, 2007) through microbial oxidation processes.

In this study, we evaluate the capability of using an atmospheric monitoring tower with a slow-response analyser to compute turbulent fluxes of CO₂, CH₄, and CO by adding a fast-response sonic anemometer. For this purpose we installed a standard eddy covariance setup for CO₂ at 100 m at the ICOS FR-SAC atmospheric tower in the south of Paris, collecting 4 months of data from July 2023 until October 2023. The chosen site is located in a sub-urban region surrounded by a mix of agriculture, forest, wetlands, roads, and areas with buildings. We computed net CO₂ flux for slow- and fast-response analysers and compared them. The high-frequency loss was determined and the correction procedure evaluated. Recommendations are made for sites interested in measuring fluxes with the background concentration setup. The storage flux was computed using three-point profile concentrations routinely measured at the ICOS tower. The seasonal variations and variations with wind directions of the CO₂, CH₄, and CO fluxes were then briefly discussed.

2 Requirements and constraints of slow- and fast-response analysers

The benefits of expanding geographical coverage and increasing the range of measurable gases using slow-response analysers are clear. The benefits of using existing networks for background concentration measurements are clear as well. However, the constraints of concentration and flux measurements are not the same, and so not all towers may be suitable. For any atmospheric tower, a couple of adversities must be addressed first. Discarding atmospheric stations on mountains that have unsuitable conditions for flux measurements, we focus on tall towers over reasonably flat landscapes.

2.1 Time response and high-frequency attenuation

Measuring flux with a 3 s response time analyser is challenging. High-frequency measurements contain information that is overlooked when measuring at slower frequencies. Fast-response analysers, typically with 100 ms response time, are needed for flux measurements to capture the small and fast eddies (turbulent fluctuations) that carry most of the flux signal in the surface layer (Kaimal and Finnigan, 1994). Using slow-response analysers on short towers would mean losing most of the signal. Fortunately, the contribution of higher frequencies to the EC flux is inversely proportional to height (Horst, 1997).

For instance, using a simple estimator for the attenuation (Horst, 1997), we find that using a slow analyser with a time response of 3 s would attenuate the flux by 60 % in unstable conditions with wind speed of 4 m s⁻¹ for a 4 m tower, decreasing to 27 % on a 30 m tower and further decreasing to 8 % on a 100 m tower. Under stable conditions, the attenuation may increase to more than 90 %, 60 %, and 25 % respectively. Thus, at 100 m we expect only a limited amount of correction despite using a slow-response gas analyser.

High-frequency (HF) corrections based on predefined or experimental cospectra profiles are well established and routinely applied to correct for tube attenuations in ICOS and other flux networks (Horst, 1997; Massman and Lee, 2002; Ibrom et al., 2007; Fratini et al., 2012). We therefore expect sampling with slow-response analysers at tall towers may be suitable because the peak of the covariance cospectrum would be caught well and could be corrected with standardized approaches (Massman, 2000).

2.2 On the consequences of the height

Firstly, measurement height can impact footprint heterogeneity and non-stationarity. Height affects the source area of flux measurements, with taller towers having larger footprints. Data points are comparable only if their footprints are similar, which is less of an issue in homogeneous ecosystems. However, the ever-changing footprints due to wind conditions have a significant impact on the measured flux in heterogeneous environments. Furthermore, when these shifts happen during the averaging period, typically 30 min, they can lead to data being flagged as non-stationary and unusable with the standard eddy covariance method. This issue may be more pronounced for certain gases.

Secondly, height can have consequences for the terms of surface flux. Surface flux includes three terms: storage, advection, and turbulent transport (Finnigan et al., 2003). Storage flux arises from the accumulation or release of the compound below the measurement height. Advection is the transport by an organized flow, and it can occur when there is a non-zero mean wind velocity. Turbulent transport involves mixing caused by the eddies near the surface. While for short towers of only a few metres, the turbulent term predominates, taller towers require storage and vertical advection to be accounted for (Aubinet et al., 2005). As the measurement height increases, so does the volume below the measuring point, allowing bigger storage fluxes, and frictional effects decrease, allowing faster wind speeds. Larger-scale atmospheric circulation patterns also become more pronounced with height. While boundary layer dynamics may impact the surface fluxes, entrainment from the top of the atmospheric boundary layer may only have significant impact on lower frequencies (Asanuma et al., 2007). Here, we assume entrainment plays a negligible role in the turbulent fluxes.

Typically, positive storage flux may result from the decoupling of surface and atmospheric dynamics. Such a decoupling may arise especially under a stably stratified surface layer, occurring at night above canopies, especially under radiative cooling conditions (Kaimal and Finnigan, 1994). Negative storage fluxes arise during the early morning when the atmospheric boundary layer rises and the stably stratified layer breaks down (Aubinet et al., 2005). At tall towers, the storage can be high and remain large in the morning when the vegetation starts photosynthesizing but the turbulence is still low (Haszpra et al., 2005). At very tall towers (300 m), the storage dominates the flux dynamics (up to 95 % of the total flux, Winderlich et al., 2014). Under ideal surface homogeneous conditions, the storage term is expected to tend to zero when averaged over a day and hence only affect the surface flux dynamics but not the integrated fluxes. Over other conditions, this might be violated. The atmospheric boundary layer expansion and contraction may contribute to vertical wind dynamics and so to vertical advection as well. In some cases, including vertical advection was sufficient to reach a closed mass balance (Mammarella et al., 2007). Horizontal advection, on the other hand, is hard to measure with the current setups that focus on single-tower measurements. During decoupling it can have consequences by moving air parcels away from the measuring point, effectively renewing the air below measuring point with air from outside the area of interest.

For the turbulent fluxes the impact is less significant. In general, fluxes which would be measured as turbulent at a short tower may first appear as storage at a taller tower while being released as turbulent flux at a later stage. In addition to this, as we go farther from the sources and sinks, the concentration gradient, and thus the instant deviation, gets smaller. Consequently, this requires more accurate instruments.

3 Material and methods

3.1 Site description

The study uses data from a 100 m tall tower in the French Alternative Energies and Atomic Energy Commission (CEA) at a research campus in Saclay, 20 km southwest of Paris (Fig. 1). The tower is part of the ICOS atmospheric network (FR-Sac) and takes part in the ICOS Cities, Pilot Applications in Urban Landscapes (PAUL), project, focused on integrated city observatories for greenhouse gases. Climatically, the area is under oceanic influence, with mild temperatures (11.5 °C annual mean) and moderate precipitation (677–700 mm annual). The surrounded landscape is dominated by artificial (buildings, roads), agriculture (mainly cereal), and forest sources. The region serves as a pathway for urban-to-suburban daily mobility, with more than 60 000 vehicles every day in 2022 according to SIREDO on the national (N118) and regional (D306, D36, D128) roads.

Figure 1 Site diagram showing tower and acquisition house and tree for scale. In the right panel, the site map and the localization in the region are shown. Land use type is indicated in the colour legend. For reference, (1) heating plant, (2) manure/composting plant, and (3) lake.

Since 2011 the site has been equipped with high-precision cavity ring-down spectroscopy (CRDS) gas analysers (CO₂, CO, CH₄, G2401; Picarro Inc., Santa Clara, CA, USA) with a varying time response of a few seconds. The analyser is placed in a ground level hut, connected to three sampling lines, 12.7 mm of diameter, collecting air at the three different heights of the tower (15, 60, 100 m above ground level, a.g.l.) alternatively every 10 min. Since 2017, a second multi-gas analyser has been measuring continuously through a parallel sampling line connected to the top of the tower (100 m a.g.l.).

The flow rate through the sampling lines is set around 12 L min⁻¹ but with no control. At the bottom of all lines connected to the CRDS analyser, the air is dried with Nafion (PermaPure, model MD-070-144S-4). The CRDS gas analysers followed the ICOS calibration procedure, aiming for a precision higher than 50, 1, and 2 ppb for CO₂, CO, and CH₄ (ICOS RI, 2020). From June to October 2023, we set up a full eddy covariance system at 100 m, consisting of a closed-path infrared (IR) gas analyser (LI-7200; Li-Cor Inc., Lincoln, NE, USA), a 0.7 m heated tube with a flow rate set to 15 L min⁻¹, and a three-dimensional sonic anemometer (Gill WindMaster; Gill Instruments Ltd, Lymington, Hampshire, UK). The anemometer kept running after October 2023. Given the thin shape of the tower and the setup position on top, no major disturbances are expected from wind measurements. The tower is also equipped with pressure (Vaisala PTB200), humidity, and temperature sensors (Vaisala HMP155) at 1.5, 60, and 100 m.

Half-hourly average dry CO₂ mixing ratio showed a high degree of comparability between instruments (R² 0.97) and no bias (slope = 1) (Fig. 2). Nonetheless we found an offset of 7.25 ppm and an average drift of −11 ppm yr⁻¹, which has no impact on eddy covariance flux.

Figure 2 Comparison between dry CO₂ mixing ratio measured by the IR (LI-7200) and the CRDS (Picarro G2401) analysers. (a) Scatter plot. (b) Mixing ratio difference (IR-CRDS) as a function of time. Dots are observations, the red line is the linear fit, and the grey line is the 1:1 line. The correlation coefficient (R²), the mean error (ME; ppm), the mean absolute error (MAE; ppm), the linear fit, and the drift are shown.

Data were not available for most of November due to maintenance in the instruments. The IR (LI-7200) was set up unconventionally, with the analyser placed horizontally and the sampling tube vertically with a U-shaped head and a rain cap turned downwards. The choice was made based on the safety for maintenance on top of the tower. Unfortunately, the IR malfunctioned during the measurement campaign and prevented all the analysers from running at the same time for the whole period.

3.2 Data processing

The mass balance equation to compute surface flux includes three terms: storage, advection, and turbulent transport (Foken et al., 2012). With the setup of this work, horizontal advection was not addressable. The vertical component of the wind (w) and the mixing ratio of a scalar s (χ_s) can then be used to calculate flux at the surface F_eco based on the turbulent covariance ($\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\chi }}_s^{\prime }}$) measured at a certain height (h_m), the storage term (${\int }_{\mathrm{0}}^{h_m}\frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\chi }}_s}}{\partial t}\mathrm{d}z$, where t is time and z is the height) and vertical advection. Here overbars stand for time averaging. The surface flux F_eco ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) can then be expressed as a function of the molar volume of dry air $V_{\mathrm{a}}=\frac{R{T}_{\mathrm{a}}}{P-e}$, where P is atmospheric pressure (Pa), e is vapour pressure (Pa), R is the ideal gas law constant (8.31 J kg⁻¹ K⁻¹), and T_a is air temperature (K):

$$\begin{array}{}\text{(1)}& F_{\mathrm{eco}}=\left(\underset{\text{storage}}{\underbrace{\underset{\mathrm{0}}{\overset{h_m}{\int }}{V}_{\mathrm{a}}^{-\mathrm{1}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\chi }}_s}}{\partial t}}\mathrm{d}z}}+\underset{\text{vertical advection}}{\underbrace{\underset{\mathrm{0}}{\overset{h_m}{\int }}{V}_{\mathrm{a}}^{-\mathrm{1}}\stackrel{\mathrm{‾}}{w}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\chi }}_s}}{\partial z}}\mathrm{d}z}}+\underset{\text{turbulent}}{\underbrace{V_{\mathrm{a}}^{-\mathrm{1}}\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\chi }}_s^{\prime }}}}\right).\end{array}$$

3.2.1 Turbulent flux calculation

The turbulent flux was calculated based on covariance, hence the name of the eddy covariance (EC) method. Pre-processing is required and was done using EddyPro 7.0.9, applying de-spiking (Mauder et al., 2013), covariance maximization for time lag, and planar fit (Wilczak et al., 2001). Note that planar fit is required to be able to estimate vertical advection. Time lag relates to the delay from sampling and measurement, and the maximization can lose reliability under noisy measurements (Langford et al., 2015). Typically, a default value and bounds are set individually for each gas and gas analyser. If an optimal value falls within the bounds, it is retained; otherwise, the default is chosen. For the Li-Cor 7200 analyser the lag time was set to 0.09 ± 0.35 s based on tube dimensions and flow rate. For the Picarro analysers, which had a 100 m line, the lag time was set to 60 ± 2 s based on comparison with the Li-Cor 7200 CO₂ concentration. This lag time is compatible with a flow rate ∼ 12.6 L min⁻¹. The ±2 s tolerance was included to account for the uncertainty over the precise travel time and possible seasonal changes linked to air viscosity dependency on temperature and filter dirtiness.

Turbulent flux was calculated using EddyPro 7.0.9. For EC flux calculation, the slow-response analyser (Picarro) was resampled to 10 Hz to synchronize with the sonic anemometer sampling rate. This was achieved by repeating each measured value until it changed.

3.2.2 Storage and vertical advection flux computation

The storage and vertical advection were calculated using the three levels of measurements done by the same instrument, alternating between the three heights by periods of 10 min. The scalars CO₂, CO, and CH₄ were measured at 15, 60, and 100 m with the CRDS analysers. The dry air volume ratio was computed at each height based on measured air relative humidity and temperature. The 10 min average measurements were linearly interpolated, and 30 min averages were calculated.

The storage flux was computed as in Aubinet et al. (2005) as the derivative over time of scalar s contained in the column below the measurement height (100 m):

$$\begin{array}{}\text{(2)}& {\text{ST}}_s=\underset{\mathrm{0}}{\overset{h_m}{\int }}{V}_{\mathrm{a}}^{-\mathrm{1}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\chi }}_s}}{\partial t}}\mathrm{d}z\sim {\displaystyle \frac{\mathrm{\Delta }{\sum }_{\mathrm{1}}^{\mathrm{3}}{V}_{{\mathrm{a}}_i}^{-\mathrm{1}}\stackrel{\mathrm{‾}}{{\mathit{\chi }}_{{\mathrm{s}}_i}}\mathrm{\Delta }{z}_i}{\mathrm{\Delta }t}},\end{array}$$

where Δt is 30 min, index i stands for the three layers (0–15, 15–60, and 60–100 m), and Δz_i is the layer depth.

The vertical advection flux was computed using the derivative over the height of scalar s and integrated over the measurement height (100 m). Calculations were done in Python using NumPy (Harris et al., 2020) and SciPy (Virtanen et al., 2020).

3.2.3 Quality flags and stability classes

Quality flags were assigned using the standard 0–1–2 flag system from FLUXNET (Mauder and Foken, 2011), involving tests for stationarity and fully developed turbulence (Foken and Wichura, 1996). The stationarity test (STA) measures the absolute relative deviation between 5 and 30 min covariances, while the integral turbulence test (ITC) assesses the deviation between measured and modelled integral turbulent characteristics. Data are considered to be of high (< 30 %), medium (30 %–100 %), or low (> 100 %) quality, based on deviation percentages for each test (worst applicable result prevails). A detailed description of the quality flags can be found in Foken and Wichura (1996).

Stability classes were defined using the stability parameter $\mathit{\zeta }=(z-d)/L$, where z is the measurement height, d the zero-plane displacement height, and L the Obukhov length. We classified stability as unstable (ζ < −0.2), near-neutral (−0.2 $>\mathit{\zeta }>$ 0.2), and stable (ζ > 0.2).

3.3 High-frequency corrections on noisy measurements

Instruments have measurement limitations, which decreases their ability to produce a true value. Closed-path gas analysers require a gas sample to pass through a tube system including filters. Longer tube lengths typically result in increased time lag and reduced high-frequency signal. The signal degradation can be represented by a transfer function, TF, which attenuates the high frequency (Ibrom et al., 2007) of the true cospectrum of w and a compound s:

$$\begin{array}{}\text{(3)}& f{\text{Sp}}_{s,\mathrm{measured}}\left(f\right)=f{\text{Sp}}_{s,\mathrm{true}}\left(f\right)\times \text{TF},\end{array}$$

where f is the frequency (Hz), and Sp_s is the spectrum between w and a scalar s. Note that we can consider the transfer function equal for the spectrum and cospectrum, as we neglect the w transfer function and spatial sensor separation for the case of this tall tower (Massman, 2000). We assume the true covariance can be estimated by multiplying the measured covariance by a correction factor, CF:

$$\begin{array}{}\text{(4)}& {\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\chi }}_s^{\prime }}}_{\mathrm{true}}=\text{CF}\times {\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\chi }}^{\prime }s}}_{\mathrm{measured}}.\end{array}$$

Acknowledging that the covariance is the integral over all frequencies of the cospectra, the correction factor (CF) can be calculated from the transfer function (TF) and a true cospectrum, which is usually taken to be wT_s (where T_s is the ultrasonic temperature). Indeed, we assume similarity of scalars in the atmospheric boundary layer and use the (co)spectrum of T_s as a proxy of the unattenuated cospectrum, as T_s is collocated to w (Ibrom et al., 2007). For CF, this yields

$$\begin{array}{}\text{(5)}& \text{CF}={\displaystyle \frac{\int C{o}_{\mathrm{w}{T}_{\mathrm{s}}}\left(f\right)\mathrm{d}f}{\int C{o}_{\mathrm{w}{T}_{\mathrm{s}}}\left(f\right)\times \text{TF}\left(f\mathrm{|}{f}_{\mathrm{c}}\right)\mathrm{d}f}},\end{array}$$

where TF can be calculated in different forms and can account for both low- and high-frequency attenuation. Experimental methods are recommended for high-frequency spectral correction (Ibrom et al., 2007; Fratini et al., 2012). We can approximate an empirical TF, explained further down, using a first-order system, as the product of a transfer function H accounting for a first-order filter's time constant, τ_c, representing the system response time (s), and a transfer function H_p, accounting for a generic phase shift φ as (Massman, 2000)

$$\begin{array}{}\text{(6)}& {\displaystyle }\text{TF}=H\times H_{\mathrm{p}}\text{(7)}& {\displaystyle }H={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left(\mathrm{2}\mathit{\pi }f{\mathit{\tau }}_{\mathrm{c}}\right)}^{\mathrm{2}}}}\text{(8)}& {\displaystyle }{H}_{\mathrm{p}}=\cos \mathit{\phi }-\mathrm{2}\mathit{\pi }f{\mathit{\tau }}_{\mathrm{c}}\sin \mathit{\phi }.\end{array}$$

Note that the cut-off frequency, f_c, equals ${\left(\mathrm{2}\mathit{\pi }{\mathit{\tau }}_{\mathrm{c}}\right)}^{-\mathrm{1}}$. Ideally H would be the measured-to-true spectra ratio for the scalar of interest. However, only the measured spectrum is known, and so Eq. (7) is fitted using the sonic temperature T_s as a proxy of the unattenuated spectrum (Ibrom et al., 2007; Fratini et al., 2012; Peltola et al., 2021):

$$\begin{array}{}\text{(9)}& H={\displaystyle \frac{{\text{Sp}}_s\left(f\right)}{{\text{Sp}}_{s,\mathrm{true}}\left(f\right)}}\approx F_{\mathrm{n}}{\displaystyle \frac{{\text{Sp}}_s\left(f\right)/{\mathit{\sigma }}_s^{\mathrm{2}}}{{\text{Sp}}_{T_{\mathrm{s}}}\left(f\right)/{\mathit{\sigma }}_{T_{\mathrm{s}}}^{\mathrm{2}}}},\end{array}$$

where F_n is a normalization factor to account for any inaccuracies in the variance.

Sometimes TF=H is used, and H_p is not considered (Ibrom et al., 2007). However, not accounting for the phase shift (e.g. using cross-covariance maximization for lag correction and solely H for cospectra correction) can bias CF (Peltola et al., 2021). Fortunately, $H_{\mathrm{p}}\approx \mathrm{1}/\sqrt{H}$, which leads to $\text{TF}=H{H}_{\mathrm{p}}\sim \sqrt{H}$ (Peltola et al., 2021). In this work we use Fratini et al. (2012), where $\text{TF}=\sqrt{H}$.

The spectra and cospectra were calculated using EddyPro 7.0.9, following Fratini et al. (2012) described here in Eqs. (5), (7), and (9). The transfer function H, accounting for the first-order filter's time constant τ_c, was estimated for each analyser and each compound through a least-squares minimization approach of the spectra (Fig. 3).

Figure 3 Transfer functions H (dotted lines) for each compound and analyser. The transfer function was fitted to the ratio of each compound's spectra to the sonic temperature spectra. All spectra are ensemble-averaged taken from EddyPro outputs, filtered for significant fluxes. Dots show the mean spectra per frequency band. The grey-shaded area shows the frequency range (2–0.0018 Hz) over which transfer functions were fitted. See Fig. S1 in the Supplement for more detailed power spectra.

From the H, TF was computed as $\sqrt{H}$, and CF was calculated with Eq. (5). For both TF optimization and CF calculation, only frequencies between 2 and 0.0018 Hz were used (see Fig. 3).

3.4 Footprint and spatial tools

For an analysis of the fluxes' footprint as a function of wind direction, we used a simple parameterization of a backward Lagrangian stochastic particle dispersion model (LPDM-B) for the footprint (Kljun et al., 2015) and computed vegetation indexes based on Sentinel 2 (ESA/Copernicus Data) and a French land use map (Thierion et al., 2022). Note that a single roughness parameter was used for the footprint calculations, while in reality the covers and thus the roughness vary according to land use and wind direction. This source of uncertainty is reduced by taking measurements at 100 m a.g.l., farther away from the surface than typical flux towers.

3.5 Performance measurements

Comparisons between instruments were carried out using mean bias and absolute error, defined as

$$\begin{array}{}\text{(10)}& {\displaystyle }\text{mean error}\left(\text{bias}\right)={\displaystyle \frac{\mathrm{1}}{N}}\sum _{n=\mathrm{1}}^N\left(X_{a,n}-X_{b,n}\right),\text{(11)}& {\displaystyle }\text{mean absolute error}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{n=\mathrm{1}}^N\left|X_{a,n}-X_{b,n}\right|,\end{array}$$

where N equals the amount of data, and X is the variable measured with instrument a and b at a time n.

In figures, 95 % confidence interval bands were calculated using the Seaborn module in Python. It uses a random sampling method with replacement strategy, bootstrapping, to construct a confidence interval (Dragicevic, 2016).

For linear fits, if not declared otherwise, the squared loss, also named ordinary least-squares method, is used. The method consists of minimizing the sum of the squares of the difference between the observed and predicted values. When robust or Huber loss is mentioned, we use a linear fit, which minimizes the squared loss for the samples where the absolute difference between the observed, y, and predicted, f(x), values is smaller than δ and the absolute loss, sum of the absolute difference, otherwise. This feature makes it less sensitive to outliers than the squared error.

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \text{Huber loss}\\ & =\left\{\begin{array}{rl}\sum \frac{\mathrm{1}}{\mathrm{2}}{\left(y_i-f\left(x_i\right)\right)}^{\mathrm{2}},& \left|y_i-f\left(x_i\right)\right|\le \mathit{\delta }\\ \sum \mathit{\delta }\left(\left|y_i-f\left(x_i\right)\right|-\frac{\mathrm{1}}{\mathrm{2}}\mathit{\delta }\right),& \text{otherwise}\end{array}\right.\end{array}\end{array}$$

By default, we arbitrarily chose δ = 5; note that very low δ values may increase the number of values considered outliers.

4 Results

4.1 Mixing ratios of CO₂, CO, and CH₄

In Fig. 4, the diel pattern shows a peak in ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}$ during the morning (07:00 CET (UTC+1)) in July and moving towards 09:00 in October and a clear valley around 15:00. The pattern disappears when moving towards winter months. CO and CH₄ both show a similar peak in the morning in autumn, although less marked. Only CO shows an afternoon peak in September, which is also the month with the clearer morning peaks for CH₄ and CO. Seasonally, CO₂, CO, and CH₄ mixing ratios are the highest in winter. This difference may be explained by a larger biogenic CO₂ sink during the daytime in summer and a higher anthropogenic CO₂ emission in winter (heating on). The difference may also be explained by larger (smaller) boundary layer thickness during the summer (winter) which can effectively dilute (concentrate) the molecules emitted at the ground.

Figure 4 Monthly dry mixing ratio diel pattern for all measured gases (CO₂, CH₄, and CO) for IR (LI-7200) and CRDS (Picarro G2401). Solid line indicates median, and the shaded region shows 95 % confidence interval.

A look into how the mixing ratios vary with wind direction reveals some spatial patterns (Fig. 5). During warmer months (July to October), the west-wind CO₂ mixing ratio was smaller than the median value, while for CH₄ and CO we can notice a value that is higher than the median for the northeast sector, especially clear for CO. In colder months (December and January), all mixing ratios were higher (also seen in Fig. 4), with easterly winds (0–180°) showing larger mixing ratios than in other directions. A peak in mixing ratios is observed for all three gases for winds coming from around 20° N, the direction from the heating plant. Interestingly, a smaller peak can be seen in the northwest direction from the lake (100 m away), bare soil fields (around 500 m away), and a regional road roundabout (around 1 km away).

Figure 5 Median dry mixing ratios by wind direction; the band shows the interquartile (25th–75th percentiles). Warmer months (JASO, July to October) in grey and colder months (DJ, December and January) in red. Extreme values in the left and right 0.1 % tails were removed. See Fig. S2 for monthly observations. CO₂ is in parts per million (ppm), while CH₄ and CO are in parts per billion (ppb). All data are shown, including non-stationary and underdeveloped turbulence.

4.2 Footprint and stationarity

4.2.1 Footprint analysis

A characterization of the site's flux footprint (Fig. 6) shows a heterogeneous landscape composition, comprising 25 % urban, 23 % agriculture, 21 % forest, and 21 % grassland areas inside the 90 % level source area. In the western part of the site (42 % forest), there is a relatively dense woodland primarily featuring deciduous trees. To the south (41 % grassland), the landscape includes a nearby golf club in the vicinity of the CEA campus. In all directions there are croplands, predominantly cultivated cereal crops (winter wheat, barley, maize) and oilseeds (rapeseed), typical of the region. In the northeast (45 % urban), the landscape aligns with the location of the CEA campus, which includes a heating plant aligned to 20° N. The 43.9 km² 80 % footprint encompasses a national road (N118) and several regional roads (D306, D36, D128), with weekly traffic of 60 000 vehicles on average (in 2022 according to SIREDO). Water ponds have a small contribution for northwest to northeast sectors (2.3 %–2.6 %). Two ponds are situated in these sectors, one approximately 100 m from the tower (northwest) and a second larger farther away (around 2.4 km northeast, visible on the map in Fig. 6a).

Figure 6 Flux footprint by land use group. (a) Footprint for the whole period, where lines indicate 10 % to 90 % (border) level source area. (b) Contribution of each land use weighted by footprint density inside 90 % level source area. (c) Monthly footprint, where lines indicate 50 % and 80 % (border) level source area. Note that for visual purposes, the urban area is coloured white-grey in the map. Footprints estimated using model in Kljun et al. (2015).

In Fig. 6c, we can see the monthly changes in composition and shape of the area contributing to the fluxes measured at the tower (the flux footprint). Some months have larger footprints (e.g. August and September), while others have smaller ones (e.g. December and January), related to changes in the dispersion conditions. This difference is explained by the largest occurrence of stable conditions during the summer nights, which leads to larger footprints than during the winter, which mainly has neutral conditions (shown by the stability ratio ζ in Fig. 7b) driven by stronger winter winds, elevated friction velocity, and cloudy conditions (Fig. 7). Note as well that the landscape is not homogeneous (Fig. 6b), and so wind direction can also change the effective profile of sources and sinks contributing to each compound flux measured at the tower. December, for instance, was the month with the least contribution from the most urban northeast sector.

The mixing layer height, and similarly the atmospheric boundary layer height, shows a clear diurnal and seasonal cycles (Fig. 7). Warmer hours of the day and months show taller boundary layer heights, implying a larger volume of developed layer in which the compounds can be diluted. The mixed-layer heights were on average several times higher than the measuring height and so limiting any effect from entrainment. During these warmer periods, the conditions are often unstable ($z/L$ < 0.2), and friction velocity is high (> 0.4 m s⁻¹). This indicates well-mixed layer and bigger eddy sizes. Conversely, colder months (December and January) showed a relatively flat diel pattern, mostly due to a shorter photoperiod, leading to a much lower boundary layer height. We also noted on-site fog was frequently observed during these periods. Concurrently, friction velocity increased on average during winter. We note that horizontal winds, over heterogeneous terrain in particular in stable and neutral conditions, would favour horizontal advection. In the scope of the present work, however, it was not quantified due to a lack of measurements (horizontal gradients of concentration and fluxes).

Figure 7 Boundary layer conditions. (a) Heights of the atmospheric boundary layer (ABLH) and the mixing layer (MLH) measured at the SIRTA in Palaiseau, 4.8 km away from the tower. Data available online (Kotthaus et al., 2023). (b) Stability parameter ($\mathit{\zeta }=(z-d)/L)$ and friction velocity (u_∗) measured at the FR-Sac tower. Absolute values of ζ bigger than 2 were ignored.

4.2.2 Stationarity and well-developed turbulence

Most of the data collected were under well-developed turbulence, 73 % if only considering high-quality (flag 0) data based on the ITC test and 99 % including medium-quality (flag 1) data. Around half of the data (42 %), with an ITC flag of 0, were also considered stationary (Stationarity flag 0), increasing to 82 % if we include a flag of 1 on both tests. The stationarity test is required for standard EC; thus the use of the latter increases the data amount by 34 % in the case that only high-quality observations are used and 55 % in the case that medium-quality data are included (Fig. 8). This saving happens more often during the day due to a higher coincidence of both flags during night. The percentages given are for the Li-Cor (IR) fast analyser but are of the same order of magnitude as for the Picarro analyser (CRDS).

Figure 8 Quality control flags for turbulence (ITC) and CO₂ stationarity (SS). Flags follow a 0–1–2 system for high, medium, and low quality. Percentage of (a) turbulence flagged data by hour of the day. Stationarity flagged CO₂ data by hour of the day for (b) the Li-Cor instrument (IR) and (d) the Picarro instrument (CRDS). Stationarity flags per ITC group are also given for the IR (b) and CRDS instrument (c). Percentages are summed to 100 % in each group and over all data (in parentheses).

4.3 Comparison of CO₂ flux between slow- and fast-response analysers

4.3.1 High-frequency spectral correction

The CRDS analysers showed significantly more high-frequency attenuation of the flux than IR analysers (Table 2 and Fig. 9), as expected, due to the much longer sampling tube of the CRDS analyser (115 m) than the IR analyser (0.7 m), as well as the slower CRDS acquisition frequency (∼ 3 s) compared to the IR (0.1 s). The difference was greater in (very) stable conditions, when higher frequencies contribute more to the flux, than in (very) unstable conditions (Fig. 9). In (very) unstable conditions, the contribution of low frequencies to the flux increased, as shown by the fact that none of the ogive levelled to 1 towards 30 min integration time (ogive slope > 0). Surprisingly CO (measured by CRDS) showed an atypical curve with stronger contribution from high frequencies, which after analysis was attributable to amplitude resolution from a signal with weak variance and a less sensitive instrument. This shape explains the higher cut-off frequency of the spectral corrections in comparison with CO₂ and CH₄ (Table 1).

Figure 9 Normalized cospectra (a) and ogives (b) of w and CO₂, CH₄, and CO covariances, for gases measured by IR and CRDS and the reference sonic temperature, T_s. Median values from July to October 2023 grouped by stability classes: ζ < −0.2 (unstable), −0.2 > ζ > 0.2 (near-neutral), and ζ > 0.2 (stable). N indicates the amount of half hourly data in each class.

Table 1 The transfer function parameters for each instrument accounting for high-frequency attenuation. Here $\text{TF}={\left(\mathrm{1}+{\left(\mathrm{2}\mathit{\pi }f{\mathit{\tau }}_{\mathrm{c}}\right)}^{\mathrm{2}}\right)}^{-\frac{\mathrm{1}}{\mathrm{2}}}$, where τ_c is the first-order filter's time constant. The cut-off frequency, f_c, equals ${\left(\mathrm{2}\mathit{\pi }{\mathit{\tau }}_{\mathrm{c}}\right)}^{-\mathrm{1}}$. F_n is a normalization factor. The optimized values correspond to optimizations as shown in Fig. 3. See Eqs. (5)–(9) and text for details.

The high-frequency attenuation was around 2 %–3 % for the fast instrument (IR), while for the CRDS instruments sampling at 100 m, it ranged from 11 % to 19 % (Table 2). We can expect larger corrections in stable conditions, characterized by a larger contribution of high frequencies to the flux, as observed for IR; contrarily CRDS shows a decrease compared to near-neutral. It is worth noting that despite the 10 Hz acquisition frequency and 100 m height, the attenuation of the IR instrument was non-negligible. Additionally, the time response of the slow CRDS analysers, estimated based on the transfer function (3.62 s), matches the acquisition frequency (ranging between 3 and 4 s), but it also matches the expected attenuation for a long tube (Fig. 13).

Table 2 Percentage of high-frequency corrections of the CO₂, CH₄, and CO fluxes per stability class for each instrument.

The median random uncertainty for CO₂ was approximately 11 % using the fast setup (IR), comparable to the slow setup (CRDS). In contrast, uncertainties were higher for CH₄ and CO, reaching 23 % and 36 % respectively. The higher uncertainty is expected due to smaller fluxes and lower instrument sensitivity.

4.3.2 Comparing CO₂ flux measured by slow- and fast-response analysers

The CO₂ fluxes computed from the IR (LI-7200) and the CRDS (Picarro) analysers were well correlated with an underestimation of 18 % of the CRDS for uncorrected fluxes that was diminished to 2 % after high-frequency corrections (Fig. 10). High-frequency correction decreased the bias, ME, by 0.18 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ and the mean absolute error (MAE) by 0.11 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, with no effect on the correlation coefficient (R²). Moreover, there was a tendency of the CRDS corrected fluxes to slightly underestimate the CO₂ fluxes under stable conditions (Fig. S3).

Figure 10 Comparison of CO₂ flux computed with the IR (LI-7200) and the CRDS (Picarro G2401) analysers, (a) before and (b) after high-frequency loss corrections. Dots are observations, red line is a robust linear relation and grey line is the 1:1 line. The correlation coefficient (R²), the mean error (ME, $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), the mean absolute error (MAE, $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). Statistics are calculated ignoring outliers from robust linear regression.

4.4 Turbulent flux dynamics

4.4.1 Diel pattern over the months

We observed a well-defined summer pattern for the CO₂ flux, with emissions during the night and sequestration during the day (Fig. 11). From summer to winter, the sink shortens in time and decreases in magnitude up to the point that during winter, the site behaves on average as a source all along the day. We note that the concentration morning peak observed in Fig. 4 does not correspond to a peak in the flux. Following the seasonal pattern of CO₂, satellite data show less green leaves from September on (Fig. S5).

Figure 11 Monthly mean turbulent flux diel pattern of CO₂, CH₄, and CO for IR (LI-7200, only CO₂) and CRDS (Picarro G2401) gas analysers. Fluxes shown after spectral correction. Data points falling within the extreme 1 % tail of the distribution were removed. Shaded area refers to standard deviation.

4.4.2 Flux by wind direction

The CH₄ fluxes showed a marked daily pattern from July to September, with higher emissions in the morning than in the afternoon, while daily emissions were overall higher in January. Looking at the CO flux, we see a marked increase in November and January but not in December, despite similar temperatures and traffic. In January, winds were relatively well distributed, while in December the most urban NE sector was rarely in the footprint, which may explain the difference between the 3 months.

Daytime and night-time turbulent fluxes reveal differences depending on wind direction and season (Fig. 12). During warmer months (January to October), the site was a daytime net sink of CO₂ when the wind was coming from the forest in the west. Nocturnal wind coming from the SW (∼ 210°) was an important source. No sector was identified as a sink for CH₄ or CO, and emissions were higher for both gases from eastern winds with a peak at around 40°.

Figure 12 CO₂, CH₄, and CO turbulent fluxes by wind direction. Warmer months (JASO) in panels (a), (b), and (c) and colder months (DJ) in panels (d), (e), and (f). The values are the median with the interquartile range in the shaded area; 0.1 % extreme values were removed. Wind direction bins with fewer than 10 observations were added to the next bin clockwise. CO₂ fluxes are in $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, and fluxes of CH₄ and CO are in nmol m⁻² s⁻¹. Note that 10–50° for CO₂ and CH₄ in December and January was plotted separately for visual purposes. See Fig. S4 for monthly values. Only stationary periods and those with well-developed turbulence are considered.

During colder months, daytime and night-time fluxes become similar. During these months when the leaves have fallen, the CO₂ sink becomes a source. CO emissions in the northeast sector increase, possibly due to the alignment with the national road, N118. Additionally, the wind direction spanning 10–45°, the direction of the local heating plant, exhibited CO₂ and CH₄ emissions significantly higher than those observed in other directions, during both daytime and night-time periods. CO also increases at around 20° but with a much smaller magnitude.

5 Discussion

5.1 Challenges of measuring on a tall tower with slow-response analysers

5.1.1 High-frequency loss corrections on the atmospheric tower configuration

Our findings revealed that an ICOS atmospheric tower configuration utilizing a CRDS gas analyser with an acquisition frequency of approximately 0.3 Hz and a tube length of 100 m exhibited a high-frequency loss correction of approximately 20 %. This correction was around 3 times more than the conventional ecosystem flux measurement setup, which employed an IR gas analyser with a 10 Hz acquisition frequency and a tube length of 0.7 m, positioned at the top of the tower. The observed transfer function (TF) for the CRDS setup closely matched the theoretical attenuation expected, as depicted in Fig. 13. Indeed, the tube and sensor attenuation together lead to a first-order time constant around 3 s, as we observed for the CRDS setup (Table 1). This outcome suggests that even with a faster measurement system or a smaller tube attenuation, only a limited reduction of the attenuation can be expected. In order to substantially decrease the high-frequency attenuation of the flux, both an increase of the acquisition frequency and a decrease in tube attenuation (decrease in tube length or increase in flow rate) would be required.

Figure 13 Transfer functions computed for (1) the sensors' acquisition frequencies (TF_acq, Horst, 1997), (2) the tube attenuation (TF_tube, Leuning and Moncrieff, 1990; Foken et al., 2012), (3) the combination of (1) and (2), and (4) the CO₂ observations, for the two setups at the tower. The Saclay atmospheric setup (CRDS) consists of a 100 m sampling line with a 9.5 mm diameter, a sampling frequency of ∼ 0.3 Hz, and a flow rate of 12.7 L min⁻¹. The conventional ecosystem setup (IR) consists of a 0.7 m sampling line with a 5.33 mm diameter, a sampling frequency of 10 Hz, and a flow rate of 15 L min⁻¹. Note that for the IR setup, the curves (1) and (3) are superposed.

We note that since the tube attenuation is higher when flow inside the tube is laminar (Lenschow and Raupach, 1991), ensuring a Reynolds number larger than ∼ 2300 is key to minimizing attenuation. We can define $\mathit{Re}=\frac{\mathrm{2}Q}{\mathit{\pi }rv}$, where Re is the Reynolds number (–), r is the tube radius (m), Q is the volumetric flow rate (m³ s⁻¹), and v is the kinematic viscosity of air. We find that for the tube in place with 9.5 mm of internal diameter, pumping ∼ 14–17 L min⁻¹ is necessary to achieve a turbulent flow. Under these conditions, the cut-off frequency would increase to more than 0.6 Hz, but the pressure would also drop from −6 to −47 mbar.

We can also notice a high-frequency peak around 10^−0.85 Hz (7 s) in the observations, which might come from resampling the CRDS concentrations to 10 Hz. The resampling by repeating the same value 30 times creates sharp corners in the time series, which might create noise and systematic problems in a Fourier transform. In the future, other high-frequency gap filling methods can be tried. However, problems related to the resampling happen at frequencies higher than setup attenuation and have limited effect on the transfer functions and corrections. It is worth mentioning that we used a first-order filter fitted on in situ data as a transfer function, following Fratini et al. (2012) and shared by other studies (Ibrom et al., 2007; Peltola et al., 2021). In the atmospheric tower configuration, where the main attenuation arises from the tube length, the transfer function may take an exponential shape, as proposed by Leuning and Moncrieff (1990) and Foken et al. (2012), and the fitting may not be perfect, as depicted in Fig. 13. The effect was however evaluated to be negligible on the correction factor.

For the same transfer function, attenuation may change based on the cospectra dependence on measuring height, wind speed, and stability parameter ($z/L$). Specifically, increases in wind speed and stability parameter, or decreases in measuring height, are expected to shift a cospectrum towards higher frequencies, thereby enhancing attenuation for a given transfer function (Horst, 1997). Theoretical expectations of the attenuation factor from Horst et al. (1997) based on empirical cospectra agree very well with our measurements under unstable to near-neutral conditions but do not entirely align with our observations for neutral and stable conditions (Fig. 14). Indeed, surprisingly, we found that the attenuation remained stable or slightly decreased for $z/L$ values over 0.2 in the case of IR and CRDS respectively. This contrasts with the prediction by Horst et al. (1997), which suggested an increase by a factor of 5 under very stable conditions. This difference needs further investigation.

Figure 14 Flux attenuation due to high-frequency losses, theoretical as lines and measured as points. The theoretical losses are computed from Horst (1997, Eq. 11), using the measured first-order time constant τ_c for the IR (0.5 s) and the CRDS (3 s).

More surprising is the HF losses found for CO₂ measured by the IR (LI-7200). On this supposedly conventional ecosystem setup (0.7 m heated tube, 5.33 inner diameter, 15 L min⁻¹ flow rate, 7 µm filter) attenuation and the transfer function were expected to be much smaller. The time constant of 0.5 s is equal to the cutting frequency (0.32 Hz) reported by Ibrom et al. (2007) for CO₂ with 50 m long tube, 8 mm diameter, and a flow rate of 20 L min⁻¹. We do not have strong evidence to explain this very large attenuation, and we are bound to speculate that this may be due to the inlet filter. Indeed, we observed a very large HF loss for H₂O, with a τ_c ranging from 0.7 s for RH < 30 % to 5 s at 70 % RH and 50 s for RH larger than 80 % (data not shown). This is the sign that the inlet did accumulate water vapour, most probably on the filter holding hygroscopic aerosols. Since CO₂ dissolves in water, the microscopic water accumulated in the tube may have buffered the CO₂, leading to a large attenuation.

The CRDS setup, however, exhibited relatively small attenuation. In comparison, Wintjen et al. (2020) reported a damping factor of around 16 %–22 % for a 48 m tube with a 6 mm diameter, measuring reactive nitrogen at 10 Hz. Despite the longer tube length and slower analysers in our study compared to Wintjen et al. (2020), their slower flow rate (2.1 L min⁻¹) and the expected stronger air–wall interactions for reactive nitrogen compounds may have contributed to the higher damping factor in their study. Correcting for high-frequency losses resulted in agreement between the IR and CRDS methods within 3 %, maintaining the elevated R² of 0.94 (Fig. 10). This demonstrates that the high-frequency correction was able to correct for the losses.

We assumed all compounds (CO₂, CH₄, CO) measured by CRDS (Picarro G2401) suffered the same attenuation and used CO₂, the best-defined curve for all three analysers (Table 1). This assumption is grounded on the fact that measurements are done by the same instrument at the same acquisition rate sampled through the same line and is backed by the proximity between CH₄ and CO₂ spectra, while the unexpected CO spectra can be explained by the CO signal noise due to the lower signal-to-noise ratio. Indeed, the noise was already larger than the signal at periods larger than 5 min (Fig. 3). Similarly, the small step increase around 4 s for the CRDS analysers corresponds to the actual measurement interval.

Furthermore, while this paper primarily examines the high-frequency aspect of the signal, it is important to note that the 30 min integration period did not allow the low frequencies of the fluxes to be captured entirely, especially under unstable conditions, as can be seen by the ogive slope being non-zero at the lowest frequencies (Fig. 9). A recent study on urban tall towers also reported low-frequency contribution for kinematic heat and CO₂, indicating the importance of low-frequency corrections (Lan et al., 2024).

5.1.2 Storage and advection terms

The current paper focused on calculating turbulent flux using a setup dedicated to measure background concentration. To estimate the surface flux, storage and advection must be considered. The typical atmospheric setup also performs profile measurements, which can provide data to estimate both storage and vertical advection terms. The profile setup is not ideal since measurements at each height are not taken simultaneously, and the 20 min delay between measurements at each height can induce errors. A typical error can happen when for instance an urban plume rich in CO₂ and CO is passing through the tower and both the asynchronism and delay in measurements can introduce artefacts in the fluxes. However, the estimations are still helpful and may give an order of magnitude and general patterns when aggregated.

The CO₂ storage flux amplitude was of the same order of magnitude as the turbulent fluxes (Fig. 15). This agrees with observations in Haszpra et al. (2005) and means that the storage term is important when looking into the diel surface flux pattern. In particular we observe a strong emission during sunset not captured by the turbulent flux and a clear negative storage flux in the morning in July. A similar diel pattern can be seen in CH₄ storage fluxes but not for CO. The early release in the case of CO₂ may account for both the early onset of photosynthesis and an artefact from the expansion of the mixed layer (Fig. 7). The difference between this early negative storage flux between CO₂ and CH₄ may indicate the weight of these two causes, as CH₄ is independent of photosynthesis. Atmospheric dilution results in both negative storage and positive vertical advection; however it is not a surface emission nor removal and cancels out if all conservation equation terms are considered.

Figure 15 Diel mean for CO₂, CH₄, and CO turbulent, storage, and vertical advection fluxes in panels (a), (b), and (c) and their sum, i.e. estimated surface flux, in panels (d), (e), and (f). Data for July 2023. Values for the CRDS (Picarro G2401). The band shows the interquartile (25th–75th percentiles). See Fig. S9 for all months.

Vertical advection depends on the mean vertical component of the wind (w) and thus on the quality of the tilt correction applied. On average, w was positive during unstable and near-neutral conditions ($z/L$ < 0) and near zero for stable conditions (Fig. S8). For vertical advection, however, monthly averages were near zero for all measured gases (Figs. 15 and S9). Indeed, in the long term, vertical advection could be neglected for very tall towers based on the assumption that synoptic-scale processes should counterbalance vertical advection (Davis et al., 2003).

In this work horizontal advection was assumed negligible, but in reality, it may not be (Aubinet et al., 2005). Indeed, daily mean storage in the site gravitates around zero during July but shows a non-negligible variability for CO₂ in August and September and for all gases in January (Fig. S6). In ideal conditions, the storage flux should average 0 over 24 h since what is stored at night should be de-stocked during the day. A non-zero daily storage may indicate advection.

Additionally, the chimney of the heating plant, ∼ 600 m away from the tower, situated at a certain height, may bias the storage fluxes, which rely on gradient profiles along the tower. Indeed, if the plume from the heating plant emissions is measured intermittently due to changes in wind directions, we may attribute or ignore storage fluxes where lateral advection is happening. Properly identifying such a process would require tracking the chimney plume with a 3D dispersion model and half-hourly resolution which was out of the scope of the present work.

5.2 Wind directional interpretations

Monitoring fluxes over a heterogeneous landscape demands spatial awareness. Turbulent fluxes were shown to vary depending on wind direction and season (Fig. 12). Notably the heating plant direction could be easily detected by the higher levels of CO₂ and CH₄ flux. Furthermore, winds originating from the deciduous forest (west) went from a daytime CO₂ sink in warmer months to a source in winter when leaves had fallen.

Comparison with a mixed deciduous forest site at 50 km SE (FR-Fon ICOS site, e.g. Delpierre et al., 2016) shows a similar seasonality of the CO₂ fluxes, storage terms being slightly larger in FR-Sac than in FR-Fon, and a turbulent term notably smaller in FR-Sac (Fig. S7). Differences may arise mainly from the proportion of the forest on the respective site's footprints. The forest represents around 40 % of the FR-Sac Western footprint and is expected to be higher during the day when unstable conditions shorten the footprint range. The remaining source area includes urban and traffic, which may add a positive component to the flux, offsetting the net flux during the day (Fig. 6). Although the footprint model does not account for variations in surface roughness, the results are still expected to remain valid despite this limitation.

Interpretations based on wind direction and comparisons with other sites serve the purpose of enlightening the potential of recovering spatially explicit fluxes from tower measurements. Nonetheless, accurate spatialized surface fluxes at finer scales remain a scientific challenge. The more complex discussions on relating these tall tower turbulence measurements to surface emissions and removals, as well as the intricate footprint analyses, should be further explored. Recent developments in the literature have shown promising results in addressing spatial variability. Indeed, research has been done on mapping fluxes from single-tower measurements by means of regressions (Crawford and Christen, 2015), Bayesian methods (Levy et al., 2020), or machine learning algorithms (Metzger et al., 2013; Metzger, 2018; Xu et al., 2018, 2020). For those efforts, capturing low-frequency fluxes as done in the “virtual control volume” in Metzger et al. (2018) could benefit the spatialization as they may be important, at least in urban tall towers (Lan et al., 2024). At the same time, maintaining a process-based approach could be advantageous for upscaling and understanding the ecophysiological processes in play. For flux attribution, a defensible partitioning method to decompose net fluxes under heterogeneous landscape could help map gases with considerable sources and sinks, such as CO₂.

5.3 Recommendations for atmospheric sites concerned by such a method

As of now, the ICOS network comprises 47 atmospheric sites, with 39 labelled class 1 or 2 and the remainder as on track for labelling. The ecosystem centre is more extensive, encompassing 103 sites, including 77 labelled class 1 or 2, and the remaining associated sites. This count may further rise when other regional networks are taken into consideration.

Not all atmospheric sites are adapted for flux measurements. EC towers prioritize flat surfaces, slim towers, and homogeneity, whilst atmospheric towers may prioritize locations based on grid redundancy to improve atmospheric inversions. Slim towers with limited topography around them are recommended for reliable measurements. Atmospheric measurements conducted in close proximity to large structures (e.g. just above domes) or in mountainous regions can introduce disturbances in the turbulence signal. This can lead to unreliable tilt angle correction and surface flux assessments. In some cases, flagging wind sectors that are not appropriate for EC measurements can be a straightforward solution. With that said, eddy covariance has successfully been used in mountainous landscapes using appropriate tilt corrections (Matthews et al., 2017).

For atmospheric tower research candidates interested in measuring flux, we recommend the following:

• Select at least one height for calculating fluxes through eddy covariance. For this decision, a footprint estimation (Kljun et al., 2015) may be relevant.

• Include a high-frequency 3D anemometer on the chosen height(s).

• Evaluate the first-order filter time for tube and sensors, as shown here, to verify that high-frequency attenuation is below an acceptable threshold, ∼ 20 %.

• Perform continuous mixing ratio measurements on the chosen height(s), either limiting profile measurements to specific hours or with a separate set of instruments. This ensures that low-frequency signal for eddy covariance and profile measurements can prioritize transition periods when fluxes may exhibit non-stationary behaviour or low turbulence.

• When performing profile measurements, alternate between levels every 5 min can improve the storage estimation.

• Evaluate the flow regime in the sampling tube, if possible increasing flow rate to guarantee turbulence (Re > 2300).

• Additional meteorological data (e.g. precipitation, short- and longwave incoming and outgoing radiation) and metadata (e.g. forest type, crops, transport counting) pertinent for flux interpretation should also be collected.

Note that adding a high-frequency 3D anemometer is enough for reanalysis. Ensuring continuous measurements is important for having a better-defined turbulent cospectrum and keeping low-frequency information. Certainly, profile measurements are often contemplated, and doubling the instruments may not be feasible. Therefore, we recommended restricting profile measurements to specific hours when the development of the boundary layer may overshadow the relevance of measuring flux close to the tower. These moments, typically during sunrise and sunset, provide valuable insights from both atmospheric and ecosystem (storage) perspectives. Furthermore, during these moments, standard covariance would typically flag and disregard measurements, as depicted in Fig. 8.

6 Conclusions

This study serves as a proof of concept for leveraging existing atmospheric towers to measure fluxes by simply adding a 3D anemometer. While eddy covariance on tall towers introduces challenges related to heterogeneity and storage effects, it mitigates concerns such as high-frequency attenuation. Comparing slow-response analysers with fast-response ones revealed very similar net fluxes across all stability conditions (R² > 0.94), indicating the viability of using slower instruments in this case.

It is important to note that our results focus solely on passive gases, as water was not considered due to air drying before measurement in the CRDS. For water, we could expect greater attenuation linked to tube length, as air–wall interactions of absorption and desorption are much stronger in water vapour (Massman and Ibrom, 2008). Similarly, we would not recommend measuring reactive gases, as their residence time might be too long (∼ 60 s) for accurate eddy covariance measurements to be made.

While many of the variables affecting attenuation are not under the researcher's control, limited choices remain for measurement height, tube dimensions, flow rate, and acquisition frequency. Thus we recommend continuous gas measurements to be systematically done with a high-frequency 3D anemometer and a flow rate sufficiently large to ensure turbulent flows in the sampling tube.

To calculate the surface flux, we estimated the storage and vertical advection terms. Although the storage term calculated using three heights provided useful estimations, caution is warranted due to potential biases from not measuring height at the same time but also due to the limited number of heights sampled. Our results underscore the significance of the storage term, which was as large as the turbulent flux at the measurement height.

Analysing fluxes by wind direction revealed distinct patterns, particularly between the forest (W) and campus site (NE). Notably, emissions from a heating plant significantly influenced CO₂ and CH₄ fluxes in colder months, highlighting the importance of considering local sources. While our findings align with anticipated patterns across various land uses, accurately attributing fluxes to land uses would necessitate additional modelling efforts, which were beyond the scope of this study.

Overall, this study demonstrates the potential of expanding flux measurements through a relatively inexpensive instrumentation addition, offering valuable insights for both ecosystem and atmospheric research. It further shows that the eddy covariance method has sufficiently matured so that we can use less-than-ideal instrumentation.

Appendix A

The theoretical approach for high-frequency loss corrections requires defining a transfer function, TF, for each of the relevant origins of frequency losses, i, and multiplying them to find the total transfer function, TF_total (Moore, 1986):

$$\begin{array}{}\text{(A1)}& {\text{TF}}_{\mathrm{total}}\left(f\right)=\prod {\text{TF}}_i\left(f\right).\end{array}$$

Note that TFs range between 0 and 1, and so the TF_ws is driven by the most restrictive function for each frequency. Considering only the attenuation from the air transport in the tube (TF_tube, Leuning and Moncrieff, 1990; Foken et al., 2012) and acquisition rate (TF_acq, Horst, 1997),

$$\begin{array}{}\text{(A2)}& {\text{TF}}_{\mathrm{tube}}=\left\{\begin{array}{rl}\exp \left\{-\mathrm{160}{\mathit{Re}}^{-\frac{\mathrm{1}}{\mathrm{8}}}\frac{{\mathit{\pi }}^{\mathrm{2}}{r}^{\mathrm{5}}{f}^{\mathrm{2}}L}{Q}\right\},& \mathit{Re}<\mathrm{2300}\\ \exp \left\{-\frac{{\mathit{\pi }}^{\mathrm{3}}{r}^{\mathrm{4}}{f}^{\mathrm{2}}L}{\mathrm{6}{D}_{\mathrm{s}}Q}\right\},& \mathit{Re}\ge \mathrm{2300}\end{array}\right.\end{array}$$

$$\begin{array}{}\text{(A3)}& {\text{TF}}_{\mathrm{acq}}={\left[\mathrm{1}+{\left(\mathrm{2}\mathit{\pi }f{\mathit{\tau }}_{\mathrm{w}}\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}/\mathrm{2}}\times {\left[\mathrm{1}+{\left(\mathrm{2}\mathit{\pi }f{\mathit{\tau }}_{\mathrm{s}}\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}/\mathrm{2}},\end{array}$$

where Re is the Reynolds number (–); r is the tube radius (m); f is the frequency (Hz); L is the tube length (m); Q is the volumetric flow rate (m³ s⁻¹); and τ is the first-order filter's time constant (s), where $\mathit{\tau }={\left(\mathrm{2}\mathit{\pi }{f}_{\mathrm{acq}}\right)}^{-\mathrm{1}}$ and f_acq is the acquisition frequency (Hz) for vertical wind speed or scalar. The Reynolds number is defined as $\mathit{Re}=\frac{\mathrm{2}Q}{\mathit{\pi }rv}$, where v is the kinematic viscosity.

Attenuation also depends on the cospectra. A theoretical approach is proposed in Horst (1997, Eq. 11), where

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{{\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}}_{\mathrm{measured}}}{{\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}}_{\mathrm{true}}}}=\mathrm{1}+{\left(\mathrm{2}\mathit{\pi }{n}_{\max }{\mathit{\tau }}_{\mathrm{c}}\stackrel{\mathrm{‾}}{u}/z\right)}^{-\mathit{\alpha }},\end{array}$$

where $\stackrel{\mathrm{‾}}{u}$ is the mean wind speed at height z; n_max is 0.085 in case $z/L$ < 0, else $\mathrm{2}-\mathrm{1.915}/(\mathrm{1}+\mathrm{0.5}z/L)$; and α is $\mathrm{7}/\mathrm{8}$ for $z/L$ < 0, else 1.