# Improvement of CO2 flux quality through wavelet-based Eddy Covariance: a new method for partitioning respiration and photosynthesis

Pedro Henrique Herig Coimbra, Benjamin Loubet, Olivier Laurent, Matthias Mauder, Bernard Heinesch, Jonathan Bitton, Nicolas Delpierre, Jérémie Depuydt, Pauline Buysse

► To cite this version:

Pedro Henrique Herig Coimbra, Benjamin Loubet, Olivier Laurent, Matthias Mauder, Bernard Heinesch, et al. Improvement of CO2 flux quality through wavelet-based Eddy Covariance: a new method for partitioning respiration and photosynthesis. 2023. hal-04272798

HAL Id: hal-04272798

https://hal.science/hal-04272798v1

Preprint submitted on 6 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Improvement of CO₂ flux quality through wavelet-based Eddy Covariance: a new method for partitioning respiration and photosynthesis

Pedro Henrique H. Coimbra<fnref n="1" />*, Benjamin Loubet<fnref n="1" />, Olivier Laurent<fnref n="2" />, Matthias Mauder<fnref n="3" />,<fnref n="4" />, Bernard Heinesch<fnref n="5" />, Jonathan Bitton<fnref n="5" />, Nicolas Delpierre<fnref n="6" />,<fnref n="7" />, Jérémie Depuydt<fnref n="1" />, Pauline Buysse<fnref n="1" />

## Abstract.

Eddy Covariance (EC) is praised for producing direct, continuous, and reliable flux monitoring for greenhouse gases. For CO₂, the method has been commonly used to derive gross primary productivity (GPP) and ecosystem respiration (R<sub>eco</sub>) from net ecosystem exchange (NEE). However, standard EC is impacted by non-stationarity, reducing data quality and consequently impacting standard partitioning methods that are constructed on simplistic assumptions.

This work proposes a new wavelet-based processing framework for EC tested over two French ICOS ecosystem sites, a mixed forest (FR-Fon) and cropland (FR-Gri), over several years. A new direct partitioning method was also developed, applying conditional sampling in wavelet decomposed signals. This new empirical method splits positive and negative parts of the wavelet decomposed product of the wind vertical component and CO₂ dry molar fraction, conditioned by water vapour flux, to compute GPP and R<sub>eco</sub>.

Results show 17 to 29 % fewer gaps in wavelet-based than with standard EC, varying between sites, day and night. A good correlation between methods was observed during turbulent and stationary periods (R<fnref n="2" />=0.99). However, wavelet-based NEE was 9% lower than in standard EC, likely related to low-frequency attenuation led by the detrending nature of wavelet transform.

The new wavelet-based direct partitioning provided daily GPP and R<sub>eco</sub> very similar to night- and day-time model-based partitioning methods, with the difference between our method and these standard methods smaller than between them. Our method did not produce positive GPP, a common error in the night-time method. It also showed R<sub>eco</sub> seasonal patterns coherent with management practices at the crop site (growing season, harvest, manure application), which was not always the case for the standard methods. The R<sub>eco</sub> diel cycle was noticeably different, whereas the standard methods are

<footnote n="1">¹ ECOSYS, INRAE, AgroParisTech, Université Paris-Saclay, Palaiseau, France</footnote>

<footnote n="2">² Laboratoire des Sciences du Climat et de l'Environnement, CEA, CNRS, Université Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote n="3">³ Institute of Meteorology and Climate Research - Atmospheric Environmental Research (IMK-IFU), Karlsruhe Institute of Technology, Garmisch-Partenkirchen, Germany</footnote>

<footnote n="4">⁴ Institute of Hydrology and Meteorology, Technische Universität Dresden, Dresden, Germany</footnote>

<footnote n="5">⁵ Faculté des Sciences Agronomiques de Gembloux, Unité de Physique, Gembloux, Belgium</footnote>

<footnote n="6">⁶ Ecologie Systématique Evolution, CNRS, AgroParisTech, Université Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote n="7">⁷ Institut Universitaire de France, France</footnote>

<footnote>* corresponding author: pedro-henrique.herig-coimbra@inrae.fr</footnote>
temperature-driven; our method had a daily pattern correlated to solar radiation and a night-time pattern correlated to soil temperature.

Eddy-covariance; wavelet; CO2 flux; partitioning; photosynthesis; respiration;

# 1. Introduction

Global surface temperature is 1.1 °C warmer compared with the pre-industrial era, even larger on land (1.6 °C) (IPCC, 2021). Climate change leads to widespread adverse impacts and related losses and damages to nature and people (IPCC, 2021b). Projections show more than 2°C warming in 2100 in part due to mismatches between implemented policies and long-term goals (IPCC, 2021a). The warming results from the increase in greenhouse gases (GHG) concentration in the atmosphere, which is, in turn, the result of anthropogenic emissions (IPCC, 2021). The largest share of these emissions (86%) comes from fossil fuel CO₂ emissions (9.6 ± 0.5 Pg C yr⁻¹) (Canadell et al., 2021). While stopping these emissions should remain the first objective, mitigating climate change will require decreasing all GHG sources. Agriculture, forestry and other land use (AFOLU) is a significant source of GHG (12.0 ± 2.9 GtCO₂eq yr⁻¹) (Jia et al., 2019). However, it has the potential to remove CO₂ from the atmosphere. Indeed, carbon uptake by vegetation has increased over the past decades, but uncertainties remain about whether this trend will continue (Canadell et al., 2021). Considering these uncertainties, soil carbon sequestration in croplands and grasslands has a considerable potential for removing CO₂ from the atmosphere (0.4–8.6 GtCO₂eq yr⁻¹) (Jia et al., 2019). Measuring and separating the different processes of land-atmosphere carbon flux is crucial to advise and monitor policies and goals effectively. Doing it, however, is not trivial.

The international scientific community is leveraging advanced techniques to produce reliable surface-atmosphere GHG flux monitoring. Eddy Covariance (EC) is praised for directly and continuously measuring surface turbulent fluxes. Since the early measurements, the method has been applied to different gases, including water vapour, CO₂, CH₄ and N₂O (Valentini et al., 1996; Moncrieff et al., 1997; Fowler et al., 1995). Active development of instrumentation and standardization of the methods and networks has made it the reference measurement for terrestrial ecosystem GHG fluxes (Pastorello et al., 2020).

Continuity is essential to compute annual GHG budgets and long-term soil carbon balance. Despite the attempt to have near-continuous observations, a fraction of the observation is either missing or non-reliable, resulting in data gaps. These gaps can sometimes be due to technical reasons but, most importantly, related to under developed turbulence and non-stationarity of the flow, both required to compute reliable fluxes with the standard EC method (Aubinet et al., 2012; Pastorello et al., 2020).

In FLUXNET2015 (Pastorello et al., 2020), with more than 1500 site-years of data worldwide, 60 % of the CO₂ flux is gap-filled. Filling these short gaps is problematic because it can significantly bias the annual GHG budgets (Du et al., 2014; Vekuri et al., 2023).

Efforts have been mobilized to find defensible methods to fill CO₂ flux gaps (Falge et al., 2001) and make it part of the standard post-processing (Wutzler et al., 2018; Pastorello et al., 2020). From various methods, artificial neural networks (ANN) (Papale and Valentini, 2003; Moffat, 2012) and look-up tables (Falge et al., 2001; Reichstein et al., 2005) seem more
promising than processed-based ones (Moffat et al., 2007). Ensuring enough high-quality data for gap-filling is crucial
during night-time, for less stationary surface fluxes (CH₄, N₂O) (Irvin et al., 2021; Mishurov and Kiely, 2011) and
heterogeneous sites (Aubinet et al., 2002).

Methods that resolve surface flux in time and frequency do not require stationarity. Continuous Wavelet Transform (CWT) (Torrence and Compo, 1998) and Discrete Wavelet Transform (DWT) (Mallat, 1989) are two of these methods that decompose a signal in time and per scale (Farge, 1992; Farge and Schneider, 2001). These methods are sought in airborne campaigns when short-time resolution is needed (Strunin and Hiyama, 2004; Mauder et al., 2007; Desjardins et al., 2018; Metzger et al., 2013) and to retrieve outbursts and non-stationary flux (Schaller et al., 2017; Göckede et al., 2019). By not requiring stationarity, they yield high-quality data with fewer gaps, hence, more available data for analysis to feed gap-filling algorithms and further partitioning.

EC measurements provide the net turbulent surface flux, which is often the primary information we are interested in. However, often enough, we need to partition CO₂ flux (Net Ecosystem Exchange, NEE) into gross primary productivity (GPP) and ecosystem respiration (R<sub>eco</sub>). The need comes because they are modelled in surface models, which we desire to calibrate under varying environmental forcing variables, mainly temperature, radiation, soil and air humidity (Duffy et al., 2021).

Standard methods for NEE partitioning involve using photosynthetically non-active periods to estimate R<sub>eco</sub> and further extrapolate it (night-time method) or use a light-response model for GPP (day-time method) (Reichstein et al., 2012). Alternative partitioning methods involve filtering raw EC data based on the correlation sign between CO<sub>2</sub> and H<sub>2</sub>O mixing ratios before calculating the eddy covariance (Thomas et al., 2008; Scanlon and Kustas, 2010; Klosterhalfen et al., 2019; Zahn et al., 2022). These methods are called conditional sampling. The Thomas et al. (2008) method assumes that, during the day time, CO<sub>2</sub> and H<sub>2</sub>O ecosystem fluxes have opposite signs when a parcel of air is coming from the plant crown, which is dominated by photosynthesis and have similar signs when air is coming from the ground where respiration dominates. This method has been developed and evaluated over a range of canopies and showed a good ability to separate plants (net primary production, NPP) from soil sources (soil respiration R<sub>soil</sub>) (Zeeman et al., 2013; Zahn et al., 2022). A similar method was developed by (Scanlon and Albertson, 2001; Scanlon and Sahu, 2008; Scanlon and Kustas, 2010), which uses the same concept and additional considerations on correlations between H<sub>2</sub>O and CO<sub>2</sub> to partition between stomatal and non-stomatal CO<sub>2</sub> and H<sub>2</sub>O fluxes. This method also showed consistent H<sub>2</sub>O and CO<sub>2</sub> partitioned fluxes (Scanlon and Kustas, 2012; Sulman et al., 2016; Wang et al., 2016; Perez-Priego et al., 2018; Rana et al., 2018; Scanlon et al., 2019). A comparison between Thomas et al. (2008) and Scanlon and Kustas (2010) (Klosterhalfen et al., 2019) showed that the latter provided overall larger soil flux components than the former.

This study presents a wavelet-based framework for turbulent flux calculation. It also proposes a new partitioning method inspired by the conditional sampling concept but uses frequency-resolved fluxes to provide GPP and R<sub>eco</sub>. We first evaluate how many gaps are generated by each flux computing method and how these impact half-hourly statistics, annual carbon cycles and standard partitioning. We then evaluate the new partitioning method by comparing it to the standard methods
(Reichstein et al., 2005). We use CO₂ flux measured in two contrasted ecosystems ICOS sites in a single climatic region, with four years (2019-2022) from a deciduous mixed forest site (FR-Fon) and two years (2021-2022) from a crop site (FR-Gri).

## 2. Material and methods

In this work we processed EC data from two ICOS sites. Both sites were treated equally and passed through the same steps (Figure 1).

<figure><img src="image_2.png" /><figcaption>Figure 1. Processing steps in standard and wavelet-based Eddy Covariance in this work.</figcaption></figure>

### 2.1. Site description

The study uses data from two French sites located in the Parisian region, and part of the ICOS network (https://www.icos-cp.eu) and FLUXNET. Climatically, the area can be described as oceanic with mild temperatures (11.2-11.5°C annual mean) and moderately wet (677-700 mm annual precipitation).

The first site, FR-Gri, is a 19-hectare crop site (Loubet et al., 2011), rotating between maize, wheat, barley, and rapeseed with intermediate crops. The measuring system is setup on a short tower that moves from 2 to 4 meters according to the crop growth. From this site, we used data from January 2021 to December 2022, consisting of winter rapeseed until 31 July 2021, followed by winter wheat from 7 October 2021 to 5 July 2022 and barley seeded on 11 October 2022. The second site, FR-Fon, is a deciduous broadleaf mixed forest mainly composed of oak and a dominant height of 25 meters at the age of
100 (Delpierre et al., 2016). The eddy covariance setup is located at 37 meters. From this site, we used data from January 2019 to December 2022.

In both sites, the Eddy Covariance setup consisted of a closed-path infrared gas analyzer (LI-7200; Li-Cor Inc., Lincoln, NE, USA) and a three-dimensional sonic anemometer (Gill HS; Gill Instruments Ltd, Lymington, Hampshire, UK). Both instruments and acquisition setups followed ICOS guidelines and protocols (Sabbatini et al., 2018).

## 2.2. EC flux processing

To compute the atmosphere-biosphere flux, we consider a virtual rectangle box extending from the ground to the location of the eddy-covariance setup of width W, length L and height h_m. The mass balance of a scalar in the virtual box is used to retrieve the expression of the overall ecosystem flux F_eco (g m⁻² s⁻¹). The mass balance includes a storage term (I), an advection transport term (II) and a turbulent diffusion term (III), which, when integrated over the three dimensions of the virtual box, equals (Foken et al., 2012; Metzger, 2018; Aubinet et al., 2005; van Gorsel et al., 2009):

$$0 = \int_{-L}^{L} \int_{-W}^{W} \int_{0}^{h_m} \left[ -S + \underbrace{\overline{\rho_d \frac{\partial \chi_s}{\partial t}}}_{i} + \underbrace{\overline{\rho_d u \frac{\partial \chi_s}{\partial x} + \rho_d v \frac{\partial \chi_s}{\partial y} + \rho_d w \frac{\partial \chi_s}{\partial z}}}_{II} + \underbrace{\frac{\partial \overline{\rho_d u' \chi'_s}}{\partial x} + \frac{\partial \overline{\rho_d v' \chi'_s}}{\partial y} + \frac{\partial \overline{\rho_d w' \chi'_s}}{\partial z}}_{III} \right] dz \, dx \, dy \quad (1)$$

Where S is the ecosystem volumetric flux (g m⁻³ s⁻¹), ρ_d the dry air density, χ_s the scalar dry mole fraction (mol mol⁻¹), t the time (s), while u, v and w are the upwind, crosswind and vertical component of the windspeed (m s⁻¹). Overbars indicate time averaging; quotation marks the instantaneous deviation from the mean. Assuming a horizontally homogeneous ecosystem (homogeneity in ecosystem functioning and structure over x and y) allows suppressing the horizontal derivatives and integrals in eq. (1) (Finnigan et al., 2003; Metzger, 2018). This assumption also leads to a zero dry air vertical flux due to continuity ρ_d w = 0 (Webb et al., 1980). Then, recognising that the integral of S over the height z is F_eco - F_soil and that similarly the integral of ∂ρ_d w'χ'_s / ∂z over z is ρ_d w'χ'_s - F_soil, eq. (1) leads to:

$$F_{eco} = \overline{\rho_d} \cdot \left( \int_{0}^{h_m} \frac{\partial \overline{\chi_s}}{\partial t} dz + \overline{w' \chi'_s} \right) \quad (2)$$

Where w'χ'_s is the turbulent flux at h_m. The ecosystem flux can hence be evaluated from w'χ'_s and the storage term (∫_0^h_m ∂w'χ'_s / ∂t dz), which may be significant in medium and tall towers but can be neglected in small ones. In practice, w'χ'_s is computed from a time series of w and χ_s and sampled at a frequency typically higher than 5 Hz to capture the smallest eddies contributing to the flux (Gu et al., 2012).

In this study, we use three methods to evaluate w'χ'_s: the standard eddy covariance method (EC_s), and two frequency resolved methods, one using a Continuous Wavelet Transform (CW-EC) and the other using a Discrete Wavelet Transform (DW-EC).
From $w$ and $\chi_s$, the standard method consists of calculating the product of the instantaneous deviations of both variables from their respective means (covariance); frequency-resolved methods work the same, using, however, a priory decomposed instantaneous deviation (Figure 2).

<figure><img src="image_3.png" /><figcaption>Figure 2: Conceptual scheme showing the main processing steps: data pre-processing (1); covariance calculation consisting of Reynolds decomposition and product of instantaneous deviation (2) and using frequency decomposition (3); post-processing (4). $w$ is the vertical component of the wind velocity, $\chi$ is the mixing ratio, $j$ represents the frequency scale, bars are for averaging and quotation marks are deviations from the mean.</figcaption></figure>

### 2.2.1. Data pre-processing

In this study, data flux pre-processing was done using EddyPro 7.0.9. For all flux calculations, we applied de-spiking (Mauder et al., 2013), covariance maximization for time lag correction, and double rotation for tilt correction (Wilczak et al., 2001). The time lag default was set to 0.08 s for a 71.1 cm tube with a 5.3 mm inner diameter and 15 L/min flow rate and was allowed to vary in its vicinity. Using closed path systems and dry mixing ratios for gas avoids compensating for density fluctuations (Kowalski and Serrano-Ortiz, 2007). No detrending was applied. Standard flux calculations require further corrections to address low and high-frequency losses (Massman and Lee, 2002). These spectral corrections are usually applied after flux calculation. Here, for simplicity's sake and since we are focusing on comparing the flux calculation, gap filling and partitioning methods rather than interpreting the fluxes themselves, these corrections were omitted.

### 2.2.2. Standard Eddy covariance (EC<sub>s</sub>)

The eddy covariance method consists of calculating the covariance $\overline{w'\chi_s'}$ (step 2 in Figure 2). Typically, the time average for fluxes is 30 minutes to 1 hour (Rebmann et al., 2018; Pastorello et al., 2020; Aubinet et al., 1999). Turbulent fluctuations $\chi_s'$ and $w'$ are formally defined as deviations from an ensemble average and not from a time-average. The ergodic assumption is required to make the ensemble and time average equivalent. In sum, the averaging period should thus be
stationary and sufficiently long to gather enough data to get a low random error (Lumley and Panofsky, 1964; Kaimal and Finnigan, 1994). In general, a 30-minute period satisfies these requirements.

### 2.2.3. Wavelet Transform methods

Wavelet transform is a bandpass filter allowing decomposition of a time series into sub-series defined for a given frequency. The following steps explain how to perform a frequency-resolved covariance using wavelets (Figure 2, panel 3). More details can be found in (Farge, 1992; Torrence and Compo, 1998; Farge and Schneider, 2001).

Any signal $f(t)$ can be decomposed into different scales, which results in the signal itself once added up. The simplest example is the Reynolds decomposition that separates a time series into its mean and its instantaneous deviation:

$$f(t) = f'(t) + \overline{f(t)} \qquad (3)$$

In Eq. 3, the mean, $\overline{f(t)}$, is the low-frequency component, with a frequency representative of $1/T$, where T is the averaging period. Similarly, a time series can be decomposed into J sub-series, each representative of a frequency j:

$$f(t) = \sum_{j=0}^{J} f(t, j) \qquad (4)$$

The wavelet transform is a way to decompose the signal using a mother wavelet $\psi$, a wave function with finite support (Eq. 5). Unlike removing an average, wavelet transform yields stationary sub-series (Torrence and Compo, 1998). Considering N discrete observations with a sampling period $\delta_t$, so that $t = n \delta_t$ where n is the time index, we can generate a family of wavelets normalized in L²-norm:

$$\psi_{n,j}(n') = s_j^{-1/2} \psi \left[ \frac{(n' - n)\delta_t}{s_j} \right] \qquad (5)$$

Where $s_j$ is the scaling factor, usually defined using a geometric progression with a maximum limited by the total sampling period $N\delta_t : s_j = s_0 2^{j\delta_j}$, for $j = 0, 1, ... J$. Here, J is the size of the set of scales, $s_0$ is the smallest resolvable scale, approximately $2\delta_t$, and $\delta_j$ is the scale factor. The convolution of the signal $f(n)$ with a scaled mother wavelet $\psi$, yields the $\tilde{f}(n, j)$, the wavelet coefficient for time series $f(n)$:

$$\tilde{f}(n, j) = \sum_{n'=0}^{N-1} f(n') \psi_{n,j}(n') \qquad (6)$$

From which one can reconstruct the signal:
$$f(n) = \frac{\delta j \delta t^{0.5}}{C_{\delta} \psi_{0}(0)} \sum_{j=0}^{\infty} \Re\{\tilde{f}(n, j)\} \qquad (7)$$

$C_{\delta}$ is a scale-independent reconstruction factor depending on the chosen mother wavelet function. Note that $f(n)$ is the sum of all its components, however at any specific frequency the wavelet transform works such that the decomposed signal $f(n, j)$ averages to zero and should be interpreted as its instantaneous deviation. Note that $C_{\delta}$ is only required for continuous wavelet decompositions in which the wavelet function is not an orthogonal base. Although this allows an arbitrary set of scales to be chosen, providing a more resolved signal spectrum (Arts and van den Broek, 2022; Torrence and Compo, 1998), a correction factor $C_{\delta}$ is required to compensate for these overlaps, not required in discrete wavelet decomposition. Applying eqns (5 and 6) to $w$ and $\chi_s$ we can calculate the total flux as:

$$F_c = \overline{\rho_d \overline{w' \chi_s'}} = \overline{\rho_d C_{\varphi}} \sum_j \overline{w' \chi_s'}(j) \qquad (8)$$

Where $\overline{w'\chi_s'}(j)$ is the mean of the product between $w$ and $\chi_s$ decomposed signals at frequency $j$ (frequency-resolved covariance), $C_{\varphi}$ is the reconstruction factor depending on the chosen mother wavelet function and determined empirically by comparing $\overline{w'\chi_s'}(j)$ to $\overline{w'\chi_s'}$, where averaging is done over the time index $n$. When continuous wavelet decomposition is used, since it is not an orthogonal base, the sum of the $\overline{w'\chi_s'}(j)$ for all $j$ is not strictly equal to $\overline{w'\chi_s'}$ as cross-correlations between scales are not zero. The empirical factor $C_{\varphi}$ is used to correct for this effect (see Supp. Mat. A). When discrete decomposition is used, $C_{\varphi} = 1$, since the orthogonality of the wavelet functions, which characterizes DWT, implies that total energy is conserved and yields independent frequencies, hence cross-correlations between scales is zero. The orthogonal base also forces the set of scales to be discrete, $\delta_j = 1$ and $s_j = 2^j$, for $j = 0,1, ... J$ (Farge, 1992). The wavelet coefficients are then:

$$\psi_{n,j}(n') = 2^{j/2} \psi(2^j n' - n) \qquad (9)$$

A great interest in the DWT is that the orthogonality and progressively smaller decomposition make it far cheaper computationally than CWT at the expense of a coarser resolution in frequency (Mallat, 1989), making it a good candidate for time series longer than a couple of weeks without significant difference (see Figure S 1 for a brief comparison).

Commonly used wavelet functions are the Morlet and the Mexican Hat for continuous decomposition for they are well-defined in the frequency and time domain (Schaller et al., 2017), and the Daubechies 6 for the discrete decomposition (Table 1).

<table><caption>**Table 1:** Mother wavelets used in this study. Mother wavelet formula, $\psi(\eta)$, empirically derived factors, $C_{\delta}$ and $\psi_0(0)$, from (Farge, 1992) and $C_{\varphi}$ (see Figure S 1).</caption><thead><tr><td>Name</td><td>Decomposition</td><td>ψ(η)</td><td>C<sub>δ</sub></td><td>ψ<sub>0</sub>(0)</td><td>C<sub>φ</sub></td></tr></thead></table>
<table><tr><td>Morlet (k<sub>ψ</sub>= 6)</td><td>Continuous</td><td>e<sup>ik<sub>ψ</sub>η</sup>e<sup>-η<sup>2</sup>/2</sup></td><td>0.776</td><td>π<sup>-1/4</sup></td><td>5.271</td></tr><tr><td>DOG (m=2), a.k.a. Marr or Mexican Hat</td><td>Continuous</td><td>-1<sup>m</sup> d<sup>m</sup>/dη<sup>m</sup>(e<sup>-η<sup>2</sup>/2</sup>)</td><td>3.541</td><td>0.867</td><td>16.568</td></tr><tr><td>Daubechies (k=6)</td><td>Discrete</td><td>(-1)<sup>k</sup>a<sub>N-1-k</sub></td><td>1</td><td>1</td><td>1</td></tr></table>

In this study, we used the discrete decomposition and the Daubechies (k=6) wavelet (Daubechies, 1988). For comparisons with the standard eddy covariance method we compute the covariance by summing scales $s_j$ smaller than 1800 seconds (30 min) in eq. (8). Calculations were done using *PyWavelets* module (Lee et al., 2019). Despiking (Mauder et al., 2013) was used on each sub-series to eliminate any unrealistic values identified.

### 2.2.4. Cone of influence

Wavelet coefficients calculated with the convolution product in eq. (6) are subject to the influence of neighbours, resulting in a time “influence cone” that grows with decreasing frequency (Torrence and Compo, 1998). This cone renders the reconstruction unusable at the edges of the dataset and for scaling factors close to the dataset duration. The cone of influence (COI) is the boundary of the wavelet spectrum, which, exterior to its edge, effects become important. It is defined as:

$$COI = f \delta t n, n = 0, 1, ..., N-1, N, N+1, ..., 1,0 \qquad (10)$$

f is the Fourier factor specific to each wavelet.

We extended the dataset over periods larger than the period of interest for every averaging time to avoid the cone of influence.

## 2.3. Timeseries flagging and gap-filling

Previous steps allowed us to calculate $\overline{w'c'}$. We still need to verify the EC's assumptions through a quality check (Figure 1). Non-stationarity data for standard EC and periods with a lack of turbulence for both standard and wavelet-based EC are considered unreliable and thus flagged and further gap-filled.

### 2.3.1 Quality flags

Quality flags followed the standard 0-1-2 flag system used in FLUXNET (Mauder and Foken, 2011). The system is based on two tests, one for stationarity and another to verify that turbulence is fully developed (Foken and Wichura, 1996). It is important to recall that standard EC cannot be used during non-stationary moments, but wavelet decomposition yields stationary sub-series that allow skipping this step.

The stationarity test (STA) calculates the absolute relative deviation between the mean of the covariances computed over 5-min intervals and the covariance computed over a 30-min period:
$$STA = \left| \frac{\frac{1}{6} \sum_{i=1}^{6} \frac{\overline{(w'c')^5}_{i} - \overline{(w'c')^{30-mn}}}{\overline{(w'c')^{30-mn}}} \right| \quad (11)$$

The turbulence test, or integral turbulence characteristics (ITC) test, identifies if eddies are fully developed by calculating the absolute relative deviation between the measured and modelled integral turbulent characteristic $\sigma_w/u_*$.

$$ (\sigma_w/u_*)^{\text{model}} = \begin{cases} 0.21 \ln\left(\frac{z+f}{u_*}\right) + 3.1, & \text{if } -0.2 < z/L < 0.4 \\ 2\left(\frac{z}{L}\right)^{1/8}, & \text{else} \end{cases} \quad (12) $$

Where $f$ is Coriolis parameter (s⁻¹), $u_*$ friction velocity (m/s), $z$ is the height (m), $L$ Obukhov length (-), and $z+$ is set to 1 meter for mathematical convention so that $\frac{z+f}{u_*}$ is dimensionless (Thomas and Foken, 2002).

$$ ITC = \left| \frac{(\sigma_w/u_*)^{\text{model}} - (\sigma_w/u_*)^{\text{measurement}}}{(\sigma_w/u_*)^{\text{measurement}}} \right| \quad (13) $$

A detailed description of the quality control procedures can be found in Foken and Wichura (1996) and Mauder and Foken (2011). Data is considered high-quality when this deviation is below 30% for all applicable tests, as shown in Table 2.

<table><caption>**Table 2:** Quality flag categories proposed by Mauder and Foken (2011), based on stationarity and integral turbulence characteristics (ITC) tests presented by Foken and Wichura (1996).</caption><thead><tr><th>Quality flag</th><th>Stationarity test (STA)</th><th>Integral turbulence characteristics test (ITC)</th></tr></thead><tbody><tr><td>0 (High)</td><td>&lt; 30 %</td><td>&lt; 30 %</td></tr><tr><td>1 (Medium)</td><td>31 – 100 %</td><td>31 – 100 %</td></tr><tr><td>2 (Low)</td><td>&gt; 100 %</td><td>&gt; 100 %</td></tr></tbody></table>

### 2.3.2. u* filtering

Further screening is necessary to discard observations below a friction velocity threshold ($u_{*crit}$) (Wutzler et al., 2018; Papale et al., 2006). Under stable stratified atmospheric conditions, the EC technique has been shown to underestimate nocturnal CO₂ respiration (Goulden et al., 1996; Baldocchi, 2003). The reason is that the turbulence is attenuated by the positive air density gradient (Kaimal and Finnigan, 1994). As biotic flux is not expected to depend on turbulence, we can define a threshold value for friction velocity ($u_{*crit}$) below which the measured ecosystem CO₂ flux starts to decrease. Below $u_{*crit}$ turbulence is not developed enough to mix the surface layer and the EC to perform well. This method provides an alternative way to determine the turbulent requirement based on an ecosystem function instead of using a physical-based as with ITC.
Once a threshold is defined, observations below this threshold are dropped and gap-filled (Gu et al., 2005; Aubinet et al., 2012).

The u\* threshold was determined using the REddyProc library in R and was free to vary among seasons (Wutzler et al., 2018).

### 2.3.3. Gap-filling

Gap-filling was performed on data flagged for medium (1) and low (2) quality or below u\*crit for ECs (stationary and turbulence flag considered) and DW-EC (only turbulence flag considered). We used the Marginal Distribution Sampling (MDS) method, the most commonly used gap-filling method (Pastorello et al., 2020). MDS consists of sampling data in the temporal vicinity of the data to be gap filled, usually a 15-day window, with similar meteorological conditions defined by the income shortwave, the air temperature, and the vapour pressure deficit. This subset yields a distribution function used to fill the gap, exploiting both the meteorological drivers and the temporal auto-correlation structure of NEE (Reichstein et al., 2005). For the calculations, we used the REddyProc library in R (Wutzler et al., 2018).

## 2.4. NEE partitioning

Flux partitioning refers to the division of the Net Ecosystem Exchange (NEE) into the ecosystem respiration ($R_{eco}$) and the gross primary productivity (GPP). Ecosystem respiration refers to the release of CO₂ by organisms during their metabolic activities, including autotrophic respiration by plants and heterotrophic respiration by micro- and macro-organisms in soil and the ecosystem. GPP represents the uptake of CO₂ by plants through photosynthesis:

$$NEE = GPP + R_{eco} \qquad (14)$$

GPP is a flux directed from the atmosphere to the ground (negative), while $R_{eco}$ is from the ground to the atmosphere (positive). In standard practice, partitioning relies on the presumed responses of GPP and $R_{eco}$ to light, water, and temperature. We applied the known night- and day-time methods on both standard and wavelet-based CO₂ fluxes, and propose here a new method for the wavelet-based flux.

### 2.4.1. Night-time partitioning method

Night-time (NT) partitioning assumes that GPP is zero at night, so NEE equals $R_{eco}$ (Reichstein et al., 2005). Reference respiration rate is then parametrized using an Arrhenius-type temperature response model for nocturnal measurements and projected into the day (Lloyd and Taylor, 1994).

$$R_{eco} = R_{ref} \cdot e^{E_0 \left( \frac{1}{T_{ref} - T_0} - \frac{1}{T_{air} - T_0} \right)} \qquad (15)$$
Where $R_{ref}$ (µmol·m⁻²·s⁻¹) is the reference respiration rate at the reference temperature ($T_{ref}$ = 15°C), $T_{air}$ is air temperature, $T_0$ is fixed at -46.02°C, $E_0$ (°C) is the temperature sensitivity, a free parameter. A constant value is estimated for $E_0$ for the whole year, while $R_{ref}$ is estimated every five days using a 15-day window (Reichstein et al., 2005).

Further references to NT estimations use the terms NT-GPP and NT-$R_{eco}$. The R code implementation for NT is available to download from https://github.com/bgctw/REddyProc (Wutzler et al., 2018).

### 2.4.2. Day-time partitioning method

Day-time (DT) partitioning differs from NT in that a light response curve (Lasslop et al., 2010) is parametrized using day-time measurements. NEE is estimated as follows:

$$NEE = \frac{\alpha\beta R_g}{\alpha R_g + \beta} + R_{eco} \qquad (16)$$

Where $R_{eco}$ is a respiration model eq. (15), $R_g$ is the global radiation (Wm⁻²), $\alpha$ (µmol·C·J⁻¹) is the initial slope of the light response curve, and $\beta$ (µmol·m⁻²·s⁻¹) is the maximum rate of CO₂ uptake of the canopy at light saturation. $\beta$ is estimated using an exponentially decreasing function of atmospheric vapour pressure deficit of air (VPD):

$$\beta = \begin{cases} \beta_0 e^{-k(VPD-VPD_0)}, & VPD > VPD_0 \\ \beta_0, & VPD < VPD_0 \end{cases} \qquad (17)$$

Note that what is physiologically more relevant in $\beta$ is the leaf-to-air VPD, which can vary from atmospheric VPD in the same direction as leaf temperatures vary from air temperature. However, flux sites measure atmospheric rather than leaf-to-air VPD.

The standard calibration procedure is done in two steps. First, $E_0$ and $R_{ref}$ are estimated using night-time observations. The remaining parameters ($\alpha$, $\beta_0$, k, and $VPD_0$) and $R_{ref}$ (now using previous estimation as a prior) are fitted using Eq. (16) on day-time data.

Recent studies have proposed the inhibition of leaf respiration in the light as a source of mismatch between EC and independent $R_{eco}$ measurements (Wehr et al., 2016). A modified version of standard partitioning has been proposed to include this mechanism (Keenan et al., 2019). The modified DT version preserves the structure of the original DT but uses $R_{ref}$ prior, fitted during night-time, for nocturnal partitioning while estimating daytime as usual. The R code implementation for DT's original and modified versions can be downloaded from https://github.com/bgctw/REddyProc (Wutzler et al., 2018). Unless specified otherwise, DT estimations follow the modified version (Keenan et al., 2019) and are referred to as DT-GPP and DT-$R_{eco}$.

### 2.4.3. A new direct wavelet-based partitioning method

Direct observations of gross primary productivity (GPP) and ecosystem respiration ($R_{eco}$) are not feasible at the field scale, thus justifying the necessity of model-based partitioning. Thomas et al. (2008) and Scanlon and Kustas (2010) have proposed
ingenious ways of incorporating prior knowledge of co-processes to compute a model-free partitioning of soil respiration ($R_{soil}$) and plant net primary productivity (NPP).

In this study, we take advantage of the frequency decomposition of $w'\chi_s'$ using wavelets to go beyond what was done by Thomas et al. (2008), thus the name discrete-wavelet based conditional sampling (DW-CS). In an empirical approach we conditionally sample the frequency decomposed products $w'(t, j) * \chi_s'(t, j)$ to separate positive and negative components of the CO₂ and H₂O fluxes in each frequency j. The underlying empirical assumption is that wavelet decomposition should allow to trap in each frequency the positive and negative “gusts” which are mixed up in the original signal. For simplicity $w'\chi_s'(t, j)$ is replaced by $\phi_s$ in the following equations:

$$ \phi_c = \phi_c^+ + \phi_c^- = \phi_c^+|\phi_v^+ + \phi_c^+|\phi_v^- + \phi_c^-|\phi_v^+ + \phi_c^-|\phi_v^- \quad (18) $$

Where c and v stand for CO₂ and H₂O, x⁺ stands for sampling x when x is positive and the opposite for x⁻, and x|y stands for sampling x when y is true. We could assume that positive CO₂ flux ($\phi_c^+$) is $R_{eco}$ and negative ($\phi_c^-$) is GPP. However, to guarantee physical meaning of GPP we took advantage of GPP's dependency to light, more precisely photosynthetic photon flux density (PPFD), and set GPP to zero during night ($PPFD \le 10 \mu mol m^{-2} s^{-1}$). We further considered that the daytime ($PPFD > 10 \mu mol m^{-2} s^{-1}$) negative CO₂ fluxes conditioned by negative water vapour fluxes ($\phi_c^-|\phi_v^-$) as non-realistic and therefore attributed it equally to $R_{eco}$ and GPP (see Figure 3), which leads to the following definition or $R_{eco}$ and GPP:

$$ \begin{array}{l} PPFD \le 10 \mu mol m^{-2} s^{-1} \left\{ \begin{array}{l} Reco = \phi_c^+ + \phi_c^- \\ GPP = 0 \end{array} \right. \hfill (19) \\ PPFD > 10 \mu mol m^{-2} s^{-1} \left\{ \begin{array}{l} Reco = \phi_c^+ + 0.5 \times \phi_c^- |\phi_v^- \\ GPP = \phi_c^- |\phi_v^+ + 0.5 \times \phi_c^- |\phi_v^- \end{array} \right. \end{array} $$

<figure><img src="image_4.png" /><figcaption>Figure 3: Conceptual scheme for wavelet-based NEE flux partitioning. $\phi_c$ stands for $w'\chi_{CO_2}'(t, j)$ and $\phi_v$ for $w'\chi_{H_2O}'(t, j)$. Quadrants and arrows in the figure show conceivable fluxes during day and night. In quadrants, grey arrows show reallocation from unlikable (question mark) and unreasonable (“x”) fluxes towards the most probable actual flux.</figcaption></figure>
## 2.5. Performance measurements

Comparisons between methods were carried out using mean bias and the annual gap-filled CO2 flux balance error. Defined as:

$$Mean Error (bias) = \frac{1}{N} \sum_{n=1}^{N} (NEE_{x,n} - NEE_{y,n}) \qquad (20)$$

$$Mean Absolute Error = \frac{1}{N} \sum_{n=1}^{N} |NEE_{x,n} - NEE_{y,n}| \qquad (21)$$

Where N equals the amount of data, NEE$_{x,n}$ is the Net Ecosystem Exchange calculated using one of the x methods among EC$_{S}$ and DW-EC at a time n.

# 3. Results

## 3.1. CO$_{2}$ flux computed by EC and DW-EC

In this section we compare the CO$_{2}$ flux (or NEE) computed by EC$_{S}$ and DW-EC, and then analyse the additional information on DW-EC's cospectra. During the photosynthetically active months (warmer months for FR-Fon and crop seasons for FR-Gri), the two sites were carbon sinks with a negative NEE (Figure 4). During winter, when the trees have lost their leaves and crop sites are bare soil, the lack of GPP transforms sites into sources with a positive NEE. As a consequence, NEE in FR-Fon showed a clear seasonality, while FR-Gri showed a more variable pattern. We observe a substantial decrease in absolute value in the NEE for short periods during summers and springs for all years. Some relate to cloudy days, others to high vapour pressure deficit, which indicates air dryness. In June and July 2019, France was hit by short heat waves (Sousa et al., 2020; Pohl et al., 2023). In the crop site (FR-Gri), we identify the crop season in the spring of both years and the intercrop in Autumn 2021 by the decrease in NEE. We can also notice that harvest is done long after NEE has become positive; this is to bring the crop to maturity after senescence. Finally, the crop site also showed an earlier growth compared to the forest site, which is expected as the trees at this site are 100 years old on average and have, therefore, a late foliar development during the year while the crops were winter crops, which are in their growth stage early during the year.

Overall, the NEE ranged from -10 to 6 µmol m-2 s-1, with stronger respiration during winter and spring at the crop site compared to the forest site. Daily mean NEE estimated by EC$_{S}$ and DW-EC were very close to each other (R$^2$= 0.97 (0.98), ME= 0.1 (0.05) µmol m-2 s-1, MAE= 0.33 (0.38) µmol m-2 s-1, ECS=1.08 (1.12) ×DW linear fit for FR-Fon (FR-Gri)).
<figure><img src="image_5.png" /><figcaption>Figure 4: (a) NEE cospectra derived by DW-EC and averaged half-hourly. Colours indicate NEE cospectra, and grey for missing data. (b) Average for NEE cospectra (black), exclusively negative (blue) and positive (red) values of NEE. (c) Daily average NEE computed from the NEE cospectra integrated over 30-minutes⁻¹ (black, DW) and ECₛ (grey).</figcaption></figure>

Looking at NEE's cospectra, we can see a peak around 6 seconds⁻¹ (0.16 Hz) frequency in FR-Gri and around 50 seconds⁻¹ (0.02 Hz) in FR-Fon (Figure 4 b). The peak frequency is related to the measurement height, being higher on the 2 to 4 meters tower in FR-Gri than on the 37 meters tower FR-Fon (around 10 meters above the canopy). The measurement height affects the frequency contribution to the CO₂ flux because the height above ground constrains the size of the eddies. Indeed, the maximum cospectral frequency is linearly dependent on height and modulated by the wind speed and the stability parameter z/L, where L is the Obukhov length (Kaimal and Finnigan, 1994). Indeed, the cospectra calculated using equations in (Horst, 1997) peaks at the same 50 seconds⁻¹ (6 seconds⁻¹) frequency for FR-Fon (FR-Gri). A secondary maximum at around 30 minutes⁻¹ can be seen in the positive and negative NEE cospectra but disappears on the NEE. The cospectra of NEE's positive (φ_c⁺) and negative (φ_c⁻) counterparts were overall mirrored with, however, slight differences: a higher contribution of higher frequencies on φ_c⁺ and of lower frequencies on φ_c⁻ in FR-Fon, suggesting large coherent structures may contribute more to GPP (defined mainly by φ_c⁻) than R_eco on average. This was not observed at the crop site.

Seasonally, NEE (φ_c) cospectra aligned with theoretical estimations, displaying peaks near the expected frequencies (Figure S 3). During months with high carbon sequestration (from April to October in FR-Fon and February to June in FR-Gri) and under neutral or unstable stratification, the negative (φc-) portion of cospectra exhibited lower-frequency peaks compared to the positive (φc+) portion. Conversely, stable conditions prompted the opposite pattern, albeit less pronounced in the crop site due to an unexpected secondary peak around 30 minutes⁻¹, which softens the distribution on the higher frequencies. This secondary peak became most evident from June (stable) and July (neutral and unstable) through November, coinciding with
the post-harvest period after the primary crop cycle and before Winter sets in. No similar seasonality was seen in the theoretical cospectra, which considers micrometeorology conditions.

Still looking at raw data, both methods show a clear daily and seasonal pattern for NEE expected for these ecosystems (Figure 5). Indeed, half-hourly DW- and EC<sub>s</sub>-NEE were very close to each other (R<sup>2</sup>= 0.98, Bias=0.14 µmol m<sup>-2</sup> s<sup>-1</sup>, MAE=0.58 µmol m<sup>-2</sup> s<sup>-1</sup>, EC<sub>s</sub>=1.08×DW linear fit, sites combined) when both were high-quality data, deteriorating when moving to medium and low-quality data (Figure S 2). In the forest site, during March and April, we can see peaks in EC<sub>s</sub>-NEE 5<sup>th</sup> and 95<sup>th</sup> percentile which are lower in DW. These two months had the highest non-stationarity in the site, yielding 60% of the observations unreliable. The same ECS-NEE peaks are seen for the crop site but are less closely related to the stationarity flag. They could be related to night-time CO<sub>2</sub> spikes due to advection from the nearby animal barns (around 600 m west).

365

370

<figure><img src="image_6.png" /><figcaption>Figure 5. Half-hourly NEE estimated monthly using EC<sub>s</sub> (blue) and DW (orange). The darker region indicates interquartile (25<sup>th</sup> and 75<sup>th</sup> percentile), and the lighter region with dotted lines indicates the 5<sup>th</sup> and 95<sup>th</sup> percentile. Below the curves, the monthly statistics are shown: the percentage of non-stationarity (STA>0) and low turbulence (ITC>0) data, the correlation coefficient (R<sup>2</sup>), the mean error (ME, µmol m<sup>-2</sup> s<sup>-1</sup>), the mean absolute error (MAE, µmol m<sup>-2</sup> s<sup>-1</sup>) and the linear fit.</figcaption></figure>

375

## 3.2. Effects of flux processing method on data coverage

In this section we analyse the effect of flux processing methods on the number of data gaps. The quality control steps related to turbulence, namely the Integral Turbulence Characteristic (ITC) test and friction velocity threshold (u\*crit), unsurprisingly discarded the most during the night (Rg < 10 Wm<sup>2</sup>) on both the forest (FR-Fon) and the crop site (FR-Gri) (Table 3). This is

380
expected since the surface cools during these periods, creating a stable layer near the surface and preventing turbulent mixing. $u_{*crit}$ impacted more DW-EC than ECs because $u_{*crit}$ is calculated at the end of the quality control; thus, for ECs, all periods with co-occurrence between non-stationarity and low $u_{*crit}$ ($u_{*crit} < u_{*crit}$) had already been dropped.

The non-stationarity test flagged a significant amount of data during the day and night. However, considering the co-occurrence between flags, it impacted more day-time observations when turbulence is usually well-developed. The difference in the total amount of discarded data between methods reflects this. During night time, ECs discarded around 17 % more than DW-EC and 20 to 30 % more during day-time (Table 3).

<table><caption>**Table 3:** Cumulative percentages of discarded data at each quality control step for the CO₂ flux at the two sites. Medium and low-quality data are replaced. Symbol * means less than 1%.</caption><thead><tr><th rowspan="2">Quality control step</th><th colspan="2">FR-Gri (2021-22)</th><th colspan="2">FR-Fon (2019-22)</th></tr><tr><th>Night</th><th>Day</th><th>Night</th><th>Day</th></tr></thead><tbody><tr><td>Missing data</td><td>5 %</td><td>5 %</td><td>10 %</td><td>10 %</td></tr><tr><td>Turbulence not fully developed (ITC medium or low)</td><td>28 %</td><td>14 %</td><td>17 %</td><td>6 %</td></tr><tr><td colspan="5">Discrete Wavelet Transform (DW-EC)</td></tr><tr><td>Friction velocity threshold ($u_{*crit}$ = 5<sup>th</sup>/50<sup>th</sup>/95<sup>th</sup> percentile)</td><td>10 / 18 / 28 %</td><td>1 / 2 / 5 %</td><td>12 / 20 / 30 %</td><td>2 / 5 / 10 %</td></tr><tr><td><strong>Total</strong></td><td><strong>49 %</strong></td><td><strong>20 %</strong></td><td><strong>42 %</strong></td><td><strong>19 %</strong></td></tr><tr><td colspan="5">Standard (ECs)</td></tr><tr><td>Non-stationarity (STA medium or low)</td><td>19 %</td><td>20 %</td><td>23 %</td><td>30 %</td></tr><tr><td>Friction velocity threshold ($u_{*crit}$ = 5<sup>th</sup>/50<sup>th</sup>/95<sup>th</sup> percentile)</td><td>6 / 13 / 21 %</td><td>*/* / 3 %</td><td>5 / 11 / 19 %</td><td>*/ 2 / 6 %</td></tr><tr><td><strong>Total</strong></td><td><strong>66 %</strong></td><td><strong>40 %</strong></td><td><strong>60 %</strong></td><td><strong>48 %</strong></td></tr><tr><td><strong>Difference between DW-EC and ECs</strong></td><td><strong>17 %</strong></td><td><strong>20 %</strong></td><td><strong>18 %</strong></td><td><strong>29 %</strong></td></tr></tbody></table>

When looking at the length of the gaps (Table 4), we found a decrease in all gap lengths when using DW-EC compared to EC. With more observations and narrower gaps, we expect that the DW-EC method would improve the accuracy of any commonly used gap-filling methods (Moffat et al., 2007), improving the annual NEE accuracy.

<table><caption>**Table 4:** Occurrence of gaps by length for each site identified by EC and DW-EC. One occurrence is a period of 1 or many points of gap. In parentheses: percentage of the data concerned by the gap length over total data length.</caption><thead><tr><th rowspan="2">Gap length (record number)</th><th colspan="2">FR-Gri</th><th colspan="2">FR-Fon</th></tr><tr><th>ECs</th><th>DW-EC</th><th>ECs</th><th>DW-EC</th></tr></thead><tbody><tr><td>1-2</td><td>2538 (9 %)</td><td>2084 (7 %)</td><td>5076 (9 %)</td><td>2192 (4 %)</td></tr><tr><td>3-5</td><td>592 (6 %)</td><td>424 (4 %)</td><td>1373 (7 %)</td><td>600 (3 %)</td></tr></tbody></table>
### 3.3 Effects of flux processing method on gap filling

In this section, NEE measurements using the standard EC ($EC_s$-NEE) were compared with the DW-EC (DW-NEE) method to assess degree of their agreement, potential biases, and the reliability of the DW-EC method. Both methods were gap-filled using MDS; however, the gap-filling was performed on a different number of data as shown by the quality control filtering excluding more data in ECs than in DW (Table 3). Daily gap-filled DW- and $EC_s$-NEE agreed well (Figure 6.b), with only marginal differences from before gap-filling, suggesting gap-filling was unbiased over a day on these sites. Unexpectedly, gap-filled NEE in the forest site had a MAE 0.1 µmol m⁻² s⁻¹ higher than raw NEE. To understand the increase in MAE, we calculated the RMSE, which showed a decrease from 3.9 µmol m⁻² s⁻¹ (raw) to 1.9 µmol m⁻² s⁻¹ (gap-filled), suggesting MAE distribution got narrower, also confirmed by looking at the quartiles (not shown here). Despite these differences, on both sites, MAE was of the same order of magnitude as the mean random uncertainty (1.03 µmol m⁻² s⁻¹ in FR-Fon; 0.73 µmol m⁻² s⁻¹ in FR-Gri).

<figure><img src="image_7.png" /><figcaption>Figure 6. (a) Daily averaged NEE was calculated using $EC_s$ (blue) and DW-EC (orange), both gap-filled with the MDS method. Dotted vertical lines show the start or end of the seasons (years for FR-Fon and crop season for FR-Gri). (b) Daily DW- and $EC_s$-NEE, in grey 1:1 line, in orange true linear relation. On the bottom the correlation coefficient (R²), the mean error (ME, µmol m⁻² s⁻¹), the mean absolute error (MAE, µmol m⁻² s⁻¹) and the linear fit.</figcaption></figure>

Overall, daily DW underestimated $EC_s$ by 9 %, and half-hourly by 4 % (8 %) in FR-Fon (FR-Gri). This is consistent with other flux studies using wavelets (Desjardins et al., 2018; Mauder et al., 2007; Metzger et al., 2013; Schaller et al., 2017), which found underestimations between 3 % and 9 %. Wavelet-based EC way of calculating instant deviation (e.g., w' and $CO_2$') works as a low-frequency filter and detrends the signal instead of simply subtracting the mean as in standard EC. Detrending has been found to lead to an underestimation of around 2 % to 15 %, depending on the running mean filtering used (Rannik and Vesala, 1999), which would explain the observed underestimations.
To disentangle the differences due to gap-filling from those due to the flux computation, we gap-filled the DW-NEE with the DW gaps (DW as previously done), and with the gaps from EC<sub>S</sub> (DW') (Figure S 9). The comparison between DW' and DW yields high correlation (R<sup>2</sup>= 0.97 (0.98) for FR-Fon (FR-Gri)), small bias (ME=-0.01 (0.05) µmol m<sup>-2</sup> s<sup>-1</sup>) and deviation (MAE=0.5 (0.2) µmol m<sup>-2</sup> s<sup>-1</sup>). Based on the small difference between DW (only discards underdeveloped turbulence) and DW' (discards underdeveloped turbulence and non-stationarity) we can conclude that the gap-reduction effect was small for these sites. Accounting for the non-stationary conditions would increase (decrease) the annual NEE by 0.01 (0.05) µmol m<sup>-2</sup> s<sup>-1</sup> or around 1 (2) % in FR-Fon (FR-Gri).

## 3.4 Effects of flux processing method on standard NEE partitioning

In this section we examine whether the gaps in NEE obtained with EC<sub>S</sub>- or DW-NEE have an impact on the partitioning of NEE in GPP and R<sub>eco</sub>. Half-hourly observations show that using EC<sub>S</sub>- or DW-NEE yields similar GPP: R<sup>2</sup>=0.94 (0.97), ME=-0.35 (0.4) µmol m<sup>-2</sup> s<sup>-1</sup>, MAE= 1.0 (0.71) µmol m<sup>-2</sup> s<sup>-1</sup>, ECS=1.07 (1.09)×DW, and similar R<sub>eco</sub>: R<sup>2</sup>=0.69 (0.74), ME=-0.24 (0.33) µmol m<sup>-2</sup> s<sup>-1</sup>, MAE= 0.48 (0.85) µmol m<sup>-2</sup> s<sup>-1</sup>, ECS=1.05 (1.05)×DW, for NT (DT) method (see daily mean statistics in Figure S 4). For all cases, DT yielded higher R<sup>2</sup> than NT; for other statistics, it depended on the site and on which variable was considered (GPP or R<sub>eco</sub>). For instance, NT yielded smaller ME and closer to 1 linear relation than DT in the forest site but the opposite in the crop site. R<sub>eco</sub> estimations using DT method on DW-NEE were higher than on EC<sub>S</sub> (by 5 %), while the opposite is true for all the other cases.

<figure><img src="image_8.png" /><figcaption>Figure 7. GPP and Reco daily average from January 2019 to December 2022 using night- (blue) and day-time (orange) partitioning on EC<sub>S</sub> (dashed) and DW-EC (solid). Note that positive values estimate R<sub>eco</sub> and negative GPP. Dotted vertical lines do not influence the data; they assign the start or end of the season (years for FR-Fon forest and crop season for FR-Gri). Inverted triangles indicate daily temperature mean (red) or cumulated rain higher than the 99th percentile of that year.</figcaption></figure>

We note that differences in DT and NT in FR-Fon in June 2022 happened during a several-weeks-long gap (see Figure 4) and should not be interpreted. While DT estimations show overall good agreement between EC<sub>S</sub> and DW-EC, NT
estimations using DW-NEE yield a smaller GPP and R<sub>eco</sub> than ECS in June 2020 (FR-Fon) and October 2021 (FR-Gri) (Figure 7). On both occasions, NT's R<sub>ref</sub> parameter for EC<sub>s</sub> was at its maximum (Figure S 5). At the end of October 2021, the FR-Gri intermediate crop was harvested, and some residues were left on the field. NT cannot distinguish the different carbon sources and calculates a single temperature response curve. The increase in R<sub>ref</sub> led to an increase in R<sub>eco</sub> exponentially with the warmer day-time temperatures, forcing a physically inconsistent positive GPP. Interestingly, NT estimations using DW-NEE were lower, and R<sub>ref</sub> did seem to follow a smoother seasonal pattern (Figure S 5). In June 2020, a bias could emerge from moving from cloudy and rainy to sunny and warmer.

DT's light-response model, would avoid the positive GPP problem but would produce estimations more sensitive to dynamic day-time conditions. Several occasions where DT estimations are higher than NT coincide with high day-time temperatures (Summer 2019, August 2021) or intense precipitations (June 2021).

**3.5. Evaluation of the new wavelet-based method for direct flux partitioning**

In this section, we compare the new DW-CS method with NT, both calculated using DW-NEE as basis to compare only the partitioning algorithm. NT was chosen as the reference method due to the relative complexity of interpreting DT's variance found in the previous section. However, the results would be similar with DT (see Figure S 6 for a comparison between all methods). Overall, the partitioning methods agreed well (Figure 8b) with a mean absolute daily error of 0.81 (0.65) µmol m⁻² s⁻¹ in FR-Fon (FR-Gri), lower than random uncertainty, 1.03 (0.73) µmol m⁻² s⁻¹ in FR-Fon (FR-Gri). Comparison between the DW-CS method using DW-NEE and NT method using EC<sub>s</sub>-NEE (Figure S 6, sites combined) yields higher bias, absolute daily error, and an increase in the underestimation. This is due to the already existing differences between DW- and EC<sub>s</sub>-NEE (see Figure 6).

<figure><img src="image_9.png" /><figcaption>Figure 8. (a) Daily averaged GPP and R<sub>eco</sub> calculated using standard night-time partitioning (NT, blue) and wavelet-based direct partitioning DW-CS (DW, orange), using DW-NEE as base data. Positive values show R<sub>eco</sub> and negative values show GPP. Dotted vertical lines show the start or end of the season (calendar years for forest site FR-Fon and cropping season for FR-Gri). (b) Daily</figcaption></figure>
NT versus DW GPP and R<sub>eco</sub> (both on the same graph), in grey 1:1, in orange linear fit. On the bottom statistics for GPP and R<sub>eco</sub> combined the correlation coefficient (R<sup>2</sup>), the mean error (ME, µmol m<sup>-2</sup> s<sup>-1</sup>), the mean absolute error (MAE, µmol m<sup>-2</sup> s<sup>-1</sup>) and the linear fit.

In the forest site, differences between DW-CS and NT during Summer 2019 (particularly warm) and 2021 (particularly rainy) fall into periods with peak temperatures (June and July 2019, August 2021) or precipitation (August 2019, June 2021).

In the crop site, DW-CS Reco estimations were higher on a few occasions than NT, when NT estimated an erroneously positive GPP. In October 2021, this happened after the intermediate crop was harvested and herbicides were used, possibly generating a pulse in Reco that was captured by the direct partitioning method. On 26/07/2021, the use of solid manure before barley seeding also generated erroneous positive GPP. In August 2021, a similar pulse was observed which happens after crop harvesting.

Unsurprisingly, given the methods equation, half-hourly results showed different R<sub>eco</sub> diel patterns between methods (Figure 9, DT included in Figure S 8). NT-R<sub>eco</sub> increased smoothly with temperature; DW-CS- R<sub>eco</sub> was flatter during the night and decreased during sunrise and sunset and showed an inverted U-shape curve during the day. Depending on the developmental stage, the day-time R<sub>eco</sub> can be larger (spring, peak season) or smaller (senescence, summer) than the night-time R<sub>eco</sub>.

<figure><img src="image_10.png" /><figcaption>Figure 9. Diel patterns of R<sub>eco</sub> and GPP estimated by wavelet-based conditional sampling (DWCS) and by standard night-time modelling (NT) during climatic seasons (FR-Fon) and the phenophases of green-up, peak season, senescence, and bare soil (excluding September 2021 due to intermediate crop), months are indicated by their first letter in parentheses. Note R<sub>eco</sub> and GPP are not in the same scale.</figcaption></figure>

A closer inspection on φ<sub>c</sub> components (Figure 10) reveals diel patterns resemblance between φ<sub>c</sub><sup>+</sup>|φ<sub>v</sub><sup>+</sup> and soil respiration (R<sub>soil</sub>) (daytime decrease during certain seasons), and φ<sub>c</sub><sup>+</sup>|φ<sub>v</sub><sup>-</sup> and plant respiration (R<sub>plant</sub>) (bimodal, with a maximum during
daytime) found in (Järveoja et al., 2020). R<sub>soil</sub> responds to soil temperature rather than air temperature and so follow the delayed warming and cooling pattern of the soil at the depth where respiration is maximum. We found here that φ<sub>c</sub><sup>+</sup>|φ<sub>v</sub><sup>+</sup> diel pattern follows quite well soil temperature around 20 cm depth (at 16 (30) cm depth for FR-Fon (FR-Gri) due to difference in available measurement depths) while φ<sub>c</sub><sup>+</sup>|φ<sub>v</sub><sup>-</sup> follows rather closely air temperature and incoming radiation.

490

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Diel patterns of φ<sub>c</sub><sup>+</sup>|φ<sub>v</sub><sup>+</sup> (associated to heterotrophic respiration, Rh) and φ<sub>c</sub><sup>+</sup>|φ<sub>v</sub><sup>-</sup> (associated to autotrophic respiration, Ra) and Rh's and Ra's main abiotic controls including air temperature (T<sub>air</sub>), soil temperature at 16 (30) cm depth in FR-Fon (FR-Gri) (ΔT<sub>soil</sub>, showed as deviation from the seasonal mean for readability), and photosynthetic photon flux density (PPFD) during climatic seasons (FR-Fon) and the phenophases of green-up, peak season, senescence, and bare soil, months indicated in parentheses. (cf. (Järveoja et al., 2020) Fig. 4)</figcaption></figure>

495

# 4. Discussion

## 4.1. On the differences between standard- and DW-EC

Results suggest that the DW method successfully captures the NEE dynamics observed by standard EC on a half-hourly (Figure 4) and annually (Figure 6) basis, highlighting its potential as a reliable alternative for flux analyses. The proposed DW-EC method obtained around 20% more high-quality data by not requiring stationarity. Naturally, reducing the gap fillings to zero is impossible since a certain amount of gap filling will always be necessary, even without technical problems and insufficient turbulence. Still, gap-filling is essential for retrieving continuous data series, but its use should be limited to the strict necessity. Even largely employed methods, such as Marginal Distribution Sampling (MDS), has shown poorer performance during night-time due to fewer observations (Moffat et al., 2007) and higher latitudes due to a skewed radiation distribution (Vekuri et al., 2023). Part of this is because standard EC method cannot handle non-stationary CO<sub>2</sub> flux, flagging

500

505
it as unreliable. Still, the increase in high-quality data can help improve the performance of gap-filling itself (Moffat et al., 2007), even if this benefit did not seen to be relevant for our two sites.

The DW underestimation most probably emerges from the detrending nature of wavelet transform. The decomposition
effectively disentangles each frequency, subtracting trends or fluctuations from the signal that span periods longer than the
frequency under consideration. In this case, low-frequency correction could help decrease the difference between EC and
DW, since in other contexts it has shown to reduce the difference between EC calculated using different detrending
strategies (including no detrending) to as low as 1 % (Mauder et al., 2021; Moncrieff et al., 2006; Rannik and Vesala, 1999).
Differently, including DW's lower frequencies (> 30 minutes⁻¹) is not as simple since although it requires to make
assumption on the cospectra in the low-frequency range as well as that lower frequencies use more 'neighbour' data making
continuity of good quality data more of an issue.

To enforce comparison between methods, using the same pre-processing corrections was important. However, some relevant
questions should be raised about this choice. Time lag and axis rotation corrections, in particular, transform w and CO₂,
creating artificial breaking points between observations of two neighbour half-hours. Both corrections were built for ECs,
where each half-hour is separated from the other. Wavelet decomposition, however, does use neighbour values at all times
and will assign these breaks to corresponding frequencies. We employed the commonly used double rotation method for axis
rotation, but planar fit, often recommended, should prevent these breaking points. Finnigan et al. (2003) have shown that 30-
min double-rotation is equivalent to high pass filtering, but may also add up part of the horizontal fluxes into the vertical
flux, thereby biasing the measured flux.

On another aspect, for sticky compounds such as ammonia (Ferrara et al., 2012) or VOCs (Loubet et al., 2022), the time lag
may be different for different concentration fluctuation frequencies. Using the wavelet decomposition to retrieve frequency
dependent time lags may be tested, especially for ammonia which have shown to lead to systematic underestimation of eddy
covariance fluxes, and show asymmetrical correlation functions (Ferrara et al., 2012 Fig 5). In sum, further wavelet-based
flux calculations may require revisiting some current pre-processing methods.

## 4.2. On the seasonal differences in the cospectra

NEE cospectra matched reasonably well with standard modelled cospectra (Horst, 1997) in peak and shape (Figure S 3). The forest site showed low-frequency attenuation, which is less visible in the crop site. Positive ($\phi_c^+$) and negative ($\phi_c^-$) parts showed seasonal patterns not seen in the theoretical curve, suggesting it shall be explained by something different than micrometeorological factors (wind speed or Obukov length) and measuring height (for the agricultural site) obtained from the modelled cospectra (Horst, 1997). Indeed, $\phi_c^-$ is expected to come mainly from the leaves, while $\phi_c^+$ is from the ground and leaves and brought up from in-canopy by injection. Intuitively, the second process moves eddies forward on the Kolmogorov cascade by encountering obstacles between ground and canopy, while the first coming from the top of the canopy is less influenced by obstacles. This difference in eddy size transporting the signal would explain why, during photosynthetically active months, $\phi_c^-$ shows lower frequencies than $\phi_c^+$ during neutral and unstable conditions (mostly
daytime). The seasonality would thus come from the absence (presence) of leaves and crops, which leads to a reduced (increased) number of “obstacles”. In stable conditions (mostly night), the same pattern is not seen because $\phi_c$ is expected to be small, and indeed, at times, it differs significantly from the characteristic NEE cospectra.

Of course, cospectra analysis is not a specificity of DW nor any frequency decomposed method; Fourier transform may be used in post-processing spectral corrections. Wavelets remain however advantageous by not requiring stationarity. Wind-velocity coordinate rotations may however be performed over longer periods than 30 min. Planar fit approaches should therefore be preferred.

**4.3. On partitioning methods and possible sources of error**

In standard practice, a modelled response of NEE to light, water, and temperature over days is used to split it into GPP and $R_{eco}$. The standard night-time method (NT) employs a nocturnal temperature response model to estimate $R_{eco}$ (Reichstein et al., 2005), while day-time (DT) incorporates a light response curve to estimate both GPP and $R_{eco}$ during day-time (Keenan et al., 2019; Lasslop et al., 2010). We could expect little difference when comparing the same partitioning method on standard EC and wavelet-based fluxes, given their similarity (Figure 6).

However, a particularly interesting result can be seen in August 2021 in FR-Gri, when the increase in high-quality data in DW-NEE made NT's estimation undoubtedly more reasonable compared to its prediction using EC<sub>s</sub>-NEE (Figure 7). On this occasion, NT relying on EC<sub>s</sub> projected an unrealistic rise in absolute GPP following intermediate crop harvest and herbicide application. Conversely, NT based on DW-EC showed the expected absolute GPP decrease, albeit with implausible positive GPP values (a known issue for the NT method). Nonetheless, standard model-based partitioning yielded somewhat different estimations at times (Figure 7), despite being informed by relatively similar input data (Figure 6), underscoring model-based uncertainty. In addition to that, by rendering partitioning contingent on distinct periods, spatial heterogeneity becomes an issue (Wehr et al., 2016). Direct partitioning methods are based on single 30-minute periods and are, therefore, relatively free from this issue, given that neighbourhood influence is limited.

Measuring directly GPP and $R_{eco}$ at the ecosystem level poses challenges, yielding inconclusive comparative studies. A more direct way of measuring Reco is by using dark chambers (Järveoja et al., 2020) or using carbon isotopes (more precisely the ratio between <sup>13</sup>C to <sup>12</sup>C) (Wehr et al., 2016). However, $R_{eco}$ by NT partitioning was found 25% higher than isotopic-derived $R_{eco}$ fluxes in a deciduous temperate forest during June-July (Wehr et al., 2016) and 16% to 22% higher than automatic dark chambers observations in a peatland (Järveoja et al., 2020). This is often attributed to models' limited capacity to replicate diel patterns (Wehr et al., 2016; Keenan et al., 2019; Griffis et al., 2004; Järveoja et al., 2020). In particular the dynamics of root-microbe-soil system is not well characterised. Some large-scale girdling experiments however show that soil respiration is highly correlated to photosynthesis supply to roots: respiration was found to be reduced by 37% within 5 days (54% in 1-2 months), after stopping the supply of photosynthesis to roots (Högberg et al., 2001). <sup>13</sup>C labelling studies further showed a 35 hours half-life of soil respiratory efflux in a forest (Högberg et al. 2008). NT's nocturnal calibration on respiration does not account for photosynthesis transfer processes. In addition to that, the inhibition of leaf respiration by light is estimated to
cause systematic overestimation of daytime ecosystem basal respiration estimated using NT (Wehr et al., 2016; Keenan et al., 2019). In sum, these references indicate that respiration dial pattern may be much more complex than what is currently include in temperature response models.

In this study, we propose a new direct partitioning method, the discrete wavelet-based conditional sampling (DW-CS), based on the conditional sampling of frequency decomposed w and χ<sub>CO₂</sub> product (Eq. (18)). Overall, daily mean NT estimations were 12 (22) % higher than DW-CS for GPP (R<sub>eco</sub>). Much noticeably, the DW-CS R<sub>eco</sub> diel cycle unveils a bimodal trend (Figure 9) which as previously observed by thorough chamber methods measurements (Järveoja et al., 2020) and on the ¹³C/¹²C derived estimations (Wehr et al., 2016). This diurnal variation was attributed in Wehr et al. (2016) to the inhibition of leaf respiration in light, and in Järveoja et al. (2020) to the differential response of R<sub>eco</sub> to soil temperature and air or plant temperature. Another noticeable result show in Figure 9, is that during the crop's peak growing season, DW-CS R<sub>eco</sub> is much higher than NT R<sub>eco</sub> during the day. This increase may be explained by an increase in autotrophic respiration concomitant with higher GPP from plant growth during that period. This feature cannot obviously be captured by the NT method which uses nocturnal calibration.

In Figure 10, the similarities between φ<sub>c</sub><sup>+</sup>|φ<sub>v</sub><sup>+</sup> (φ<sub>c</sub><sup>+</sup>|φ<sub>v</sub><sup>-</sup>) and R<sub>soil</sub> (R<sub>plant</sub>) from Järveoja et al., (2020) indicate the potential to perform even further detailed partitioning. The similarities may relate to soil evaporation being higher than plant evaporation when a respiration signal (φ<sub>c</sub><sup>+</sup>) is measured. Indeed, R<sub>plant</sub> and GPP are both dependent on stomata and so a cut in GPP, for instance due to lower incoming radiation induced by cloud or shadowing, could cause either or both stomatal closure and a decrease in plant surface temperature. Whether by physical constraint or condensation, during these moment R<sub>plant</sub> would come with a negative water flux (φ<sub>c</sub><sup>+</sup>|φ<sub>v</sub><sup>-</sup>). Evaluating the proposed partitioning method would require measuring R<sub>eco</sub>, R<sub>soil</sub> and R<sub>plant</sub> at the field scale, which require further research.

## 4.4 Perspectives on using wavelet-based EC for less gap-filling and direct partitioning

In this study, we have explored how including non-stationary fluxes, which are omitted by standard EC method, modified the computed CO₂ flux and further gapfilling. These periods carry real information on the surface flux. For instance, dynamic light environments can trigger rapid but non-coordinated photosynthesis and stomatal response (McAusland et al., 2016), possibly leading to non-stationary NEE. When filtered out, those non-stationary events effectively “blinds” the gap-filling methods and final users from these transition events.

Similarly, non-homogeneous footprints are often encountered in ecosystem monitoring sites, although everything is done to minimise these conditions. Agricultural fields and sub-urban and urban areas are especially prone to source heterogeneity, mainly in the shape of a local intensive anthropogenic source (animal grazing on the field, animal barns, tractors, roads, chimneys) that may also lead to changes in CO₂ concentrations and fluxes with wind direction and hence non-stationarity in these components (Crawford and Christen, 2015). Thus, even if the impact is arguably small on monthly and annual net flux budgets on ecosystem towers, using DW-EC becomes especially relevant in setups or situations with multiple local sources that are hard to isolate.
Intermittent turbulence was identified as the main problem for nocturnal EC, which is leads to the u\* filtering approach (Aubinet, 2008). It is important to note that in certain cases of such intermittent turbulence, a non-stationary flux may be delayed from the process that generated this flux. Indeed, when very low turbulence is followed by a burst of wind (ejections or sweeps), measured flux includes releasing accumulated stock (Katul et al., 2006; Kaimal and Finnigan, 1994). As a consequence, including non-stationary fluxes retrieved by DW in the standard model-based partitioning methods, which do not consider these peaks as night-time delayed respiration, may lead to biased gapfilling. This bias, however, should affect less direct partitioning using DW-CS, which exclusively relies on data from the same snapshot of time.

The conditional sampling method presented here could be further developed to use soil and plant, heterotrophic and autotrophic or even biogenic and anthropogenic tracers like carbon and water isotopes to improve our understanding of the carbon cycling in an ecosystem. Including co-produced gases can be the key to perform more elaborated attribution of fluxes to ecosystem compartments. Carbonyl sulphide (COS) is a known tracer of photosynthesis (Maseyk et al., 2014). Combining CO₂ flux with COS flux together with H₂O would lead to a way to further partition GPP and R eco. In another context, carbon monoxide has been used to identify fossil fuel emissions (Super et al., 2017). Wavelet-based conditional sampling emerges as a promising framework for integrating such data, either directly as proposed here or through hybrid methods in which more elaborated models could be used to refine the partitioning method.

**5. Conclusions**

Discrete Wavelet-based Eddy Covariance (DW-EC) yielded around 17 % and 29 % fewer gaps than standard Eddy Covariance (ECs) over four years of data in a forest site (FR-Fon) and two years in a crop site (FR-Gri) in the French Parisian region, respectively. We can expect even larger gap differences in perturbated environments (topography, inhomogeneous areas). The half-hour high-quality NEE (stationary and within well-developed turbulence) computed by wavelets were highly correlated to standard eddy covariance (R²= 0.98 for both FR-Fon and FR-Gri), worsening for medium (0.73 and 0.52) and low-quality (0.03 and 0.0) data. At the daily scale, this correlation was kept (R²= 0.97), but with a slight bias with DW around 9% lower than ECs (mean error = 0.1 and 0.18 µmol m s⁻¹, mean absolute error = 0.42 and 0.39 µmol m⁻² s⁻¹, ECs=1.09×DW by linear fitting for FR-Fon and FR-Gri). This effect is likely related to the detrending nature of wavelet decomposition which leads to low-frequency attenuation of the flux. The supposed advantage of reduced gaps for DW-EC lead to sensibly similar NEE budgets: +2 (−1) % in FR-Fon (FR-Gri) when compared to DW-EC forced to have the same gaps as ECs. This suggests that for standard sites (mostly homogeneous and flat) moving towards DW-EC would not significantly improve the annual budget. However, partitioning using ECs- and DW-NEE yielded different GPP and R eco, particularly for the night-time method (NT), where more high-quality observations made estimations arguably more credible in the crop site.

A new partitioning method is proposed, combining discrete wavelet transform and conditional sampling (DW-CS). The method splits positive and negative parts of the product of the wavelet decomposed vertical component of the wind, w\*(j),
and a scalar, $\chi_s'(j)$. The underlying empirical assumption is that wavelet decomposition should allow to trap in each frequency the positive and negative “gusts” which are mixed up in the original signal. Further including PPFD and $w' \chi_{H_2O}'$, to attribute unrealistic CO₂ fluxes, led to a method for estimating $R_{eco}$ and GPP. Compared to DW-CS, Night-time partitioning (NT) showed better correlation and smaller errors than day-time methods. Mean absolute errors between NT and DW-CS (0.8 and 0.65 µmol m⁻² s⁻¹ in FR-Fon and FR-Gri) were lower than the NEE random uncertainty (1.03 and 0.73 µmol m⁻² s⁻¹ in FR-Fon and FR-Gri). But most noticeably, DW-CS led to a different $R_{eco}$ diel pattern compared to temperature only driven models, with a daily respiration pattern that follows radiation (and hence GPP) and a night-time pattern that follows soil respiration. This diel pattern was already observed using chambers and has some ground to be more realistic than the standard NT and DT approaches: this pattern may reflect either a differentiated temperature response from soil and plants, a light inhibition response from plants, or a time shift between photosynthesates production and their transport to roots. Our DW-CS approach was not validated by field measurements of net flux components to confirm that this respiration pattern was really happening in the observed sites. This study however strongly suggests to further evaluate $R_{eco}$ diel pattern as it may have strong impacts on how global CO₂ cycle is modelled. The DW-CS we present here should be further tested and refined as it has the benefit of integrating at the field scale without needing extra measurements which also allows reprocessing of old data. We also note that DW-CS could be developed to incorporate other tracers like COS to better partition the CO₂ fluxes between ecosystem compartments.

Eddy covariance has improved observations and, indirectly, models for the last decades. This study shows that standardising wavelets for EC measurements can be operational using discrete wavelets decomposition. This would be very beneficial as it includes non-stationary data, and hence reduce gaps, and allows a look into transitory process. The simplicity and flexibility of DW-EC make it also very easy for (re)analysis. The method would be powerful for CH₄ and N₂O fluxes which are highly non-stationary and mode difficult to gap-fill, as well as in urban setups for the same reasons. The new direct partitioning method shows great promises in providing fully-observation-based partitioning at the field scale. However, partitioning methods, in general, and our new wavelet-based method in particular, need further validation experiments across ecosystems and environmental conditions.

## Acknowledgements

We acknowledge the ICOS ERIC for providing the FR-Gri and FR-Fon data. This project has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No 101037319 (PAUL ICOS-Cities project).
## Data availability

ICOS data for the FR-Gri site can be downloaded from the carbon portal: https://meta.icos-cp.eu/resources/stations/ES_FR-Gri.

ICOS data for the FR-Fon site can be downloaded from the carbon portal: https://meta.icos-cp.eu/resources/stations/ES_FR-Fon.

## Code

Code used in the analysis presented in this paper is available online and can be accessed at https://github.com/pedrohenriquecoimbra/coimbra-et-al-wavelet-based-partitioning

## References

Arts, L. P. A. and van den Broek, E. L.: The fast continuous wavelet transformation (fCWT) for real-time, high-quality, noise-resistant time-frequency analysis, Nat. Comput. Sci., 2, 47–58, https://doi.org/10.1038/s43588-021-00183-z, 2022.

Aubinet, M.: Eddy Covariance CO2 Flux Measurements in Nocturnal Conditions: An Analysis of the Problem, Ecol. Appl., 18, 1368–1378, https://doi.org/10.1890/06-1336.1, 2008.

Aubinet, M., Grelle, A., Ibrom, A., Rannik, Ü., Moncrieff, J., Foken, T., Kowalski, A. S., Martin, P. H., Berbigier, P., Bernhofer, Ch., Clement, R., Elbers, J., Granier, A., Grünwald, T., Morgenstern, K., Pilegaard, K., Rebmann, C., Snijders, W., Valentini, R., and Vesala, T.: Estimates of the Annual Net Carbon and Water Exchange of Forests: The EUROFLUX Methodology, in: Advances in Ecological Research, vol. 30, edited by: Fitter, A. H. and Raffaelli, D. G., Academic Press, 113–175, https://doi.org/10.1016/S0065-2504(08)60018-5, 1999.

Aubinet, M., Heinesch, B., and Longdoz, B.: Estimation of the carbon sequestration by a heterogeneous forest: night flux corrections, heterogeneity of the site and inter-annual variability, Glob. Change Biol., 8, 1053–1071, https://doi.org/10.1046/j.1365-2486.2002.00529.x, 2002.

Aubinet, M., Berbigier, P., Bernhofer, Ch., Cescatti, A., Feigenwinter, C., Granier, A., Grünwald, Th., Havrankova, K., Heinesch, B., Longdoz, B., Marcolla, B., Montagnani, L., and Sedlak, P.: Comparing CO2 Storage and Advection Conditions at Night at Different Carboeuroflux Sites, Bound.-Layer Meteorol., 116, 63–93, https://doi.org/10.1007/s10546-004-7091-8, 2005.

Aubinet, M., Feigenwinter, C., Heinesch, B., Laffineur, Q., Papale, D., Reichstein, M., Rinne, J., and Van Gorsel, E.: Nighttime Flux Correction, in: Eddy Covariance: A Practical Guide to Measurement and Data Analysis, edited by: Aubinet, M., Vesala, T., and Papale, D., Springer Netherlands, Dordrecht, 133–157, https://doi.org/10.1007/978-94-007-2351-1_5, 2012.

Baldocchi, D. D.: Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: past, present and future, Glob. Change Biol., 9, 479–492, https://doi.org/10.1046/j.1365-2486.2003.00629.x, 2003.

Canadell, J. G., Monteiro, P. M. S., Costa, M. H., Cotrim da Cunha, L., Cox, P. M., Eliseev, A. V., Henson, S., Ishii, M., Jaccard, S., Koven, C., Lohila, A., Patra, P. K., Piao, S., Rogelj, J., Syampungani, S., Zaehle, S., and Zickfeld, K.: Global
700

Carbon and other Biogeochemical Cycles and Feedbacks, edited by: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Clim. Change 2021 Phys. Sci. Basis Contrib. Work. Group Sixth Assess. Rep. Intergov. Panel Clim. Change, 673–816, https://doi.org/10.1017/9781009157896.007, 2021.

705

Crawford, B. and Christen, A.: Spatial source attribution of measured urban eddy covariance CO2 fluxes, Theor Appl Climatol., 119, 733–755, https://doi.org/10.1007/s00704-014-1124-0, 2015.

Daubechies, I.: Orthonormal bases of compactly supported wavelets, Commun. Pure Appl. Math., 41, 909–996, https://doi.org/10.1002/cpa.3160410705, 1988.

710

Delpierre, N., Berveiller, D., Granda, E., and Dufrêne, E.: Wood phenology, not carbon input, controls the interannual variability of wood growth in a temperate oak forest, New Phytol., 210, 459–470, https://doi.org/10.1111/nph.13771, 2016.

Desjardins, R. L., Worth, D. E., Pattey, E., VanderZaag, A., Srinivasan, R., Mauder, M., Worthy, D., Sweeney, C., and Metzger, S.: The challenge of reconciling bottom-up agricultural methane emissions inventories with top-down measurements, Agric. For. Meteorol., 248, 48–59, https://doi.org/10.1016/j.agrformet.2017.09.003, 2018.

Du, Q., Liu, H., Feng, J., and Wang, L.: Effects of different gap filling methods and land surface energy balance closure on annual net ecosystem exchange in a semiarid area of China, Sci. China Earth Sci., 57, 1340–1351, https://doi.org/10.1007/s11430-013-4756-5, 2014.

715

Duffy, K. A., Schwalm, C. R., Arcus, V. L., Koch, G. W., Liang, L. L., and Schipper, L. A.: How close are we to the temperature tipping point of the terrestrial biosphere?, Sci. Adv., 7, eaay1052, https://doi.org/10.1126/sciadv.aay1052, 2021.

Falge, E., Baldocchi, D., Olson, R., Anthoni, P., Aubinet, M., Bernhofer, C., Burba, G., Ceulemans, R., Clement, R., Dolman, H., Granier, A., Gross, P., Grünwald, T., Hollinger, D., Jensen, N.-O., Katul, G., Keronen, P., Kowalski, A., Lai, C. T., Law, B. E., Meyers, T., Moncrieff, J., Moors, E., Munger, J. W., Pilegaard, K., Rannik, Ü., Rebmann, C., Suyker, A., Tenhunen, J., Tu, K., Verma, S., Vesala, T., Wilson, K., and Wofsy, S.: Gap filling strategies for defensible annual sums of net ecosystem exchange, Agric. For. Meteorol., 107, 43–69, https://doi.org/10.1016/S0168-1923(00)00225-2, 2001.

720

Farge, M.: Wavelet transforms and their applications to turbulence, Annu. Rev. Fluid Mech., 24, 395–458, 1992.

Farge, M. and Schneider, K.: Analysing and Computing Turbulent Flows Using Wavelets, in: New trends in turbulence Turbulence: nouveaux aspects: 31 July – 1 September 2000, edited by: Lesieur, M., Yaglom, A., and David, F., Springer, Berlin, Heidelberg, 449–504, https://doi.org/10.1007/3-540-45674-0_9, 2001.

725

Ferrara, R. M., Loubet, B., Di Tommasi, P., Bertolini, T., Magliulo, V., Cellier, P., Eugster, W., and Rana, G.: Eddy covariance measurement of ammonia fluxes: Comparison of high frequency correction methodologies, Agric. For. Meteorol., 158–159, 30–42, https://doi.org/10.1016/j.agrformet.2012.02.001, 2012.

Finnigan, J. J., Clement, R., Malhi, Y., Leuning, R., and Cleugh, H. A.: A Re-Evaluation of Long-Term Flux Measurement Techniques Part I: Averaging and Coordinate Rotation, Bound.-Layer Meteorol., 107, 1–48, https://doi.org/10.1023/A:1021554900225, 2003.

730

Foken, T., Aubinet, M., and Leuning, R.: The Eddy Covariance Method, in: Eddy Covariance: A Practical Guide to Measurement and Data Analysis, edited by: Aubinet, M., Vesala, T., and Papale, D., Springer Netherlands, Dordrecht, 1–19, https://doi.org/10.1007/978-94-007-2351-1_1, 2012.
735 Foken, Th. and Wichura, B.: Tools for quality assessment of surface-based flux measurements, Agric. For. Meteorol., 78, 83–105, https://doi.org/10.1016/0168-1923(95)02248-1, 1996.

Fowler, D., Hargreaves, K., Skiba, U., Milne, R., Zahniser, M., Moncrieff, J., Beverland, I., and Gallagher, M.: Measurements of Ch4 and N2o Fluxes at the Landscape Scale Using Micrometeorological Methods, Philos. Trans. R. Soc. - Math. Phys. Eng. Sci., 351, 339–355, https://doi.org/10.1098/rsta.1995.0038, 1995.

740 Göckede, M., Kittler, F., and Schaller, C.: Quantifying the impact of emission outbursts and non-stationary flow on eddy-covariance CH4 flux measurements using wavelet techniques, Biogeosciences, 16, 3113–3131, https://doi.org/10.5194/bg-16-3113-2019, 2019.

745 van Gorsel, E., Delpierre, N., Leuning, R., Black, A., Munger, J. W., Wofsy, S., Aubinet, M., Feigenwinter, C., Beringer, J., Bonal, D., Chen, B., Chen, J., Clement, R., Davis, K. J., Desai, A. R., Dragoni, D., Etzold, S., Grünwald, T., Gu, L., Heinesch, B., Hutyra, L. R., Jans, W. W. P., Kutsch, W., Law, B. E., Leclerc, M. Y., Mammarella, I., Montagnani, L., Noormets, A., Rebmann, C., and Wharton, S.: Estimating nocturnal ecosystem respiration from the vertical turbulent flux and change in storage of CO2, Agric. For. Meteorol., 149, 1919–1930, https://doi.org/10.1016/j.agrformet.2009.06.020, 2009.

750 Goulden, M. L., Munger, J. W., Fan, S.-M., Daube, B. C., and Wofsy, S. C.: Measurements of carbon sequestration by long-term eddy covariance: methods and a critical evaluation of accuracy, Glob. Change Biol., 2, 169–182, https://doi.org/10.1111/j.1365-2486.1996.tb00070.x, 1996.

Griffis, T. J., Black, T. A., Gaumont-Guay, D., Drewitt, G. B., Nesic, Z., Barr, A. G., Morgenstern, K., and Kljun, N.: Seasonal variation and partitioning of ecosystem respiration in a southern boreal aspen forest, Agric. For. Meteorol., 125, 207–223, https://doi.org/10.1016/j.agrformet.2004.04.006, 2004.

755 Gu, L., Falge, E. M., Boden, T., Baldocchi, D. D., Black, T. A., Saleska, S. R., Suni, T., Verma, S. B., Vesala, T., Wofsy, S. C., and Xu, L.: Objective threshold determination for nighttime eddy flux filtering, Agric. For. Meteorol., 128, 179–197, https://doi.org/10.1016/j.agrformet.2004.11.006, 2005.

760 Gu, L., Massman, W. J., Leuning, R., Pallardy, S. G., Meyers, T., Hanson, P. J., Riggs, J. S., Hosman, K. P., and Yang, B.: The fundamental equation of eddy covariance and its application in flux measurements, Agric. For. Meteorol., 152, 135–148, https://doi.org/10.1016/j.agrformet.2011.09.014, 2012.

Heskel, M. A., Atkin, O. K., Turnbull, M. H., and Griffin, K. L.: Bringing the Kok effect to light: A review on the integration of daytime respiration and net ecosystem exchange, Ecosphere, 4, art98, https://doi.org/10.1890/ES13-00120.1, 2013.

765 Högberg, P., Högberg, M. N., Göttlicher, S. G., Betson, N. R., Keel, S. G., Metcalfe, D. B., Campbell, C., Schindlbacher, A., Hurry, V., Lundmark, T., Linder, S., and Näsholm, T.: High temporal resolution tracing of photosynthate carbon from the tree canopy to forest soil microorganisms, New Phytologist, 177, 220–228, https://doi.org/10.1111/j.1469-8137.2007.02238.x, 2008.

770 Högberg, P., Nordgren, A., Buchmann, N., Taylor, A. F. S., Ekblad, A., Högberg, M. N., Nyberg, G., Ottosson-Löfvenius, M., and Read, D. J.: Large-scale forest girdling shows that current photosynthesis drives soil respiration, Nature, 411, 789–792, https://doi.org/10.1038/35081058, 2001.

Horst, T. W.: A simple formula for attenuation of eddy fluxes measured with first-order-response scalar sensors, Bound.-Layer Meteorol., 82, 219–233, https://doi.org/10.1023/A:1000229130034, 1997.
IPCC: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, In Press, https://doi.org/10.1017/9781009157896, 2021.

IPCC: Climate Change 2022: Mitigation of Climate Change, https://doi.org/10.1017/9781009157926, 2022a.

IPCC: Climate Change 2022: Impacts, Adaptation and Vulnerability, https://doi.org/10.1017/9781009325844, 2022b.

Irvin, J., Zhou, S., McNicol, G., Lu, F., Liu, V., Fluet-Chouinard, E., Ouyang, Z., Knox, S. H., Lucas-Moffat, A., Trotta, C., Papale, D., Vitale, D., Mammarella, I., Alekseychik, P., Aurela, M., Avati, A., Baldocchi, D., Bansal, S., Bohrer, G., Campbell, D. I., Chen, J., Chu, H., Dalmagro, H. J., Delwiche, K. B., Desai, A. R., Euskirchen, E., Feron, S., Goeckede, M., Heimann, M., Helbig, M., Helfter, C., Hemes, K. S., Hirano, T., Iwata, H., Jurasinski, G., Kalhori, A., Kondrich, A., Lai, D. Y., Lohila, A., Malhotra, A., Merbold, L., Mitra, B., Ng, A., Nilsson, M. B., Noormets, A., Peichl, M., Rey-Sanchez, A. C., Richardson, A. D., Runkle, B. R., Schäfer, K. V., Sonnentag, O., Stuart-Haëntjens, E., Sturtevant, C., Ueyama, M., Valach, A. C., Vargas, R., Vourlitis, G. L., Ward, E. J., Wong, G. X., Zona, D., Alberto, Ma. C. R., Billesbach, D. P., Celis, G., Dolman, H., Friborg, T., Fuchs, K., Gogo, S., Gondwe, M. J., Goodrich, J. P., Gottschalk, P., Hörtnagl, L., Jacotot, A., Koebsch, F., Kasak, K., Maier, R., Morin, T. H., Nemitz, E., Oechel, W. C., Oikawa, P. Y., Ono, K., Sachs, T., Sakabe, A., Schuur, E. A., Shortt, R., Sullivan, R. C., Szutu, D. J., Tuittila, E.-S., Varlagin, A., Verfaillie, J. G., Wille, C., Windham-Myers, L., Poulter, B., and Jackson, R. B.: Gap-filling eddy covariance methane fluxes: Comparison of machine learning model predictions and uncertainties at FLUXNET-CH4 wetlands, Agric. For. Meteorol., 308–309, 108528, https://doi.org/10.1016/j.agrformet.2021.108528, 2021.

Järveoja, J., Nilsson, M. B., Crill, P. M., and Peichl, M.: Bimodal diel pattern in peatland ecosystem respiration rebuts uniform temperature response, Nat. Commun., 11, 4255, https://doi.org/10.1038/s41467-020-18027-1, 2020.

Jia, G., Shevliakova, E., Artaxo, P., De Noblet-Ducoudré, N., Houghton, R., House, J., Kitajima, K., Lennard, C., Popp, A., and Sirin, A.: Land-climate interactions, in: Climate Change and Land: an IPCC special report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems, 2019.

Kaimal, J. C. and Finnigan, J. J.: Atmospheric Boundary Layer Flows: Their Structure and Measurement, Oxford University Press, 304 pp., 1994.

Kaimal, J. C., Wyngaard, J. C., Izumi, Y., and Coté, O. R.: Spectral characteristics of surface-layer turbulence, Q. J. R. Meteorol. Soc., 98, 563–589, https://doi.org/10.1002/qj.49709841707, 1972.

Katul, G., Porporato, A., Cava, D., and Siqueira, M.: An analysis of intermittency, scaling, and surface renewal in atmospheric surface layer turbulence, Physica D: Nonlinear Phenomena, 215, 117–126, https://doi.org/10.1016/j.physd.2006.02.004, 2006.

Keenan, T. F., Migliavacca, M., Papale, D., Baldocchi, D., Reichstein, M., Torn, M., and Wutzler, T.: Widespread inhibition of daytime ecosystem respiration, Nat. Ecol. Evol., 3, 407–415, https://doi.org/10.1038/s41559-019-0809-2, 2019.

Klosterhalfen, A., Graf, A., Brüggemann, N., Drüe, C., Esser, O., González-Dugo, M. P., Heinemann, G., Jacobs, C. M. J., Mauder, M., Moene, A. F., Ney, P., Pütz, T., Rebmann, C., Ramos Rodríguez, M., Scanlon, T. M., Schmidt, M., Steinbrecher, R., Thomas, C. K., Valler, V., Zeeman, M. J., and Vereecken, H.: Source partitioning of H₂O and CO₂ fluxes based on high-frequency eddy covariance data: a comparison between study sites, Biogeosciences, 16, 1111–1132, https://doi.org/10.5194/bg-16-1111-2019, 2019.
810

Kowalski, A. S. and Serrano-Ortiz, P.: On the relationship between the eddy covariance, the turbulent flux, and surface exchange for a trace gas such as CO2, Bound.-Layer Meteorol., 124, 129–141, https://doi.org/10.1007/s10546-007-9171-z, 2007.

815

Lasslop, G., Reichstein, M., Papale, D., Richardson, A. D., Arneth, A., Barr, A., Stoy, P., and Wohlfahrt, G.: Separation of net ecosystem exchange into assimilation and respiration using a light response curve approach: critical issues and global evaluation, Glob. Change Biol., 16, 187–208, https://doi.org/10.1111/j.1365-2486.2009.02041.x, 2010.

Lee, G. R., Gommers, R., Waselewski, F., Wohlfahrt, K., and O’Leary, A.: PyWavelets: A Python package for wavelet analysis, J. Open Source Softw., 4, 1237, https://doi.org/10.21105/joss.01237, 2019.

Lloyd, J. and Taylor, J. A.: On the Temperature Dependence of Soil Respiration, Funct. Ecol., 8, 315–323, https://doi.org/10.2307/2389824, 1994.

820

Loubet, B., Laville, P., Lehuger, S., Larmanou, E., Fléchard, C., Mascher, N., Genermont, S., Roche, R., Ferrara, R. M., Stella, P., Personne, E., Durand, B., Decuq, C., Flura, D., Masson, S., Fanucci, O., Rampon, J.-N., Siemens, J., Kindler, R., Gabrielle, B., Schrumpf, M., and Cellier, P.: Carbon, nitrogen and Greenhouse gases budgets over a four years crop rotation in northern France, Plant Soil, 343, 109, https://doi.org/10.1007/s11104-011-0751-9, 2011.

825

Loubet, B., Buysse, P., Gonzaga-Gomez, L., Lafouge, F., Ciuraru, R., Decuq, C., Kammer, J., Bsaibes, S., Boissard, C., Durand, B., Gueudet, J.-C., Fanucci, O., Zurfluh, O., Abis, L., Zannoni, N., Truong, F., Baisnée, D., Sarda-Estève, R., Staudt, M., and Gros, V.: Volatile organic compound fluxes over a winter wheat field by PTR-Qi-TOF-MS and eddy covariance, Atmospheric Chem. Phys., 22, 2817–2842, https://doi.org/10.5194/acp-22-2817-2022, 2022.

Lucas-Moffat, A. M., Schrader, F., Herbst, M., and Brümmer, C.: Multiple gap-filling for eddy covariance datasets, Agric. For. Meteorol., 325, 109114, https://doi.org/10.1016/j.agrformet.2022.109114, 2022.

830

Mallat, S. G.: A theory for multiresolution signal decomposition: the wavelet representation, IEEE Trans. Pattern Anal. Mach. Intell., 11, 674–693, https://doi.org/10.1109/34.192463, 1989.

Maseyk, K., Berry, J. A., Billesbach, D., Campbell, J. E., Torn, M. S., Zahniser, M., and Seibt, U.: Sources and sinks of carbonyl sulfide in an agricultural field in the Southern Great Plains, Proc. Natl. Acad. Sci., 111, 9064–9069, https://doi.org/10.1073/pnas.1319132111, 2014.

835

Massman, W. J. and Lee, X.: Eddy covariance flux corrections and uncertainties in long-term studies of carbon and energy exchanges, Agric. For. Meteorol., 113, 121–144, https://doi.org/10.1016/S0168-1923(02)00105-3, 2002.

Mauder, M. and Foken, T.: Documentation and Instruction Manual of the Eddy-Covariance Software Package TK3, 60, 2011.

840

Mauder, M., Foken, T., Aubinet, M., and Ibrom, A.: Eddy-Covariance Measurements, in: Springer Handbook of Atmospheric Measurements, edited by: Foken, T., Springer International Publishing, Cham, 1473–1504, https://doi.org/10.1007/978-3-030-52171-4_55, 2021.

Mauder, M., Desjardins, R. L., and MacPherson, I.: Scale analysis of airborne flux measurements over heterogeneous terrain in a boreal ecosystem, J. Geophys. Res. Atmospheres, 112, https://doi.org/10.1029/2006JD008133, 2007.

845

Mauder, M., Cuntz, M., Drüe, C., Graf, A., Rebmann, C., Schmid, H. P., Schmidt, M., and Steinbrecher, R.: A strategy for quality and uncertainty assessment of long-term eddy-covariance measurements, Agric. For. Meteorol., 169, 122–135, https://doi.org/10.1016/j.agrformet.2012.09.006, 2013.
Metzger, S.: Surface-atmosphere exchange in a box: Making the control volume a suitable representation for in-situ observations, Agric. For. Meteorol., 255, 68–80, https://doi.org/10.1016/j.agrformet.2017.08.037, 2018.

Metzger, S., Junkermann, W., Mauder, M., Butterbach-Bahl, K., Trancón y Widemann, B., Neidl, F., Schäfer, K., Wieneke, S., Zheng, X. H., Schmid, H. P., and Foken, T.: Spatially explicit regionalization of airborne flux measurements using environmental response functions, Biogeosciences, 10, 2193–2217, https://doi.org/10.5194/bg-10-2193-2013, 2013.

Mishurov, M. and Kiely, G.: Gap-filling techniques for the annual sums of nitrous oxide fluxes, Agric. For. Meteorol., 151, 1763–1767, https://doi.org/10.1016/j.agrformet.2011.07.014, 2011.

Moffat, A. M.: A new methodology to interpret high resolution measurements of net carbon fluxes between terrestrial ecosystems and the atmosphere, 2012.

Moffat, A. M., Papale, D., Reichstein, M., Hollinger, D. Y., Richardson, A. D., Barr, A. G., Beckstein, C., Braswell, B. H., Churkina, G., Desai, A. R., Falge, E., Gove, J. H., Heimann, M., Hui, D., Jarvis, A. J., Kattge, J., Noormets, A., and Stauch, V. J.: Comprehensive comparison of gap-filling techniques for eddy covariance net carbon fluxes, Agric. For. Meteorol., 147, 209–232, https://doi.org/10.1016/j.agrformet.2007.08.011, 2007.

Moncrieff, J., Clement, R., Finnigan, J., and Meyers, T.: Averaging, Detrending, and Filtering of Eddy Covariance Time Series, in: Handbook of Micrometeorology, Atmospheric and Oceanographic Sciences Library, vol. 29, 7–31, https://doi.org/10.1007/1-4020-2265-4_2, 2006.

Moncrieff, J. B., Massheder, J. M., de Bruin, H., Elbers, J., Friborg, T., Heusinkveld, B., Kabat, P., Scott, S., Soegaard, H., and Verhoef, A.: A system to measure surface fluxes of momentum, sensible heat, water vapour and carbon dioxide, J. Hydrol., 188–189, 589–611, https://doi.org/10.1016/S0022-1694(96)03194-0, 1997.

Papale, D. and Valentini, R.: A new assessment of European forests carbon exchanges by eddy fluxes and artificial neural network spatialization, Glob. Change Biol., 9, 525–535, https://doi.org/10.1046/j.1365-2486.2003.00609.x, 2003.

Papale, D., Reichstein, M., Aubinet, M., Canfora, E., Bernhofer, C., Kutsch, W., Longdoz, B., Rambal, S., Valentini, R., Vesala, T., and Yakir, D.: Towards a standardized processing of Net Ecosystem Exchange measured with eddy covariance technique: algorithms and uncertainty estimation, Biogeosciences, 3, 571–583, https://doi.org/10.5194/bg-3-571-2006, 2006.

Pastorello, G., Trotta, C., Canfora, E., Chu, H., Christianson, D., Cheah, Y.-W., Poindexter, C., Chen, J., Elbashandy, A., Humphrey, M., Isaac, P., Polidori, D., Reichstein, M., Ribeca, A., van Ingen, C., Vuichard, N., Zhang, L., Amiro, B., Ammann, C., Arain, M. A., Ardö, J., Arkebauer, T., Arndt, S. K., Arriga, N., Aubinet, M., Aurela, M., Baldocchi, D., Barr, A., Beamesderfer, E., Marchesini, L. B., Bergeron, O., Beringer, J., Bernhofer, C., Berveiller, D., Billesbach, D., Black, T. A., Blanken, P. D., Bohrer, G., Boike, J., Bolstad, P. V., Bonal, D., Bonnefond, J.-M., Bowling, D. R., Bracho, R., Brodeur, J., Brümmer, C., Buchmann, N., Burban, B., Burns, S. P., Buysse, P., Cale, P., Cavagna, M., Cellier, P., Chen, S., Chini, I., Christensen, T. R., Cleverly, J., Collalti, A., Consalvo, C., Cook, B. D., Cook, D., Coursolle, C., Cremonese, E., Curtis, P. S., D’Andrea, E., da Rocha, H., Dai, X., Davis, K. J., Cinti, B. D., Grandcourt, A. de, Ligne, A. D., De Oliveira, R. C., Delpierre, N., Desai, A. R., Di Bella, C. M., Tommasi, P. di, Dolman, H., Domingo, F., Dong, G., Dore, S., Duce, P., Dufrêne, E., Dunn, A., Dušek, J., Eamus, D., Eichelmann, U., ElKhidir, H. A. M., Eugster, W., Ewenz, C. M., Ewers, B., Famulari, D., Fares, S., Feigenwinter, I., Feitz, A., Fensholt, R., Filippa, G., Fischer, M., Frank, J., Galvagno, M., et al.: The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data, Sci. Data, 7, 225, https://doi.org/10.1038/s41597-020-0534-3, 2020.

Perez-Priego, O., Katul, G., Reichstein, M., El-Madany, T. S., Ahrens, B., Carrara, A., Scanlon, T. M., and Migliavacca, M.: Partitioning Eddy Covariance Water Flux Components Using Physiological and Micrometeorological Approaches, J. Geophys. Res. Biogeosciences, 123, 3353–3370, https://doi.org/10.1029/2018JG004637, 2018.
Pohl, F., Rakovec, O., Rebmann, C., Hildebrandt, A., Boeing, F., Hermanns, F., Attinger, S., Samaniego, L., and Kumar, R.: Long-term daily hydrometeorological drought indices, soil moisture, and evapotranspiration for ICOS sites, Sci. Data, 10, 281, https://doi.org/10.1038/s41597-023-02192-1, 2023.

Rana, G., Palatella, L., Scanlon, T. M., Martinelli, N., and Ferrara, R. M.: CO2 and H2O flux partitioning in a Mediterranean cropping system, Agric. For. Meteorol., 260–261, 118–130, https://doi.org/10.1016/j.agrformet.2018.06.007, 2018.

Rannik, Ü. and Vesala, T.: Autoregressive filtering versus linear detrending in estimation of fluxes by the eddy covariance method, Bound.-Layer Meteorol., 91, 259–280, https://doi.org/10.1023/A:1001840416858, 1999.

Rebmann, C., Aubinet, M., Schmid, H., Arriga, N., Aurela, M., Burba, G., Clement, R., De Ligne, A., Fratini, G., Gielen, B., Grace, J., Graf, A., Gross, P., Haapanala, S., Herbst, M., Hörtnagl, L., Ibrom, A., Joly, L., Kljun, N., and Franz, D.: ICOS eddy covariance flux-station site setup: A review, Int. Agrophysics, 32, 471–494, https://doi.org/10.1515/intag-2017-0044, 2018.

Reichstein, M., Falge, E., Baldocchi, D., Papale, D., Aubinet, M., Berbigier, P., Bernhofer, C., Buchmann, N., Gilmanov, T., Granier, A., Grünwald, T., Havránková, K., Ilvesniemi, H., Janous, D., Knohl, A., Laurila, T., Lohila, A., Loustau, D., Matteucci, G., Meyers, T., Miglietta, F., Ourcival, J.-M., Pumpanen, J., Rambal, S., Rotenberg, E., Sanz, M., Tenhunen, J., Seufert, G., Vaccari, F., Vesala, T., Yakir, D., and Valentini, R.: On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm, Glob. Change Biol., 11, 1424–1439, https://doi.org/10.1111/j.1365-2486.2005.001002.x, 2005.

Reichstein, M., Stoy, P. C., Desai, A. R., Lasslop, G., and Richardson, A. D.: Partitioning of Net Fluxes, in: Eddy Covariance: A Practical Guide to Measurement and Data Analysis, edited by: Aubinet, M., Vesala, T., and Papale, D., Springer Netherlands, Dordrecht, 263–289, https://doi.org/10.1007/978-94-007-2351-1_9, 2012.

Sabbatini, S., Mammarella, I., Arriga, N., Fratini, G., Graf, A., Hörtnagl, L., Ibrom, A., Longdoz, B., Mauder, M., Merbold, L., Metzger, S., Montagnani, L., Pitacco, A., Rebmann, C., Sedlák, P., Šigut, L., Vitale, D., and Papale, D.: Eddy covariance raw data processing for CO2 and energy fluxes calculation at ICOS ecosystem stations, Int. Agrophysics, 32, 495–515, https://doi.org/10.1515/intag-2017-0043, 2018.

Scanlon, T. M. and Albertson, J. D.: Turbulent transport of carbon dioxide and water vapor within a vegetation canopy during unstable conditions: Identification of episodes using wavelet analysis, J. Geophys. Res. Atmospheres, 106, 7251–7262, https://doi.org/10.1029/2000JD900662, 2001.

Scanlon, T. M. and Kustas, W. P.: Partitioning carbon dioxide and water vapor fluxes using correlation analysis, Agric. For. Meteorol., 150, 89–99, https://doi.org/10.1016/j.agrformet.2009.09.005, 2010.

Scanlon, T. M. and Kustas, W. P.: Partitioning Evapotranspiration Using an Eddy Covariance-Based Technique: Improved Assessment of Soil Moisture and Land–Atmosphere Exchange Dynamics, Vadose Zone J., 11, vzj2012.0025, https://doi.org/10.2136/vzj2012.0025, 2012.

Scanlon, T. M. and Sahu, P.: On the correlation structure of water vapor and carbon dioxide in the atmospheric surface layer: A basis for flux partitioning, Water Resour. Res., 44, https://doi.org/10.1029/2008WR006932, 2008.

Scanlon, T. M., Schmidt, D. F., and Skaggs, T. H.: Correlation-based flux partitioning of water vapor and carbon dioxide fluxes: Method simplification and estimation of canopy water use efficiency, Agric. For. Meteorol., 279, 107732, https://doi.org/10.1016/j.agrformet.2019.107732, 2019.
925

Schaller, C., Göckede, M., and Foken, T.: Flux calculation of short turbulent events &ndash; comparison of three methods, Atmospheric Meas. Tech., 10, 869–880, https://doi.org/10.5194/amt-10-869-2017, 2017.

Sousa, P. M., Barriopedro, D., García-Herrera, R., Ordóñez, C., Soares, P. M. M., and Trigo, R. M.: Distinct influences of large-scale circulation and regional feedbacks in two exceptional 2019 European heatwaves, Commun. Earth Environ., 1, 1–13, https://doi.org/10.1038/s43247-020-00048-9, 2020.

930

Strunin, M. A. and Hiyama, T.: Applying wavelet transforms to analyse aircraft-measured turbulence and turbulent fluxes in the atmospheric boundary layer over eastern Siberia, Hydrological Processes, 18, 3081–3098, https://doi.org/10.1002/hyp.5750, 2004.

Sulman, B. N., Roman, D. T., Scanlon, T. M., Wang, L., and Novick, K. A.: Comparing methods for partitioning a decade of carbon dioxide and water vapor fluxes in a temperate forest, Agric. For. Meteorol., 226–227, 229–245, https://doi.org/10.1016/j.agrformet.2016.06.002, 2016.

935

Super, I., Denier van der Gon, H. A. C., Visschedijk, A. J. H., Moerman, M. M., Chen, H., van der Molen, M. K., and Peters, W.: Interpreting continuous in-situ observations of carbon dioxide and carbon monoxide in the urban port area of Rotterdam, Atmospheric Pollution Research, 8, 174–187, https://doi.org/10.1016/j.apr.2016.08.008, 2017.

Thomas, C. and Foken, T.: Re-evaluation of integral turbulence characteristics and their parameterisations, in: 15th Symposium on Boundary Layers and Turbulence, Boston, 129–132, 2002.

940

Thomas, C., Martin, J. G., Goeckede, M., Siqueira, M. B., Foken, T., Law, B. E., Loescher, H. W., and Katul, G.: Estimating daytime subcanopy respiration from conditional sampling methods applied to multi-scalar high frequency turbulence time series, Agric. For. Meteorol., 148, 1210–1229, https://doi.org/10.1016/j.agrformet.2008.03.002, 2008.

Torrence, C. and Compo, G. P.: A Practical Guide to Wavelet Analysis, Bull. Am. Meteorol. Soc., 79, 61–78, https://doi.org/10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2, 1998.

945

Twardosz, R., Walanus, A., and Guzik, I.: Warming in Europe: Recent Trends in Annual and Seasonal temperatures, Pure Appl. Geophys., 178, 4021–4032, https://doi.org/10.1007/s00024-021-02860-6, 2021.

Valentini, R., De ANGELIS, P., Matteucci, G., Monaco, R., Dore, S., and Mucnozza, G. E. S.: Seasonal net carbon dioxide exchange of a beech forest with the atmosphere, Glob. Change Biol., 2, 199–207, https://doi.org/10.1111/j.1365-2486.1996.tb00072.x, 1996.

950

Vekuri, H., Tuovinen, J.-P., Kulmala, L., Papale, D., Kolari, P., Aurela, M., Laurila, T., Liski, J., and Lohila, A.: A widely-used eddy covariance gap-filling method creates systematic bias in carbon balance estimates, Sci. Rep., 13, 1720, https://doi.org/10.1038/s41598-023-28827-2, 2023.

Wang, W., Smith, J. A., Ramamurthy, P., Baeck, M. L., Bou-Zeid, E., and Scanlon, T. M.: On the correlation of water vapor and CO2: Application to flux partitioning of evapotranspiration, Water Resour. Res., 52, 9452–9469, https://doi.org/10.1002/2015WR018161, 2016.

955

Warm Winter 2020 Team, & ICOS Ecosystem Thematic Centre: Warm Winter 2020 ecosystem eddy covariance flux product for 73 stations in FLUXNET-Archive format—release 2022-1 (Version 1.0), ICOS Carbon Portal, https://doi.org/10.18160/2G60-ZHAK, 2020.

Webb, E. K., Pearman, G. I., and Leuning, R.: Correction of flux measurements for density effects due to heat and water vapour transfer, Q. J. R. Meteorol. Soc., 106, 85–100, https://doi.org/10.1002/qj.49710644707, 1980.

960
Wehr, R., Munger, J. W., McManus, J. B., Nelson, D. D., Zahniser, M. S., Davidson, E. A., Wofsy, S. C., and Saleska, S. R.: Seasonality of temperate forest photosynthesis and daytime respiration, Nature, 534, 680–683, https://doi.org/10.1038/nature17966, 2016.

Wilczak, J. M., Oncley, S. P., and Stage, S. A.: Sonic Anemometer Tilt Correction Algorithms, Bound.-Layer Meteorol., 99, 127–150, https://doi.org/10.1023/A:1018966204465, 2001.

Wutzler, T., Lucas-Moffat, A., Migliavacca, M., Knauer, J., Sickel, K., Šigut, L., Menzer, O., and Reichstein, M.: Basic and extensible post-processing of eddy covariance flux data with REddyProc, Biogeosciences, 15, 5015–5030, https://doi.org/10.5194/bg-15-5015-2018, 2018.

Zahn, E., Bou-Zeid, E., Good, S. P., Katul, G. G., Thomas, C. K., Ghannam, K., Smith, J. A., Chamecki, M., Dias, N. L., Fuentes, J. D., Alfieri, J. G., Kwon, H., Caylor, K. K., Gao, Z., Soderberg, K., Bambach, N. E., Hipps, L. E., Prueger, J. H., and Kustas, W. P.: Direct partitioning of eddy-covariance water and carbon dioxide fluxes into ground and plant components, Agric. For. Meteorol., 315, 108790, https://doi.org/10.1016/j.agrformet.2021.108790, 2022.

Zeeman, M. J., Eugster, W., and Thomas, C. K.: Concurrency of Coherent Structures and Conditionally Sampled Daytime Sub-canopy Respiration, Bound.-Layer Meteorol., 146, 1–15, https://doi.org/10.1007/s10546-012-9745-2, 2013.
Supplementary materials for: Improvement of CO₂ flux quality through wavelet-based Eddy Covariance: a new method for partitioning respiration and photosynthesis

Pedro Henrique H. Coimbra<fnref n="1" />*, Benjamin Loubet<fnref n="1" />, Olivier Laurent<fnref n="2" />, Matthias Mauder<fnref n="3" />,<fnref n="4" />, Bernard Heinesch<fnref n="5" />, Jonathan Bitton<fnref n="5" />, Nicolas Delpierre<fnref n="6" />,<fnref n="7" />, Jérémie Depuydt<fnref n="1" />, Pauline Buysse<fnref n="1" />

A Demonstrating covariance can be calculated using decomposed signals

Assuming two time series variables, x and y, the sum of the $\overline{x'y'}(j)$ for all frequencies, j is not strictly equal to the covariance $\overline{x'y'}$. Indeed, we have:

$$ \overline{x'y'} = \sum_{j=0}^{J} x'_{j} \sum_{j=0}^{J} y'_{j} = \sum_{j=0}^{J} \overline{x'_{j}y'_{j}} + \sum_{k \neq j} \overline{x'_{k}y'_{k}} \quad (A.1) $$

For discrete wavelet transform DWT, the orthogonality of the wavelets base implies independent frequencies, i.e. $\overline{x'_{k}y'_{j}} = 0$ for $k \neq j$, hence $\overline{x'_{j}y'_{j}} = \sum_{j=1}^{J} \overline{x'_{j}y'_{j}}$. For continuous wavelet transform CWT, a coefficient $C_{\phi}$ is introduced to ensure energy conservation and correct for cross-correlations of x and y between scales j, leading to:

$$ \overline{x'y'} = C_{\phi} \sum_{j} \overline{x'_{j}y'_{j}} \quad (A.2) $$

$C_{\phi}$ depends on the wavelet chosen (Table 1) Alternatively, a direct formulation of the covariance was proposed by (Torrence and Compo, 1998) based on $\tilde{f}_{x}(n, j)$ and $\tilde{f}_{y}(n, j)$ the wavelet coefficients for time series x and y:

$$ \overline{x'y'} = \frac{\delta j \delta t}{C_{\delta} N} \sum_{n=0}^{N-1} \sum_{j=0}^{J} \tilde{f}_{x}(n, j) \tilde{f}_{y}(n, j) \quad (A.3) $$

Where $C_{\delta}$ is a scale-independent reconstruction factor depending on the chosen mother wavelet function (Table 1).

<footnote n="1">¹ ECOSYS, INRAE, AgroParisTech, Université Paris-Saclay, Palaiseau, France</footnote>

<footnote n="2">² Laboratoire des Sciences du Climat et de l'Environnement, CEA, CNRS, Université Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote n="3">³ Institute of Meteorology and Climate Research - Atmospheric Environmental Research (IMK-IFU), Karlsruhe Institute of Technology, Garmisch-Partenkirchen, Germany</footnote>

<footnote n="4">⁴ Institute of Hydrology and Meteorology, Technische Universität Dresden, Dresden, Germany</footnote>

<footnote n="5">⁵ Faculté des Sciences Agronomiques de Gembloux, Unité de Physique, Gembloux, Belgium</footnote>

<footnote n="6">⁶ Ecologie Systématique Evolution, CNRS, AgroParisTech, Université Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote n="7">⁷ Institut Universitaire de France, France</footnote>

<footnote>* corresponding author: pedro-henrique.herig-coimbra@inrae.fr</footnote>
<figure><img src="image_12.png" /></figure>

Figure S 1. Covariance is calculated using the standard equation (w'c' in black), DWT considering cross-correlation (Σw'jΣc_j' in blue), and ignoring it (Σw'jc_j' in orange), idem for CWT plus without C_φ (dotted orange), using covariance equation from Torrence and Compo (1998) (CWT<sub>TC1998</sub>). Top 20 Hz data before time averaging, bottom half-hour average. Data for FR-Gri 03/05/2022.
<figure><img src="image_13.png" /><figcaption>Figure S 2. Half-hourly NEE calculated using EC<sub>S</sub> and DW-EC grouped by quality flags. high: well developed turbulence (ITC < 30 %) and stationary (STA < 30 %), medium: at least one of the tests higher than 30 % but both lower than 100%, low: at least one of the tests higher than 100 %. In grey 1:1 and in orange linear relation. No gap filling was used.</figcaption></figure>

1010
<figure><img src="image_14.png" /><figcaption>Figure S 3. Monthly mean cospectra of w'CO₂' (φ_c), its positive and negative parts, and modelled cospectra following Horst (1997) grouped by stratification status. Cospectra curves sum to 1. Horst (1997) cospectra are calculated using measured mean wind speed, displacement height and Obukhov length. Peak frequencies are shown with an arrow.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure S 4. Correlation Matrix for GPP (orange, bottom-left) and R<sub>eco</sub> (blue, upper-right) estimations, sites combined. NT and DT partitioning methods were used with NEE calculated using standard EC (ECNT and ECDT, respectively) and discrete wavelets (DWNT and DWDT).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure S 5. (a) Parameters estimated for night-time partitioning method (NT). (b) Parameters estimated for day-time partitioning method (DT).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure S 6. Correlation Matrix for GPP (orange, bottom-left) and Reco (blue, upper-right) estimations, sites combined, using the following partitioning methods: NT, DT (with light inhibition), DW-CS (DWCS), φ<sub>c</sub>, and ECCS (same partitioning as DWCS but used for w'χ<sub>CO₂</sub>' without wavelet decomposition). The diagonal shows flux distribution.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

Figure S 7. Ratios of daily ecosystem respiration ($R_{eco}$) estimated by a standard night-time modeling approach (NT) and measured by discrete wavelet conditional sampling (DW). Symbols indicate ratios of daily NT and DW for the years 2019–2022; the red line represents the block-average (window size = 14 days) with shaded bands indicating ±1 standard error. Horizontal line represents unity of the ratio. Vertical dotted lines represent the start or end of the season (calendar years for forest site FR-Fon and crop season for FR-Gri).

1035
<figure><img src="image_19.png" /><figcaption>Figure S 8. Monthly averaged half-hourly GPP and R<sub>eco</sub> estimations using NT, DT and DW methods. The darker region indicates interquartile (25<sup>th</sup> and 75<sup>th</sup> percentile), and the lighter region with dotted lines indicates the 5<sup>th</sup> and 95<sup>th</sup> percentile.</figcaption></figure>
<figure><img src="image_20.png" /></figure>

Figure S 9. Seasonally cumulated sum for NEE in Figure 6 (a) including additionally discrete-wavelet-based NEE forcing same gaps as EC<sub>S</sub> (DW'-EC).