Diurnal variability of atmospheric O₂, CO₂, and their exchange ratio above a boreal forest in southern Finland

The exchange ratio (ER) between atmospheric O₂ and CO₂ is a useful tracer for better understanding the carbon budget on global and local scales. The variability of ER (in mol O₂ per mol CO₂) between terrestrial ecosystems is not well known, and there is no consensus on how to derive the ER signal of an ecosystem, as there are different approaches available, either based on concentration (ER_atmos) or flux measurements (ER_forest). In this study we measured atmospheric O₂ and CO₂ concentrations at two heights (23 and 125 m) above the boreal forest in Hyytiälä, Finland. Such measurements of O₂ are unique and enable us to potentially identify which forest carbon loss and production mechanisms dominate over various hours of the day. We found that the ER_atmos signal at 23 m not only represents the diurnal cycle of the forest exchange but also includes other factors, including entrainment of air masses in the atmospheric boundary layer before midday, with different thermodynamic and atmospheric composition characteristics. To derive ER_forest, we infer O₂ fluxes using multiple theoretical and observation-based micro-meteorological formulations to determine the most suitable approach. Our resulting ER_forest shows a distinct difference in behaviour between daytime (0.92 ± 0.17 mol mol⁻¹) and nighttime (1.03 ± 0.05 mol mol⁻¹). These insights demonstrate the diurnal variability of different ER signals above a boreal forest, and we also confirmed that the signals of ER_atmos and ER_forest cannot be used interchangeably. Therefore, we recommend measurements on multiple vertical levels to derive O₂ and CO₂ fluxes for the ER_forest signal instead of a single level time series of the concentrations for the ER_atmos signal. We show that ER_forest can be further split into specific signals for respiration (1.03 ± 0.05 mol mol⁻¹) and photosynthesis (0.96 ± 0.12 mol mol⁻¹). This estimation allows us to separate the net ecosystem exchange (NEE) into gross primary production (GPP) and total ecosystem respiration (TER), giving comparable results to the more commonly used eddy covariance approach. Our study shows the potential of using atmospheric O₂ as an alternative and complementary method to gain new insights into the different CO₂ signals that contribute to the forest carbon budget.

1 Introduction

To understand how the increasing carbon dioxide (CO₂) levels in the atmosphere will change our climate, we need to know the sources and sinks of CO₂ separately. The main sources are fossil fuel combustion and land-use change, and the main sinks are the net uptake by the terrestrial biosphere and the oceans (Friedlingstein et al., 2022). The net terrestrial biospheric sink (net ecosystem exchange, NEE) results from many fluxes of which the two largest are typically gross primary production (GPP) and the total ecosystem respiration (TER). Knowing these gross fluxes separately will allow for better estimates of the changing behaviour of the biosphere carbon sink, as GPP and TER respond differently to climate change and increasing atmospheric CO₂ levels (Cox et al., 2013; Ballantyne et al., 2012).

Using tracers in addition to CO₂ allows us to gain further insights into GPP and TER, without relying on a temperature-based function to parameterize TER as is used for eddy covariance (EC) measurements (e.g. Reichstein et al., 2005). Tracers such as atmospheric O₂ (Keeling and Manning, 2014), as well as COS, δ¹³C, or Δ¹⁷O, have the important advantage of sharing a process or pathway with CO₂ directly (Wehr et al., 2016; Whelan et al., 2018; Peters et al., 2018; Koren et al., 2019; Kooijmans et al., 2021). This allows one to use numerical models to test formulations of processes, such as stomatal and mesophyll exchange, photosynthesis, pool-specific respiration, and even turbulent canopy exchange. Atmospheric O₂ is directly coupled to CO₂ in several processes through the so-called exchange ratio (ER) (Keeling and Manning, 2014; Manning and Keeling, 2006; Keeling et al., 1993). This ER indicates the number of moles of O₂ that are consumed per moles of CO₂ that are produced (or vice versa) and gives a process-specific signature (Keeling, 1988).

On the global scale, the O₂:CO₂ molar ratio ER has been used to derive the global oceanic CO₂ sink and determine the global carbon budget (Stephens et al., 1998; Rödenbeck et al., 2008; Tohjima et al., 2019). This is done by solving the atmospheric budgets of O₂ and CO₂ with the following equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{d}t}}=F-O-B,\text{(2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{\mathrm{O}}_{\mathrm{2}}}{\mathrm{d}t}}=-{\mathit{\alpha }}_{\mathrm{F}}F+{\mathit{\alpha }}_{\mathrm{B}}B+Z_{{\mathrm{O}}_{\mathrm{2}}},\end{array}$$

where F is the fossil fuel CO₂ emissions, O is ocean CO₂ uptake, B is the net terrestrial biosphere sink of CO₂, and $Z_{{\mathrm{O}}_{\mathrm{2}}}$ indicates the ocean O₂ outgassing. Symbols α_F and α_B indicate the global ERs for fossil fuel combustion and the net terrestrial biosphere sink, respectively. In these global studies simplified global average values are used for α_F and α_B, where α_F is determined from the global mixture of fuels burned, which results in 1.38 [mol mol⁻¹] (Keeling and Manning, 2014), and α_B was determined by laboratory measurements and a literature study of different plant and soil materials, which resulted in 1.1 [mol mol⁻¹] (Severinghaus, 1995). Furthermore, α_B is used to combine O₂ and CO₂ into atmospheric potential oxygen (APO) (Stephens et al., 1998), which is used in determining the ocean carbon sink, and recently it has also shown to be a suitable tracer to detect fossil fuel emission reductions during the COVID-19 pandemic (Pickers et al., 2022). For these larger-scale applications using APO, it is important to have good estimates for the terrestrial biosphere ERs.

Figure 1 Schematic overview of the different O₂:CO₂ exchange ratio (ER) signals, measured and analysed in and above a forest, influenced by the different O₂ and CO₂ fluxes and meteorological processes (a), together with a more detailed look on which processes influence the different ER signals (b). Panel (a) shows the direction of the surface fluxes during the day in the surface layer, which includes the roughness sublayer and the inertial sublayer. During the night the direction of the O₂ and the CO₂ surface fluxes are the other way around. The ER of the atmosphere (ER_atmos) is determined from the change over time (Δ_(t)) in the O₂ and CO₂ concentration measurements, and the ER of the forest (ER_forest) is calculated from the surface fluxes of O₂ and CO₂, which are inferred from (∼) the vertical gradient (Δ_(z)). ER_a represents assimilation processes that influence the gross primary production (GPP) flux, and ER_r represents respiration processes that influence the total ecosystem respiration (TER) flux. Panel (b) shows the connections between the processes, measurements, and the ERs. Dotted lines indicate smaller influences of the processes that are connected to it compared to solid lines.

On local/ecosystem scales, previous studies have shown that this terrestrial biosphere ER is not a constant value of 1.1 as used on the global scale and that it shows a certain degree of temporal and spatial variability. These studies either measured the oxidative ratios (ORs) from elemental composition analysis (Worrall et al., 2013; Randerson et al., 2006; Gallagher et al., 2017) or derived the ER from atmospheric concentration measurements (Battle et al., 2019; Seibt et al., 2004; van der Laan et al., 2014). Note that there is a distinction in the terminology between ER and OR. The OR indicates the stoichiometry of specific materials, whereas the ER indicates the exchange between the atmosphere and organisms or ecosystems. By using elemental composition analysis, the OR reflects the relationship between O₂ and CO₂ over a longer timescale of years or decades and only reflects the OR from the materials that are sampled. By using atmospheric concentration measurements for the ER, the ER reflects a shorter timescale compared to the OR of hourly to daily time periods, and it also reflects a different spatial scale, as the ER includes all processes that are originating from the footprint. The spatial scale that is covered by the ER signal depends on the type of measurements or modelling, i.e. leaf, canopy, or ecosystem. Both the OR- and ER-based studies showed that O₂:CO₂ molar ratio of the biosphere changes per ecosystem and over different time periods. The ER from the gas exchange experiments can furthermore be used for the separation of GPP and TER, using a specific ecosystem ER, which are determined with two alternative approaches (see Fig. 1) (Seibt et al., 2004; Stephens et al., 2007; Ishidoya et al., 2013, 2015; Battle et al., 2019). The first is the ER of the atmosphere (ER_atmos), which is the ratio of the evolution of the atmospheric O₂ and CO₂ concentration measurements over time, and the second is the ER of the forest (ER_forest), which is the ratio of the net surface fluxes of O₂ and CO₂ above the canopy, including all processes occurring below the canopy, including both vegetation and soil exchange. First attempts to estimate ER_forest were made using one-box models (Seibt et al., 2004; Ishidoya et al., 2013). More accurate estimates of ER_forest would be based on in situ measured O₂ and CO₂ surface fluxes; however, O₂ currently cannot yet be measured accurately using EC techniques. Ishidoya et al. (2015) showed the first surface fluxes of O₂ using vertical gradients of O₂, an alternative technique to EC, and CO₂ measurements at two heights above the canopy in the surface layer in a temperate forest in Japan. Their results showed that the ER_forest signal could be used to separate the NEE signal into GPP and TER, consistent with the separation method for EC measurements using an empirical function of air temperature.

When using O₂ to separately estimate GPP and TER fluxes, it is important to use the value for ER that represents ecosystem exchange. Seibt et al. (2004) showed that the signal of ER_atmos cannot be directly linked to the exchange of carbon in the terrestrial biosphere, because in addition to the biosphere ER_atmos is also affected by advection, boundary layer dynamics, and entrainment (Fig. 1). In contrast, Ishidoya et al. (2015) found similar values for ER_atmos and ER_forest. So far, there is no clear consensus on which signal should be used to indicate the ER of the ecosystem. Furthermore, since atmospheric O₂ measurements are challenging to make, only a few studies exist that measured atmospheric O₂ from flasks (Seibt et al., 2004) or continuously (Ishidoya et al., 2015; Stephens et al., 2007; Battle et al., 2019) above an ecosystem and derive ER signals. The uncertainty and spatial and temporal variability of O₂:CO₂ molar ratio of the biosphere are therefore not well known (Manning and Keeling, 2006; Keeling and Manning, 2014), and knowledge about the difference between ER_forest and ER_atmos, its variability across difference regions and ecosystems, and how ER_forest can be used on both the local and global scale to advance our understanding of the carbon cycle is still limited. Therefore, more and longer in situ time series of atmospheric O₂ measurements are needed, and further understanding of O₂ and CO₂ exchange above and below the canopy is crucial to continue the pioneering work by Seibt et al. (2004), Stephens et al. (2007), Ishidoya et al. (2015), and Battle et al. (2019) and to improve the application of the global biosphere ER, resulting in a better understanding of the carbon balance on local, regional, and global scales.

The aim of this study is to improve upon existing methods to calculate ER_forest and get a better comparison of the ER_atmos and ER_forest signals. We carried out a measurement campaign in Hyytiälä, Finland, for two short periods in spring/summer 2018 and 2019 where both O₂ and CO₂ were measured at two heights with a setup including a differential fuel cell analyser for O₂. We used our measurements to determine the diurnal behaviour of the relation between the concentrations and the fluxes of O₂ and CO₂, by using either one or both measurement heights on the tower. The objectives of this study are the following: (1) to extend the existing continuous O₂ records, (2) to calculate the O₂ surface fluxes in a boreal forest for the first time, (3) to combine the O₂ and the CO₂ fluxes (to calculate ER_forest from these fluxes) and to compare the ER_atmos and ER_forest signals, and (4) to use ER_forest to estimate GPP and TER fluxes.

In this paper, we first describe the measurement site, experimental setup, and methods used to derive O₂ fluxes and the different ER signals (Sect. 2). We present the measurements for the whole campaign, and we select a representative day to determine the most suitable approach for deriving O₂ fluxes and to determine ER_forest (Sect. 3). A detailed evaluation and discussion of our ER_atmos and ER_forest signals is given in Sect. 4. We finalize with our conclusion about the diurnal variability of the ER signals for a representative day of a boreal forest (Sect. 5).

2 Methods

To determine ER_atmos and ER_forest (and its diurnal variability), we measured O₂ and CO₂ continuously at two heights above a boreal forest during two short campaigns at Hyytiälä. These OXHYYGEN (oxygen at Hyytiälä) campaigns took place in the spring/summer of 2018 (3 June through 2 August) and 2019 (10 June through 17 July). In this section, we describe the measurement site and instrumental setup, as well as the methods used to determine the O₂ and CO₂ fluxes from the measured vertical gradient and the ER signals.

2.1 Measurement site

The measurements were made at Hyytiälä SMEAR II forestry station of the University of Helsinki in Finland (61^∘51^′ N, 24^∘17^′ E; +181 $\mathrm{m}\mathrm{a}.\mathrm{m}.\mathrm{s}.\mathrm{l}.$; time zone: UTC+2); this site is described in more detail in, for example, Hari et al. (2013). The SMEAR II station is a boreal site within the European Integrated Carbon Observation System (ICOS) network with atmospheric and ecosystem measurements. The SMEAR II station is located inside a homogeneous forest of Scots pine trees (Pinus Sylvestris) with a dominant canopy height of 18 m and some silver birch and aspen trees. The forest floor is covered with mosses and herbs. The soils are podzols on top of glacial till. A large lake is located close (around 550 m) to the measurement site and has a fetch of 250 m over the dominant wind direction of 230^∘. The footprint of the site is mostly influenced by natural sources, with the atmospheric signal dominated by forest exchange (Carbon Portal ICOS RI, 2022). The measurement site includes several towers, including a 128 m tall tower and a 23 m high walk-up tower, where atmospheric variables and gas concentrations are continuously measured. The operational data from this tower are publicly available online at http://avaa.tdata.fi/web/smart/smear/ (last access: 5 January 2022). Our O₂ and CO₂ measurement setup was installed in a cabin at the bottom of the 23 m high tower, and air was sampled from aspirated inlets (Blaine et al., 2006), installed at 23 m in the smaller tower and at 125 m in the tall tower, which are 5 and 107 m above the canopy height, respectively. We used both levels to calculate the vertical gradient for the flux calculations (Sect. 2.3).

Figure 2 Schematic overview of the measurement setup used at Hyytiälä. The setup includes an Oxzilla O₂ fuel cell analyser and a URAS26 NDIR CO₂ analyser. The system measured air sampled from two heights of either 23 or 125 m.

2.2 Experimental setup

The measurement setup is based on the instrument used in van Leeuwen and Meijer (2015), following the methods in van der Laan-Luijkx et al. (2010) and Stephens et al. (2007). O₂ is measured with a Sable Systems Oxzilla II fuel-cell-based instrument, and CO₂ is measured with an ABB continuous gas analyser URAS26, which is a non-dispersive infrared (NDIR) photometer. The gas-handling schematic is shown in Fig. 2.

Air was pumped from either 23 or 125 m height to the measurement system at the base of the tower. Both inlet lines were continuously flushed, where either one of the heights is measured by the system with a sample flow of around 120 mL min⁻¹ and the other flushed to the room with a higher flow rate of around 2 L min⁻¹, which allows for fast switching between the two heights. We switched between the inlets every half hour to match the EC measurements of ICOS that were already present in the tower and to get a more stable signal of O₂. The air of the selected inlet was first cooled to −60 ^∘C with a cryogenic cooler to remove water vapour from the air, before entering the system. Second-stage drying of the air streams was done with magnesium perchlorate (Mg(ClO₄)₂) traps. The sample air was continuously measured against a reference gas (differentially for O₂ and alternatively for CO₂), and the pressure in both sample and reference lines was matched to be the same using a pressure control system (MKS Instruments, types 223B, 248A, and 250E for the pressure transducers, regulating valve, and control system, respectively). The reference and sample lines were switched every 2 min between the two fuel cells in the Oxzilla analyser. We measured a set of three calibration cylinders and one target cylinder every 23 h for half an hour per cylinder.

The measurements of these calibration gases allowed for calibration of our measurements against the international Scripps Institution of Oceanography (SIO) scale for $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$. We did that by using cylinders that are filled in the laboratory at the University of Groningen, where they were calibrated with the primary Scripps cylinders (Nguyen et al., 2022). The O₂ measurements are normally expressed as $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratios in “per meg” instead of mole fraction (ppm), since O₂ is not a trace gas because of its high abundance of 20.95 %; therefore, the mole fraction varies due to changes in other gases, such as CO₂ (Keeling et al., 1998). So $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ is defined as

$$\begin{array}{}\text{(3)}& \mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})=\left({\displaystyle \frac{({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}{)}_{\text{sample}}}{({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}{)}_{\text{reference}}}}-\mathrm{1}\right)\cdot {\mathrm{10}}^{\mathrm{6}}\left[\mathrm{per}\mathrm{meg}\right].\end{array}$$

For simplicity, in this paper we use the term O₂ instead of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, and we use the term “concentration” rather than “mole fraction” when discussing both CO₂ and O₂. Equation (3) indicates a change compared to a reference level. Negative values therefore indicate concentrations of O₂ lower than the reference value. To allow for comparison of changes in CO₂ and O₂ directly, we converted the units of O₂ from per meg to ppm equivalents (ppmEq), where a change of 1 ppm CO₂ corresponds to a 4.77 per meg change in O₂ (Tohjima et al., 2005; Kozlova and Manning, 2009).

We modified the method described in van der Laan-Luijkx et al. (2010) to calibrate the measurements. The raw CO₂ measurements have a frequency of one measurement per 6 s, the raw O₂ measurements have a frequency of one measurement per second, and both give one value every 4 min in the form of ΔCO₂ and Δ(Δ)O₂, respectively. CO₂ is measured on a single cell instrument; therefore, ΔCO₂ is the difference between the 2 min averages of the sample air (S) and the reference cylinder (R), giving (S−R). For the 2 min averaged CO₂ measurements, the last 78 s of each 2 min period were used. Note that for CO₂, the NDIR system is different compared to other systems used and therefore does not need a zero gas (Pickers et al., 2017). O₂ is measured on a double-cell instrument and therefore gives a double differential signal. The Δ(Δ)O₂ value is the difference between the 2 min averaged difference between S and R and the 2 min averaged difference between R and S ($(S-R)-(R-S)$). For the 2 min averaged O₂ values, the last 100 s of each 2 min period are used. In 2019, the MKS pressure control valve was not functioning optimally, which led to a small instability in the differential pressure between the sample and reference lines. We therefore corrected the 4 min values of Δ(Δ)O₂ for this deviation measured by the MKS differential pressure sensor (PMKS) by multiplying Δ(Δ)O₂ with (0.095 × PMKS), which we derived based on the measurements of the calibration cylinders. In 2018, there was no instability in the pressure control valve; therefore, no correction was applied in that year. The PMKS deviations correlated with temperature and increased towards the end of the 2019 campaign. Figure B1 in Appendix B shows that the highest corrections were made during midday at the end of the campaign. The O₂ vertical gradient is hardly affected by the correction as it is the difference between measurements at two heights that are both undergoing the same deviation.

For both CO₂ and O₂, the 4 min values were subsequently used to calculate half-hourly means, where we excluded the first 4 min value after the heights are switched, together with the measurements that did not fall inside the boundary based on the median absolute deviation (MAD) (Rousseeuw and Verboven, 2002). For every half-hourly mean, a standard error is calculated (see Eq. 4), which is used in further analysis to determine the uncertainty of our measurements.

The linear calibration response functions for both O₂ and CO₂ were calculated for every measurement period of the calibration cylinders, which was about every 23 h. For the response functions, we used a constant slope based on the mean of all the calibration slopes measured in the specific year. The y intercept values of the response functions were interpolated to the time of the measurement, based on the two calibrations bracketing the measurement time. To facilitate the comparison of the O₂ and CO₂ measurements of the two heights and allow for flux calculations based on the vertical gradient, we interpolated the data to one measurement for every 30 min for each height. Based on the target cylinders, measured during the calibration period, the stability of the long-term measurements was determined (Table 1). A different target cylinder, with different composition of air for 2019 compared to 2018, was used, which resulted in different outcomes for the standard deviation (SD) and the mean difference for these periods. The mean difference is calculated from the target measurements at Hyytiälä compared to the calibrated values using the SIO cylinders in Groningen. The measurement period of 2018 was also longer and so more points were included for the SD and mean difference calculations. The long-term measurement precision of this device throughout the duration of the two measurement campaigns compared to the recommendations of the World Meteorological Organization (WMO) will be further discussed in Sect. 4.1.

Table 1 The mean difference and the standard deviation (SD) of the target cylinder measurements of O₂ and CO₂ for the 2018 and 2019 periods separately, together with the number of data points used to calculate these specific values.

2.3 Data analysis

For the analyses presented in this paper, we needed representative diurnal cycles of O₂ and CO₂. We looked for a representative day in 2019 when little to no clouds were present; no unexpected behaviour in the diurnal cycles for potential temperature, specific humidity, or CO₂ occurred (e.g. caused by advection); and when the O₂ data showed a clear difference between the two measurement heights. We used data from 2019 instead of 2018 as 2018 saw a large-scale drought in Europe, and 2019 was less extreme and closer to a typical boreal summer (Peters et al., 2020). However, no single representative day could be found in our 2019 record, when the O₂ data showed a clear negative vertical gradient during the day and positive during the night, in combination with the above-mentioned meteorological criteria. We therefore choose a sequence of days to create an aggregate day based on the average of several days, which is representative for this time of the year in Hyytiälä, following the same method used by Ishidoya et al. (2015). The main criterion was that the vertical O₂ gradient had to be negative during the day, and the negative relationship between the change of O₂ and CO₂ concentrations over time at 23 m was present during the entire day. This resulted in selecting the period of 7 through 12 July 2019 to create the representative day, which we used in all subsequent analyses. The half-hourly values for the representative day are the averages of the data points of the individual half-hourly values for each day in the selected period. Each time step has an uncertainty that is based on the error propagation of the standard error (SE) of the 30 min averages for each day in the aggregate and is calculated for each time step with

$$\begin{array}{}\text{(4)}& {\mathrm{SE}}_{\mathrm{aggr}}={\displaystyle \frac{\sqrt{\sum {{\mathrm{SE}}_{\mathrm{day}}}^{\mathrm{2}}}}{n}},\end{array}$$

where n is the number of days included in the aggregate, SE_day is the standard error of the 30 min average of each individual day, and the SE_aggr is the resulting standard error of a 30 min value for the representative aggregate day.

For the representative day, the two O₂:CO₂ exchange ratio (ER) signals, ER_atmos and ER_forest, were determined. ER_atmos is based on O₂ and CO₂ concentrations and is expressed as

$$\begin{array}{}\text{(5)}& {\mathrm{ER}}_{\mathrm{atmos}}=-{\displaystyle \frac{{\mathrm{\Delta }}_{\left(t\right)}{\mathrm{O}}_{\mathrm{2}}}{{\mathrm{\Delta }}_{\left(t\right)}{\mathrm{CO}}_{\mathrm{2}}}},\end{array}$$

where both Δ_(t)O₂ and Δ_(t)CO₂ are the change in concentration over a selected time period (t). This is a unitless quantity as it represents mol O₂ per mol CO₂. ER_atmos was determined by the slope of a linear regression between the concentration of O₂ and CO₂ at the same height over a specific time period (Seibt et al., 2004; Stephens et al., 2007; Ishidoya et al., 2013; Battle et al., 2019). The selected time periods were based on the period when O₂ and CO₂ had the highest negative correlation. Throughout the day, this could be divided into three periods when different processes dominate (Fig. 1). It starts with the period during the night when the atmosphere is stable and when respiration becomes the dominant surface flux (P1); therefore, the CO₂ concentration increases and the O₂ concentration decreases. Subsequently, when the sun starts to rise, the boundary layer height starts to grow, and entrainment of air from the free troposphere influences the surface measurements (P2) (Vilà-Guerau de Arellano et al., 2004). Here the CO₂ concentration decreases rapidly and the O₂ concentration increases rapidly. Finally, the period starts when the effect of boundary layer dynamics and entrainment decreases, and the assimilation flux dominates (P3); here, the CO₂ concentration decreases less rapidly and the O₂ concentration increases less rapidly. We calculated an ER_atmos signal with Eq. (5) for the nighttime (P1), the daytime (by either focusing on only P3 or both P2 and P3), and the complete day (P1 + P2 + P3). The exact boundaries of these periods have to be estimated. To be certain about the exact times that should be taken as the boundaries for each period, an atmospheric model is needed.

ER_forest is based on O₂ and CO₂ fluxes and is expressed as

$$\begin{array}{}\text{(6)}& {\mathrm{ER}}_{\mathrm{forest}}=-{\displaystyle \frac{F_{{\mathrm{O}}_{\mathrm{2}}}}{F_{{\mathrm{CO}}_{\mathrm{2}}}}},\end{array}$$

where $F_{{\mathrm{O}}_{\mathrm{2}}}$ and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ are the net mean turbulent surface fluxes above the canopy of O₂ and CO₂ over a selected time period (Seibt et al., 2004; Ishidoya et al., 2015). We derive the fluxes of O₂ and CO₂ using the vertical gradient (see next paragraph and Eq. 7). The selected time periods for ER_forest were chosen such that the transition periods between the nighttime with a stable atmosphere (when the respiration flux dominates) and the daytime with a well-mixed atmosphere (when assimilation dominates) were excluded. By excluding the transition periods, we removed the periods when the gradients of both CO₂ and O₂ were close to zero. This was done because a very small gradient makes it difficult to calculate a flux and therefore the ER_forest and also because during this period entrainment is the most dominant process. The exact duration of the transition periods was based on the maximum and minimum of both the friction velocity and the height of 27 m (z) divided by the Monin–Obukov length (L). The friction velocity and ($z/L$) indicate the measure of turbulence of the atmosphere (Stull, 1988). The mean of the remaining data points of the CO₂ and O₂ flux during the stable atmosphere period was used to calculate the ER_forest signal of the night, and the mean of the remaining data points of the CO₂ and O₂ flux during the mixed atmosphere period was used to calculate the ER_forest signal of the day. The ER_forest for the entire day is taken as the average CO₂ and O₂ flux over the entire day. For this average, no periods are excluded, and all the data points over the 24 h are taken into account. Taking the average daily fluxes to derive ER_forest is a slightly different approach compared to the study by Ishidoya et al. (2015), who use the regression line between Δ_(z)O₂ and Δ_(z)CO₂ to determine ER_forest.

Currently, unlike for CO₂, the O₂ flux cannot be measured directly with an EC system. Instead, the flux can be inferred from the flux-gradient method. To calculate the flux of a certain scalar (ϕ) with the flux-gradient method, the following equation was used (Stull, 1988):

$$\begin{array}{}\text{(7)}& F_{\mathit{\varphi }}=-K_{\mathit{\varphi }}\cdot {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial z}},\end{array}$$

where F_ϕ is the surface flux of ϕ, K is the exchange coefficient, and $\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}/\partial z$ is the vertical gradient of $\stackrel{\mathrm{‾}}{\mathit{\varphi }}$. To determine the O₂ flux with Eq. (7) (where $\stackrel{\mathrm{‾}}{\mathit{\varphi }}=\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{2}}}$), the exchange coefficient of O₂ ($K_{{\mathrm{O}}_{\mathrm{2}}}$) needs to be determined. Ishidoya et al. (2015) assumed that $K_{{\mathrm{O}}_{\mathrm{2}}}=K_{{\mathrm{CO}}_{\mathrm{2}}}$ and determined $K_{{\mathrm{CO}}_{\mathrm{2}}}$ by dividing the CO₂ flux, measured with EC, by the CO₂ vertical gradient between two measurement levels. However, the exchange coefficient can also be determined with other methods that, for example, only need two measurement heights for the vertical gradient. In this study, we explore these different options for calculating $K_{{\mathrm{O}}_{\mathrm{2}}}$. The EC measurements of the CO₂ flux were used as a reference to determine the most suitable approach. The most suitable approach to infer the O₂ flux is then used for both $K_{{\mathrm{CO}}_{\mathrm{2}}}$ and $K_{{\mathrm{O}}_{\mathrm{2}}}$. During this study, we derive the surface flux in the surface layer (Fig. 1), and we assume that the surface flux stays constant in this surface layer, which consists of the roughness sublayer and the inertial sublayer.

We categorized the methods to determine the most suitable K into two groups: the observation-based approach (also called the K-theory (Stull, 1988) or the modified Bowen ratio method (Meyers et al., 1996)) and the theoretical approach (following the similarity theory (Dyer, 1974)). For the observation-based methods, we determined the exchange coefficient (K) in Eq. (7) by dividing a flux measured at 27 m, using an EC system, by a three-height (16, 67, and 125 m) vertical gradient of a specific scalar. Ishidoya et al. (2015) used this approach to calculate their O₂ flux, using the CO₂ flux and vertical gradient of two levels. Next to CO₂, we also calculated K using potential temperature (θ) for the observation-based approach.

For the theoretical approach, the K in Eq. (7) is determined with the Monin–Obukov similarity theory (MOST) (Dyer, 1974), where logarithmic surface layer scaling applies for K, and empirical similarity functions are used to describe the effect of atmospheric stability. In addition, we used a correction which takes into account the effect of the roughness sublayer (see Appendix A for details). The SMEAR II data at 27 m were used for the calculations with MOST. When only two heights for the gradient calculations are available, there is an option to integrate Eq. (7) (de Ridder, 2010). We tested both the application with and without integration in this study. We used the ICOS data, available at the SMEAR II station, for the K calculations. For the CO₂ EC measurements, we used the gap-filled data to correct for the storage below the measurement height of the EC. Gap filling was applied when the friction velocity (u*) was below 0.4 m s⁻¹ (Kulmala et al., 2019). Appendix A gives a more elaborate explanation and provides equations of the different methods used to determine the exchange coefficients used in this study.

Finally, we select the K_ϕ from either the observation-based or the theoretical approach that produced CO₂ flux results from our CO₂ vertical gradient measurements that showed the best comparison to the EC CO₂ flux measurements. This K was used to calculate the O₂ and CO₂ fluxes, together with the vertical gradient from measurements collected during our campaigns. For our campaigns, we only have O₂ and CO₂ measurements at two heights (23 m and 125 m), which means that $\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}/\partial z$ changes into $\mathrm{\Delta }\stackrel{\mathrm{‾}}{\mathit{\varphi }}/\mathrm{\Delta }z$, and the gradient was calculated with finite differences.

After both the CO₂ and O₂ fluxes were determined, resulting in ER_forest, we subsequently calculated the O₂ : CO₂ exchange ratio signals for the assimilation processes (ER_a) and the respiration of the ecosystem (ER_r) with the following equations (Seibt et al., 2004; Ishidoya et al., 2015):

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }\mathrm{NEE}=-\mathrm{GPP}+\mathrm{TER},\text{(9)}& {\displaystyle }& {\displaystyle }\mathrm{NEE}\cdot {\mathrm{ER}}_{\mathrm{forest}}=-\mathrm{GPP}\cdot {\mathrm{ER}}_{\mathrm{a}}+\mathrm{TER}\cdot {\mathrm{ER}}_{\mathrm{r}},\end{array}$$

where NEE is the net ecosystem exchange, GPP is the gross primary production, and TER is the total ecosystem respiration. GPP and TER are always positive by definition, representing uptake and release by the ecosystem, respectively. Therefore, when GPP is larger than TER the resulting negative NEE values represent carbon uptake by the ecosystem. First, we assumed that nighttime NEE is equal to TER, which meant that the nighttime ER_forest signal is equal to ER_r. We assumed that the processes that contributed to the ER_r keep the same ratio between O₂ and CO₂ during the entire day and therefore we used a constant ER_r for the entire day. We base this assumption on studies that showed that the variability of ER_r highly depends on the bulk soil respiration (Hilman et al., 2022; Angert et al., 2015). No large changes occur in the soil temperature or the soil moisture during our (representative) diurnal cycle; therefore, we assume that the ER_r of the bulk soil respiration stays relatively constant, and with that the ER_r of the ecosystem also stays constant over the entire day. Subsequently, we calculated ER_a for both the entire diurnal cycle and the daytime using Eq. (9) with the corresponding ER_forest and the constant ER_r. We used ICOS NEE EC measurements from the SMEAR II station at a level of 27 m in the 128 m high tower. The GPP fluxes at Hyytiälä are calculated with either of the following two approaches: (1) when NEE EC measurements are available, GPP is calculated as the difference between the NEE EC measurements and the respiration flux, which is calculated using a temperature function; or (2) when NEE EC measurements are not available, GPP is calculated using an equation that is based on the air temperature and light (photosynthetically active radiation, PAR). A more detailed description of these calculations is given by Kulmala et al. (2019) and Kohonen et al. (2022).

Figure 3 The half-hourly average O₂ (a) and CO₂ (b) concentrations at Hyytiälä for spring/summer of 2018 and 2019 for the 125 and 23 m height levels, together with the vertical gradient (Δ_(z)) between these two heights (c) for both O₂ and CO₂. The shaded areas indicate the dates that were selected for the aggregate representative day (7 through 12 July 2019: yellow) and the second representative day to test the O₂ method (13 through 15 June 2019: grey). The selected days for the aggregate representative days are shown in more detail for the 23 m measurements for the gradients during 13 through 15 June (d) and (e) and during 7 through 12 July (f) and (g) for both O₂ and CO₂.

By estimating ER_r and ER_a of this boreal forest, we created the opportunity to apply atmospheric O₂ measurements to separate NEE into GPP and TER (the O₂ method). We calculated ER_r and ER_a for the representative day using Eqs. (8) and (9), and we use these to calculate GPP and TER for another representative day. We selected 13 through 15 June to create a new second aggregate day and to calculate a new ER_forest signal for the entire day (see Fig. 3d and e for a detailed view on the measurements of those days). These three days were chosen because in 2019 they showed the clearest diurnal cycle of O₂ and a negative O₂ gradient, aside from 7 through 12 July, used above. We assume here that the ER_r and ER_a calculated for the period from 7 through 12 July are representative for the period from 13 through 15 June. Studies show that the ER_r (Hilman et al., 2022) and ER_a (Bloom, 2015; Fischer et al., 2015) values vary with changing soil and atmospheric conditions. The periods for both representative days are relatively close in time and therefore have similar conditions in the soil and the atmosphere, and we can therefore assume that the ER_r and ER_a values based on the 7 through 12 July data can also be applied to the 13 through 15 June period. By using the ER_r and ER_a values determined for the first representative day (7–12 July) and ER_forest and NEE for the second representative day (13–15 June), we calculated GPP and TER from NEE for this second representative day. By comparing the GPP and TER fluxes of the O₂ method to the GPP and TER fluxes of the temperature-based function of ICOS (EC method), we could demonstrate how accurate the O₂ method is (Sect. 3.4). Both Seibt et al. (2004) and Ishidoya et al. (2015) also applied the O₂ method; however, both of these studies used chamber measurements to first determine ER_a and ER_r and then used Eqs. (8) and (9) to infer GPP and TER. Unfortunately, we did not have chamber measurements of both O₂ and CO₂ available and therefore we used Eqs. (8) and (9) to calculate ER_a and ER_r. This means that these two equations can be used in two ways: to determine the ER_a and ER_r signal or to separate NEE into GPP and TER.

The footprint of the calculated O₂ and CO₂ surface fluxes that also represents the footprint of the ER_forest, ER_r, and ER_a signals for the representative aggregate day is shown in Fig. B2 in Appendix B. The footprint is based on the method by Kljun et al. (2015), where for the height the geometric mean between 125 and 23 m is used. The footprint analysis shows that the surface fluxes are mainly influenced by the forest surrounding the tower and that the lake located close to the tower is not influencing the signal. The footprint of the O₂ and CO₂ concentrations and therefore the footprint of the ER_atmos signal can be found in the document by Carbon Portal ICOS RI (2022). This concentration footprint analysis shows that with an average wind direction of north to northeast during 7 through 12 July 2019, the concentrations measured are mainly originated from forest exchange, with hardly any influence of urban sources.

Table B1 in Appendix B gives a complete overview of which data are used for each part of this research for the two different aggregate days.

3 Results

3.1 O₂ and CO₂ time series

The calibrated half-hourly measurements of O₂ and CO₂ for 2018 and 2019 are shown in Fig. 3, together with the vertical gradients between the two measurement heights. The O₂ measurements are shown here converted from per meg to ppmEq, which is to allow for comparison of the diurnal variability for CO₂ and to calculate the ER signals. The differences between the 23 and 125 m measurements are observable for both CO₂ and O₂. During both campaigns in 2018 and 2019, the diurnal behaviour of the O₂ concentrations are anticorrelated with the CO₂ concentrations. This anticorrelation between O₂ and CO₂ is also visible from the gradient measurements, despite the relatively high uncertainty of the O₂ measurements as described in Sect. 2.2 and further elaborated on in Sect. 4.1. The period 7 through 12 July 2019 shows the most clear negative relationship between the O₂ gradient and the CO₂ gradient, and it also had the most suitable meteorological conditions and was therefore selected for the aggregate representative day (Sect. 2.3). The period 13 through 15 June shows a less clear anticorrelation between the vertical gradients of O₂ and CO₂ (Fig. 3d and e) but with clear diurnal cycles of O₂ and CO₂ suitable for the purpose of our second aggregate day (see Sect. 3.4).

3.2 Diurnal cycles

Figure 4 Diurnal cycles in local wintertime (LT; time zone UTC+2) (all times in this paper are given in local time unless stated otherwise) of the O₂ and CO₂ concentrations for the 23 and 125 m height levels (a) and the vertical gradient between both levels with the uncertainty of both O₂ and CO₂ of the representative day, taken as the average values of 7 through 12 July 2019 (b). The CO₂ measurements of the ICOS setup are shown in (a) for comparison with the CO₂ setup measured during our campaigns. The shaded colours indicate the selected different periods when the most dominant processes are the following: stable atmosphere and respiration (00:00–04:00, P1); entrainment, boundary layer growth, and assimilation (04:00–09:00, P2); convective conditions and assimilation (09:00–13:00, P3a); and the same as P3a plus a remaining artefact for the O₂ measurements after the pressure correction as explained in the text (13:00–20:00, P3b). The vertical dotted lines indicate the sunrise (03:00) and sunset (22:00). The error bars in panel (b) are half-hourly standard errors based on the error propagation of the standard errors of the data points in (a) (not shown), which were based on Eq. (4).

The measurements of O₂ and CO₂ and their vertical gradient for the representative day, are shown in Fig. 4. There are no measurements between 20:00 and 22:00, because the calibration cylinders were measured during this period. For 7 through 12 July, we used a fixed calibration time, as radiosondes were launched (not shown) during this period, and we wanted to make sure we captured the morning transition to compare with these radiosondes. Note that the daylight length at Hyytiälä is long at this time of the year, with sunrise at 03:00 and sunset at 22:00. We compared our CO₂ observations with ICOS CO₂ measurements at the same height, which shows that both instruments compare well overall, with a mean difference of 0.70 ± 0.65 ppm during the period 7 through 12 July. The comparison between the two devices was a bit difficult because of the different timing of the measurements. The diurnal cycles of O₂ and CO₂ (Fig. 4a) clearly show anticorrelated behaviour between CO₂ and O₂, which is especially visible during nighttime (23:00–04:00) and the morning transition (05:00–13:00).

Figure 4 shows four different periods that can be linked to the periods to calculate ER_atmos, described in Sect. 2.3. P1 is visible between 23:00–04:00, where respiration starts to dominate the signal and therefore the O₂ concentration decreases and the CO₂ concentration increases, in a decreasing boundary layer height dominated by thermal stratification. P2 becomes visible around 04:00 and stops around 09:00, where entrainment, the growing boundary layer, and the onset of photosynthesis causes a steep increase in the O₂ concentration and a steep decrease in the CO₂ concentration. P3 can be divided into P3a and P3b and is visible between 09:00–20:00. Between 09:00–13:00 (P3a), the photosynthesis flux starts to dominate, and both the O₂ and CO₂ concentration increase and decrease less rapidly. Between 13:00–20:00 (P3b), the O₂ concentration starts to decrease while the assimilation flux still dominates, which is a remaining artefact from the pressure correction that we applied due to the instability of the MKS pressure transducer (see Sect. 2.2). As shown in Fig. B1 in Appendix B, higher daytime temperatures cause larger PMKS deviation and therefore the effect of the pressure correction is largest during midday, leading to a larger uncertainty in the observations in that time period. The boundary of 20:00 between P3b and P1 was difficult to determine as we missed some measurements due to the calibration period, and the remaining measurements around this time have a deviation caused by the pressure transducer. Measurements at both levels show this same diurnal behaviour; however, it is more pronounced closer to the vegetation (the 23 m level).

Figure 5 The O₂ concentration plotted against the CO₂ concentration for the representative day in local wintertime (LT; time zone UTC+2), with the 23 m level in coloured points per period representing different dominant process and with the 125 m level in grey points. The dominant processes are the following: respiration (00:00–04:00), entrainment (04:00–09:00), assimilation (09:00–13:00), and a remaining artefact after the pressure correction due to the instability of the MKS pressure transducer becomes visible (13:00–20:00). The linear regression lines indicate the exchange ratio of the atmosphere (ER_atmos) during the time with a specific dominant process.

The difference between the two heights results in a vertical gradient (Fig. 4b). Similar to the diurnal cycle of the concentrations, the diurnal cycles of the gradients of O₂ and CO₂ also show anticorrelated behaviour. At 08:00, the CO₂ gradient changes from negative to positive, and the O₂ gradient changes from positive to mostly negative, reflecting CO₂ being transported downwards and O₂ upwards, respectively. The magnitude of the gradient depends on the degree of vertical mixing. The sign of the gradients changes during the day, because the lowest level (23 m) is more directly influenced by forest carbon exchange compared to the highest level (125 m). Around the time of sunset, the CO₂ gradient changes from positive to negative, and the O₂ gradient changes from negative to positive, because the lowest measurement level (23 m) is now influenced more by respiration processes of the forest and soils compared to the highest measurement level (125 m). The error bars are based on the error propagation of the standard errors of each half-hourly data point that were calculated with Eq. (4). The gradient of O₂ is hardly affected by the PMKS correction (see Fig. B1), as measurements at both heights are affected similarly.

By using Eq. (5), we calculated four distinct ER_atmos signals for different periods throughout the day at 23 m, as well as to a smaller degree at 125 m (Fig. 5 and Table 3). The same periods as shown in Fig. 4 are visible in Fig. 5. This results in an ER_atmos during the night (P1) of 1.22 ± 0.02 and two different possibilities for the ER_atmos signal during daytime. By combining both P2 and P3a, we get a signal of 2.28 ± 0.01, and by focusing only on P3a, which excludes the entrainment and the boundary layer dynamics, we get a signal of 1.10 ± 0.12. Last, by combining all the periods (P1, P2, P3), we get a signal for the complete day of 2.05 ± 0.03. The uncertainties given here only represent the uncertainty of the slopes from the regression lines in Fig. 5. The high values for the ER_atmos signal of the entire day and the daytime signal that includes entrainment and the boundary layer dynamics are not very realistic to represent an ER for the forest, and this shows that we should be careful when using ER_atmos. This will be elaborated on in Sect. 4.2.

3.3 Flux calculations for CO₂ and O₂

We explored four alternative methods to derive the O₂ flux from the vertical gradient of the two measurement levels, as described in Sect. 2.3. Figure 6 shows both the theoretical and the observation-based approaches that were used to calculate the CO₂ flux, in comparison to the ICOS EC CO₂ flux measurements at 27 m on the tower. By comparing these approaches to the EC measurements, we determined which method is most suitable to calculate the O₂ flux. The CO₂ flux measured by the EC system stays positive until around 05:00, when the respiration fluxes are the most dominant and the nocturnal boundary layer is shallower. After 05:00, the CO₂ flux of the EC system becomes negative, and the forest begins to take up CO₂ instead of emitting it. The assimilation fluxes increase and exceed the respiration fluxes, the boundary layer starts to grow, and air with lower CO₂ concentrations is entrained from the free troposphere. After 20:00, the CO₂ flux of the EC system becomes positive again as the assimilation fluxes decrease, and the respiration signal begins to dominate again while the boundary layer height decreases. We expect to find this diurnal pattern and the sign change in our calculations of the CO₂ flux from the vertical gradient method as well.

Figure 6 The CO₂ flux (a) calculated with different methods for the representative day, as described in Sect. 2.3, compared to the CO₂ flux of the ICOS EC measurements. (b) The comparison between the O₂ and CO₂ flux calculated using the method that gave the best results for the CO₂ flux calculations (using the exchange coefficient K with CO₂) for the representative day. The shaded colours indicate the regions that were selected for the following: the night signal (21:00–04:00), the day signal (09:00–17:00), and the remaining regions (04:00–09:00 and 17:00–21:00), with the time in local wintertime (LT; time zone UTC+2). The error bars of (b) are based on the error propagation of the standard error of the 30 min values for the representative day, which are based on Eq. (4).

First, we discuss the theoretical methods that are indicated in Fig. 6 with “K with MOST” and “Integrated” approaches (see Sect. 2.3). The MOST and the integrated methods both overestimate the CO₂ flux during the night, between 0:00 and 05:00. In addition the resulting CO₂ flux decreases and becomes negative too late in the day compared to the EC measurements. The CO₂ fluxes of the MOST and integrated methods evolve from a positive flux to a negative flux around 8:00. This is 3 h later than the CO₂ flux from the EC measurements. During the day, between 08:00 and 15:00, the MOST method underestimates the CO₂ uptake and the integrated method overestimates it. Table 2 shows that both MOST and the integrated methods have the highest mean difference and root mean square error (RMSE) compared to the observation-based approaches. We discuss this further in Sect. 4.3. As a result of this analysis, we decided to not use the theoretical approach to calculate the O₂ flux.

Table 2 The mean difference and the root mean square error (RMSE) of the comparison between the EC CO₂ flux measurements at 27 m in the tower and the CO₂ flux calculated with different methods for the exchange coefficient K, based on the ICOS data, each using the vertical gradient of CO₂ at 23 m and 125 m of our campaign data.

Secondly, we analyse the observation-based approaches that are indicated in Fig. 6 with “K with θ” (where K is established using ICOS vertical gradients of potential temperature and the sensible heat flux) and “K with CO₂” (where K is established using ICOS CO₂ vertical gradients and CO₂ EC data). The observation-based approaches showed a better comparison with the EC observations in determining the CO₂ flux compared the to theoretical approach. Both the θ and the CO₂ methods represent satisfactorily the nocturnal CO₂ flux between 00:00 and 05:00. After 05:00, the fluxes calculated by both methods start to decrease and change sign around the correct time (05:00) from a positive to a negative flux. During the day between 08:00 and 15:00, both the θ and the CO₂ methods underestimate the CO₂ flux but not as much as the theoretical methods. Table 2 also shows that both the θ and the CO₂ methods have the lowest mean difference and RMSE. Based on the smaller mean difference and RMSE, as well as the direct link of CO₂ with O₂, we decided to proceed with the method where K is calculated with the ICOS data of CO₂ instead of the ICOS θ data. This K was then multiplied with our measured O₂ vertical gradient between 23 and 125 m to calculate the O₂ flux. Section 4.3 presents a more complete discussion on the different methods to determine the most suitable K.

The resulting O₂ flux calculated with the exchange coefficient K based on the ICOS CO₂ data is shown in Fig. 6b. The uncertainties are based on the error propagation of the standard errors of the O₂ and CO₂ data per time step as calculated with Eq. (4), in Eq. (7). We do not calculate an uncertainty for K, as this is not the dominating term contributing to the total uncertainty. The daytime flux values have a high variability, but the inferred fluxes appear physically realistic and promising for one of the first attempts to calculate O₂ fluxes. During the night, between 0:00 and 5:00, the O₂ flux signal has a relatively stable negative value, because the forest consumes O₂ for the respiration processes. Similarly, CO₂ is released during the night, leading to a positive CO₂ flux. After 5:00, the O₂ flux becomes positive and shows a higher variability. Overall, the O₂ flux is positive during the day, which indicates that the forest produces O₂ as the assimilation rate is higher than the respiration rate. The high variability of the O₂ flux compared to the CO₂ flux is caused by the less precise measurements of the O₂ vertical gradient compared to the CO₂ gradient (Fig. 4). The measurement precision needed to measure the difference between the two levels is very high and therefore impacts the measurement of the gradient of O₂. The nighttime values of the O₂ flux are therefore more reliable than the daytime values, since the difference between the two heights is larger and due to the more stable atmospheric conditions at night.

Figure 7 The half-hourly exchange ratio of the forest (ER_forest) and the resulting averaged ER_forest for the entire day (black line), the night between 21:00–04:00 (dark blue line), and the day between 09:00–17:00 (light blue line) of the representative day (a) with the time in local wintertime (LT; time zone UTC+2). The size of the dots indicates the size of the absolute O₂ flux, and the shaded bands indicate the uncertainties of the different ER_forest signals. Note that the ER_forest lines do not match with the average of the dots in the specific time period, because the lines are based on the averaged fluxes. These different ER signals are presented in a vector diagram format with the carbon fluxes, gross primary production (GPP), total ecosystem respiration (TER), net ecosystem exchange (NEE), the ER of the assimilation processes (ER_a), and the ER of the respiration processes (ER_r) (b).

By using Eq. (6), we find three different ER_forest signals throughout the day (Fig. 7 and Table 3). The selected time periods based on the criteria described in Sect. 2.3 are between 09:00–17:00 for the daytime and between 21:00–04:00 for the nighttime (Fig. 6). This results in a nighttime ER_forest signal of 1.04 ± 0.04, a daytime ER_forest signal of 0.92 ± 0.17, and an ER_forest signal for the entire 24 h of 0.83 ± 0.24. Note that this 24 h value is not the average of the day and night ER_forest signals or from all the 30 min ER_forest signals, because we used the averaged fluxes. This means that the ER_forest signals based on high flux values, indicated in Fig. 7 with larger symbols, contribute more to the averaged ER_forest signals compared to the lower flux values. Figure 7b illustrates that when combining surface fluxes with different sign, we cannot just average the corresponding ER signals (see Sect. 4.4). The individual ER_forest values of every 30 min show a clear difference between the daytime and nighttime. The ER_forest values during the nighttime are relatively stable. The ER_forest values during the daytime show more variability, caused by the high variability of the O₂ flux during daytime (Fig. 6). The uncertainty of the ER_forest signals is determined by the propagation of the standard error of the aggregate 30 min data (based on Eq. 4), in Eqs. (7) and (6).

Table 3 The exchange ratio for the atmosphere (ER_atmos: Sect. 3.2), the forest (ER_forest: Sect. 3.3), and assimilation and respiration (ER_a and ER_r: Sect. 3.3) for different time periods of the representative day. The time periods used to calculate the signals are the following: 09:00–13:00 for day and 23:00–04:00 for night of ER_atmos and 09:00–17:00 for day and 21:00–04:00 for night of ER_forest, ER_r, and ER_a. Note that the uncertainty for ER_atmos does not represent the same uncertainty as for ER_forest, since the first is the error of the fit, and the second is based on error propagation of the half-hourly measurements.

3.4 GPP and TER calculations

We found the ER signals for assimilation (ER_a) and respiration (ER_r) by using Eq. (9) (Fig. 7b and Table 3). The assumption that ER_r stays constant throughout the day seems reasonable, because the ER_forest values stay stable during the night. The ER_r signal therefore becomes 1.03 ± 0.05. A more elaborate discussion of this assumption can be found in Sect. 4.5. ER_a of the daytime is 0.96 ± 0.11, which indicates the ER_a signal of the boreal forest when the surface fluxes are the highest. The ER_a signal of the entire diurnal cycle is 0.95 ± 0.11, which also includes the assimilation processes during sunrise and sunset. Figure 7b shows all these ER signals and how they change throughout the day, together with their carbon fluxes. ER_a, ER_r, and the resulting ER_forest signals are more realistic compared to the ER_atmos signals. We will elaborate on these differences in Sects. 4.4 and 4.5.

By using Eqs. (8) and (9) for a second representative day (13 through 15 June), with the ER_a and ER_r signals determined from the representative day, we show in Fig. 8 that the O₂ method compares well with the EC method. This means that the O₂ method could potentially be used to separate NEE into GPP and TER on any day when good simultaneous CO₂, O₂, and NEE measurements are available. The difference between the CO₂ fluxes determined with the O₂ method and the EC method of both the GPP and the TER flux are around 0.5 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, which is less that 6 % of the total gross flux. The difference is relatively small, which means that the O₂ method compares well with the EC methods to separate NEE into GPP and TER. The different error bars in Fig. 8 show how sensitive the O₂ method is to the accuracy of ER_forest. By changing ER_forest by 0.2, the GPP estimation by the O₂ method changes by 4 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, and by changing ER_forest by only 0.01, the GPP estimation changes by 0.2 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. The effect of changing ER_forest on TER has the same effect on GPP. This shows that the O₂ method is quite sensitive to ER_forest and should be measured accurately, with a suggested precision of around 0.05. With a precision of 0.05 for ER_forest, the GPP and TER fluxes derived with the O₂ method stay in the same range as the GPP and TER fluxes determined with the EC method. The application of the O₂ method will be further discussed in Sect. 4.5.

Figure 8 The CO₂ fluxes of a second representative day (13 through 15 June) for net ecosystem exchange (NEE), gross primary production (GPP), and total ecosystem exchange (TER) based on two different methods: the EC method and the O₂ method. The different error bars indicate an increase/decrease of 0.2, 0.1, or 0.01 for the exchange ratio of the forest (ER_forest) used in the O₂ method.

4 Discussion

We aimed to advance understanding of the O₂:CO₂ exchange ratio and its diurnal variability over a boreal forest by continuously measuring both O₂ and CO₂ concentrations at two heights above the canopy. These measurements gave us the possibility to compare the ER_atmos and ER_forest signal of an aggregate representative day and compare the boreal forest signals to previous studies in different ecosystems. Our ER_atmos signal changed between the day (2.28) and the night (1.22) and had an overall diurnal signal of 2.05. For the ER_forest signal, we needed to determine the O₂ and CO₂ surface fluxes based on the two heights. Different flux-calculating methods were compared. The O₂ flux was calculated with the method that resulted in the best comparison to EC fluxes for CO₂, where we found that the exchange coefficient K based on the CO₂ data was most suited. The resulting ER_forest signal showed again differences between the day (0.92) and night (1.04), and the overall diurnal ER_forest was 0.83. For these differences and variability in the ER signals, different aspects of the uncertainty have to be taken into account, on which we elaborate in the next sections.

4.1 Measurement uncertainty

By analysing the mean difference and standard deviation of the target cylinder values between 16 June 2019 and 17 July 2019 (Table 1), we see that the values are relatively high. Previous studies that used a fuel cell analyser for continuous atmospheric O₂ measurements (Battle et al., 2019; Ishidoya et al., 2013; van der Laan-Luijkx et al., 2010; Popa et al., 2010; Pickers et al., 2022) achieved measurement precision of around 5 per meg. The WMO recommends a compatibility goal of 2 per meg; however, this is difficult to achieve and so the extended compatibility goal is 10 per meg for the worldwide O₂ monitoring network (Crotwell et al., 2020), which shows that our long-term measurement precision of 19 per meg is relatively poor. This poor measurement precision could have been caused by several reasons. The O₂ values of the calibration cylinders that were used were relatively far apart, making it more difficult to measure the values around the target cylinder value. For 2018, we used calibration cylinders with the following values (on the SIO scale): −628.53, −816.17, and −1208.28 per meg, and for 2019 we used cylinders with values of −729.96, −816.17, and −1208.28 per meg. The cabin in which the instrument and cylinders were located was not well insulated, which created unstable temperature conditions that might have affected the stability of the cylinders (Keeling et al., 2007). The calibrations of our representative aggregate day took place during the night; therefore, large temperature changes during the day might have affected daytime stability of the reference cylinder. Furthermore, tiny leakages in the setup might have influenced the measurements. Due to the relatively short period for these campaigns and the remote location, it is not possible to trace back the cause of this large uncertainty. This high uncertainty resulted in a larger uncertainty of the vertical gradient of the two heights of the O₂ measurements. However, in this study we are mostly interested in the diurnal variability of the ER signal and differences between ER_atmos and ER_forest; therefore, the long-term stability of the measurements are less relevant here compared to other O₂ studies.

To reduce the effect of the high measurement uncertainty and derive a more statistically robust signal of the vertical gradient, we created an aggregate representative day based on days with similar weather and atmospheric conditions. The increased statistics of this representative aggregate day decrease the effect of the low measurement precision. We also move away from the reality of one specific day but rather focus on an average situation and variability of the ER signal above a boreal forest based on O₂ and CO₂ measurements at two levels. Given that very few previous studies focused on deriving forest ER signals globally, our analysis helps to gain further understanding of the diurnal variability and the difference between ER_atmos and ER_forest, which will be discussed in the following sections.

4.2 ER_atmos signal in comparison to previous studies

Despite the uncertainty in our measurements, there are clear differences between the slopes of O₂ and CO₂ throughout the diurnal cycle (Fig. 5). Three different ER_atmos signals are visible, with two signals for the day (2.28 ± 0.01 and 1.10 ± 0.12) and one for the night (1.22 ± 0.02) slope (Table 3). Note that the uncertainty of these values is based on the slope of the fitted line in Fig. 5 and does not represent the uncertainty in the stability of our measurements indicated in Table 1. The difference between day and night values of ER_atmos was expected, because different processes (i.e. respiration, assimilation and entrainment) with different ER signals play a role at different times during the diurnal cycle. To exclude as much as possible the effect of entrainment and the boundary layer dynamics during the morning transition, we will from now on refer to the 1.10 value as the day ER_atmos signal, which is the signal derived form period P3a. ER_atmos for the complete day results in 2.05 ± 0.03.

Table 4 The different exchange ratio (ER) signals of previous studies are given, with the ER of the atmosphere (ER_atmos), the ER of the forest (ER_forest), the ER of the respiration processes (ER_r), and the ER of the assimilation processes (ER_a). Bat, 2019 is short for Battle et al. (2019), Ish, 2015, represents Ishidoya et al. (2015), Ish, 2013, represents Ishidoya et al. (2013), and Sei, 2004, represents Seibt et al. (2004).

^a An ER signal is classified as ER_atmos when the ER signal is based on one concentration measurement of O₂ and CO₂. ^b An ER signal is classified as ER_forest when the ER signal is based on surface fluxes from either a one-box model or vertical gradient flux calculations. ^c The ER signals of the location Griffin Forest by Seibt et al. (2004) are used here. ^d The ER signals of the location Harvard Forest by Seibt et al. (2004) are used here.

When comparing our ER_atmos signals to those from Battle et al. (2019), Ishidoya et al. (2013), and Seibt et al. (2004) (Table 4), we note several similarities but also some differences regarding the specific values of the ER_atmos signals. Our daytime signal of 1.10 is similar to 1.02, 0.87, and 1.14 from the previous studies, respectively, as is our nighttime signal of 1.22 compared to 1.12 (Battle et al., 2019), 1.03 (Ishidoya et al., 2013), and 1.16 (Seibt et al., 2004). However, our 24 h ER_atmos signal of 2.05 shows an unrealistically high number which clearly does not indicate the ER of the forest only. A typical ER_atmos signal for a 24 h period lies around 1, as is shown in Table 4 and by Stephens et al. (2007) and Manning (2001). Our 24 h ER_atmos value includes the measurement points of the period that is influenced by entrainment and boundary layer dynamics (P2), for which period we found an ER signal of 2.28. The large influence of entrainment and boundary layer dynamics made it difficult to be very precise about the specific time periods to choose for P3. Moving the selected time boundaries of P3a from 09:00 to 09:30 or from 13:00 to 12:30 leads to ER_atmos values of 0.88 or 1.75, respectively. The large changes in the daytime ER_atmos due to small changes in the time boundaries show the high uncertainty of the daytime ER_atmos. Therefore, our measurements provide a confirmation of earlier indications (Seibt et al., 2004) that ER_atmos is an unreliable estimate for the ER of a forest, and we recommend to use ER_forest.

Instead, ER_atmos also represents how O₂ and CO₂ are influenced by the boundary layer dynamics and entrainment (Fig. 1). The high ER_atmos values cannot be explained by signals from other sources, such as fossil fuel combustion or exchange with the lake, as both are not represented in the footprint of our measurements (see Sect. 2.3). Furthermore, we have shown that these high values are not an artefact from the instability of the pressure stabilization, as preliminary analysis of the ER_atmos values from our 2018 measurements also show values higher than 2.0 (not shown). Although we cannot fully rule out remaining artefacts in the calibration due to, for example, temperature changes in the measurement cabin, we suggest that the more plausible explanation is that ER_atmos is highly influenced by atmospheric processes, such as entrainment. The entrainment of air from either the residual layer (early in the morning transition) or the free troposphere (after the residual layer is dissolved) could impact the ER_atmos as different sources of air are mixed. The residual layer contains air from the day before and could be affected by horizontal advection, whereas the air in the free troposphere originates from different background sources. These difference sources can have different ER signals and therefore create a mixture of air where O₂ and CO₂ are influenced differently. These air masses affect O₂ differently compared to CO₂ in the boundary layer, and an ER_atmos signal will arise that cannot be linked directly to one specific process. Even though entrainment processes also occur at locations of previous studies, we still find differences in ER_atmos. We suggest that this can be explained by difference in measurement height compared to the canopy height and different sources of background air in the free troposphere at the measurement location. For the ER_atmos signal during P2 at 125 m, we find a value of 3.40, even higher than the ER_atmos signal of 2.28 at 23 m, which indicates that the influence of entrainment increases when measuring further away from the canopy, and as a result the ER_atmos signals show higher values. Further insights into the contributions of each process to ER_atmos cannot be estimated from the measurements alone and would require using an atmospheric model.

4.3 Uncertainties in the CO₂ and O₂ flux calculations

By comparing the theoretical and observation-based methods, we determined that the most suitable method to calculate both the CO₂ and O₂ fluxes was to use the observation-based method with CO₂ data (Sect. 3.3). Figure 6 and Table 2 show that the theoretical methods (MOST and integrated) resulted in a change of the CO₂ flux that was late compared to the EC measurement. This delay has been described before and is caused by the time it takes before the turbulence can mix the CO₂ gradient driven by stable nocturnal stratification conditions and establish the corresponding gradient to how turbulent the atmosphere is (Casso-Torralba et al., 2008). When the heights of the gradient are closer together, the delay is less pronounced. However, the measurement heights used during our campaign are relatively far apart (125 m and 23 m), and the EC flux is measured at 27 m. The 125 m measurement is even located outside the surface layer during the morning transition. This made the flux-gradient method (as described in Eq. 7) less applicable, which assumes that the surface flux stays constant in the surface layer (Dyer, 1974).

Since during our campaign we only measured at two heights, we missed information on the logarithmic profile originating from the canopy top, which resulted in an underestimation of the flux using the K with MOST method. This was solved by integrating the MOST equation (integrated method). With the integrated method, the gradient is assumed to be logarithmic, and the total flux increases compared to the MOST calculation (Paulson, 1970). However, with the large difference between the two measurement heights, the integrated approach still overestimated the CO₂ flux compared to the EC measurements during both the day and the night. Also, the delay in the timing of the sign change of the gradient cannot be solved with this integrated method. We also explored the effect of adding a roughness surface layer (RSL) in the flux calculations of the theoretical methods, by adding an extra factor that accounts for this layer (not shown in the results) (de Ridder, 2010). The contribution of the RSL did not improve our results, because it also includes the delay of the gradient which was causing the largest deviation in the theoretical methods (Table 2).

By applying both observation-based methods, using either θ or CO₂ to infer the exchange coefficient K, we did not find this delay in the timing of the gradient, and the observation-based methods therefore resulted in derived fluxes close to the EC measurements. Here it has to be noted that the ICOS EC measurements of CO₂ that we used as a benchmark for the most suitable flux calculation approach were also used in calculating K with CO₂, which makes the comparison of these approaches to the CO₂ flux not fully independent. Note that we first derive K with the vertical CO₂ gradients calculated from ICOS CO₂ observations at three vertical levels, and we apply this to our own measurements of the CO₂ vertical gradient with an independent instrument (Table B1). As a result, there is not a full circularity when comparing the obtained fluxes to the EC CO₂ measurements to select which method for calculating K we use. Most previous studies that determined fluxes based on the gradient approach used θ to calculate K (Stull, 1988; Mayer et al., 2011; Wolf et al., 2008; Bolinius et al., 2016; Brown et al., 2020), because θ is the driver of convective turbulence. However, because O₂ is directly linked to CO₂ and because our statistics (Table 2) indicated that the CO₂ method resulted in a better comparison to the EC fluxes, we decided to use the ICOS CO₂ data at three levels and the CO₂ EC measurements to calculate K. This K together with the measurements of two heights by our instrument during our campaign were used to calculate both the CO₂ and the O₂ fluxes used in our study. We also tested the impact of using only two vertical levels of the ICOS CO₂ concentrations to calculate K (not shown), which was also the case in the only previous study that derived O₂ fluxes. Ishidoya et al. (2015) derived O₂ fluxes for a temperate forest in Japan using two vertical levels at 18 and 27 m height for both O₂ and CO₂ concentrations. Our comparison of deriving K based on two vertical levels (23 m and 125 m) resulted in an underestimation of the gradient and thus an overestimation of K, and as a consequence the calculated CO₂ flux was overestimated. Therefore, the three levels of ICOS CO₂ concentration measurements proved to be vital in our flux calculations here. We still missed the logarithmic profile at the surface with only the two vertical campaign measurements, and as a result we slightly underestimated the final CO₂ and O₂ flux. Therefore, we recommend to always measure at least three heights of CO₂ and O₂ inside the surface layer when they are meant to be used for flux calculations.

Our final O₂ flux (Fig. 6) shows a clear diurnal cycle, with the expected behaviour of negative values in the night (O₂ consumption for respiration) and a positive flux during the day (O₂ release during assimilation). The nighttime fluxes are more stable and give a clear signal due to the larger vertical gradient. K is more difficult to determine during the night as the EC measurements are less representative due to the low level of turbulence. However, the largest contributor to the uncertainty is our own O₂ measurements, and the larger gradient allows us to better establish the O₂ flux. The larger variability of the daytime O₂ fluxes is caused by the smaller gradient of the O₂ concentration measurements during the day (Fig. 3), when the atmosphere is more well mixed and when the difference between the two heights becomes smaller. The relatively large measurement uncertainty made it difficult to measure these small differences between the two heights and increased the noise in the fluxes. The measurement noise resulted in O₂ gradient variations that were not tied to the CO₂ gradient variations, and this degraded the correlation between the two fluxes. Despite this larger variability, we still find a clear diurnal behaviour, which allowed us to calculate ER_forest. Note that the uncertainties of the surface fluxes of O₂ and CO₂ are only based on the measurements from our campaigns, and we did not include the uncertainties that are related to the calculations of K based on the ICOS data. However, the uncertainty in K is relatively small compared to the other terms in the calculation, and the final uncertainty of estimates is dominated by the measurement uncertainty of O₂. Omitting the uncertainty associated with K therefore leads to a minor underestimate of the full uncertainty.

4.4 ER_forest signal compared to previous studies

Our resulting ER_forest signal changes throughout the diurnal cycle, with specific daytime (0.92 ± 0.17), nighttime (1.03 ± 0.05), and overall (0.84 ± 0.26) values (Fig. 7 and Table 3). The individual nighttime values show a smaller uncertainty due to the already explained effect of the larger gradient during the stable atmospheric conditions of the night. In contrast, the individual daytime values show a larger uncertainty due to the smaller gradient during the unstable atmospheric conditions of the day. We therefore used averaged values for the daytime and nighttime signals to derive the ER_forest values. While the daytime signal excludes the entrainment and the boundary layer dynamics during the morning transition, these effects are still included in the overall ER_forest signal. Note that the overall 24 h signal is not the average of the daytime and nighttime signal. The nighttime ER_forest signal represents a negative O₂ flux and a positive CO₂ flux, whereas the daytime ER_forest signal represents a positive O₂ flux and a negative CO₂ flux. This means that the daytime and nighttime surface fluxes influence the atmosphere differently; therefore, these ER_forest values cannot be averaged to calculate the overall ER_forest signal. By first calculating the average overall O₂ and CO₂ fluxes and then dividing these, we derive the overall ER_forest signal correctly.

When comparing our ER_forest signals to previous studies by Seibt et al. (2004), Ishidoya et al. (2013), and Ishidoya et al. (2015) (Table 4), we notice that the difference between the daytime and the nighttime values that we found and the specific values of the different ER_forest has some similarities and some differences. Our results, along with those by Ishidoya et al. (2013, 2015) (night: 1.11 and day: 0.98), show that the ER_forest signal of the nighttime is higher than the daytime signal, whereas Seibt et al. (2004) (day: 1.24 and night: 1.01) show the opposite behaviour. Our results are most similar to the signals of both Ishidoya et al. (2013) and Ishidoya et al. (2015), especially if we take our uncertainty range into account. The 24 h signals are difficult to compare as we used a different method to determine the overall ER_forest signal compared to Ishidoya et al. (2015). In this study we use average fluxes instead of a linear regression through either the O₂ and CO₂ fluxes or vertical gradient, and we thereby take into account the size of the fluxes that contributes most to the ER signal, which results in a flux-weighted average ER_forest. We note again that we need to distinguish between daytime and nighttime signals, and we cannot just average them. Figure 7 illustrates the need to take averages in consistent meteorological and biological periods that are characterized by similar turbulence regimes and similar signs of the O₂ and CO₂ exchange. For example, combining a small negative O₂ flux with a high ER with a large positive O₂ flux with a lower ER results in a smaller O₂ flux compared to when the ERs of both fluxes would have been the same. When we take into account our uncertainty, the complete day signal of 0.84 ± 0.26 comes close to the globally used average ER of the biosphere of 1.1 (Severinghaus, 1995). However, the specific value suggests that the overall ER_forest signal of this boreal forest lies somewhat lower than 1.1, i.e. closer to 1.0. The difference in ER_forest signals between studies can be explained with the different ER_a and ER_r signals, which we discuss in Sect. 4.5.

The ER_forest and ER_atmos signals are not identical; therefore, they do not represent the same information (Table 3). The ER_atmos signals are higher compared to the ER_forest signals; the 24 h signals especially show a large difference. Despite the higher numbers, the day and night signals of ER_atmos and ER_forest both show the same pattern, where the daytime signal is lower compared to the nighttime signal. When comparing these differences to previous studies, we find that not all studies find the same results. The difference between ER_forest and ER_atmos was not found by Ishidoya et al. (2013). In contrast, Seibt et al. (2004) found a difference between ER_forest and ER_atmos (Table 4). A reason for this could be the measurement height of ER_atmos. When ER_atmos is determined closer to the canopy and inside the roughness sublayer, it will be more influenced by the surface processes compared to measurements at higher levels, which are seeing more integrated signals of all processes that influence the concentrations inside the atmospheric boundary layer (i.e. forest exchange and non-local processes like entrainment). To get a clear answer to this question, we should further investigate to what extent ER_atmos is influenced by entrainment and boundary layer dynamics and under which conditions they can come close to ER_forest. We already show that excluding the morning transition (P2) helps to improve the ER_atmos signal. However, as already stated, it is difficult from the measurements alone to determine if the ER_atmos signal is influenced only by the surface during this period. An atmospheric model would therefore be needed to find how ER_atmos can be derived from a single measurement height and to allow for comparison with previous studies that measured at one height to determine the ER of the forest (Battle et al., 2019; van der Laan et al., 2014; Stephens et al., 2007). We are currently applying a specific mixed-layer atmospheric model to further investigate this.

4.5 The ER_a and ER_r signals

To further understand the relationship between O₂ and CO₂, we cannot use the ER_forest signal alone. To look in more detail into the processes driving the variations, we calculated the exchange ratios of respiration (ER_r) and assimilation (ER_a) (Table 3 and Fig. 7). ER_r was taken as the ER_forest nighttime signal (1.03 ± 0.05), by assuming that only respiration influences the ER_forest signal during the night and that the ER_r signal stays constant throughout the entire day. This means that both the heterotrophic and autotrophic respiration are included in ER_r, and the same components are respired in the same ratios throughout the day to keep ER_r a constant value. The studies that looked at the ER_r of an ecosystem (Hilman et al., 2022; Angert et al., 2015; Hicks Pries et al., 2020) only focused on longer timescales than the diurnal cycle. It is therefore not possible to derive diurnal variability of ER_r from these previous studies. We would expect some changes in ER_r as temperatures change during the day, and the respiration of plants involves photorespiration during daytime and dark respiration during nighttime. However, as the study by Hilman et al. (2022) showed, the variability of ER_r mainly depends on the bulk soil respiration and therefore depends on the soil temperature and soil moisture. No large changes in temperatures or soil moisture were detected during the period of the representative aggregate day; therefore, it is unlikely that the ER_r significantly changed in that period. To get a more detailed view on how the ER_r of an ecosystem changes throughout the day, more research is needed on the variability of ER_r including from plant respiration by chamber measurements.

The variability of ER_r between locations highly depends on the soil properties (Angert et al., 2015), which makes it difficult to compare with the few studies available (Seibt et al., 2004; Ishidoya et al., 2013) that measured ER_r with chamber measurements on a brown soil. The soil in our study area is a podzol, which is characterized by a high acidity with little organic matter (Buurman and Jongmans, 2005). The OR of podzols is around 1.08 (Worrall et al., 2013), and the ER of acid soils is expected to be around this OR, because carbon cannot easily dissolve into the groundwater (Angert et al., 2015), and we therefore conclude that our ER_r value of 1.04 is realistic.

We looked at two options to calculate ER_a: the ER_a based only on the daytime measurements (between 09:00 and 17:00: 0.96 ± 0.12) and ER_a based on all the measurements (throughout the 24 h period: 0.96 ± 0.11). Both numbers are close to 1, which is often assumed as a standard value for ER_a (Ishidoya et al., 2015; Severinghaus, 1995). Next to that, a value of ER_a close to 1 means that ammonium is used as a source for nitrogen instead of nitrate (Bloom et al., 1989, 2012). Ammonium is indeed a larger source for nitrogen compared to nitrate in Hyytiälä (Korhonen et al., 2013). The OR of needle leaves (and plant material in general) appears to be always close to 1.0 (Jürgensen et al., 2021), which again confirms our ER_a signals. We did not observe differences between the two ER_a signals. The transition periods between the night and the daytime were difficult to measure, because the gradient then becomes close to zero, which means there could be a possibility that next to ER_forest, ER_a also has a diurnal cycle. Again, there are only a few studies that looked at the variability of ER_a (Fischer et al., 2015; Bloom et al., 1989, 2012; Bloom, 2015). The available studies show that ER_a depends on light (Fischer et al., 2015) and the source of nitrogen in the soil (Bloom, 2015). These changes in ER_a happen when the changes in the atmosphere and the soil are sudden and persist for a longer time compared to a diurnal cycle. We can therefore say that the ER_a also does not change drastically during the day. To get a more detailed overview of ER_a, more precise measurements are needed with uncertainties lower than 0.1 for the ER signals. However, the similar values of ER_a that we find for the daytime and 24 h measurements show that ER_a is hardly affected by entrainment during the morning transition, and it would suggest that the morning transition is less of an issue for ER_forest than for ER_atmos.

By applying the O₂ method to a new aggregate day, we showed that the O₂ method gives results similar to the EC method to partition NEE and derive the GPP and TER fluxes (Fig. 8), with estimates of the uncertainties of the O₂ method. The EC method also contains uncertainties in its approach because of the reliance on a function of temperature and should therefore not necessarily be assumed to be the truth (Reichstein et al., 2005). Despite the uncertainty of both the O₂ method and the EC method, both methods give similar results for the CO₂ flux of GPP and TER. In our comparison of the O₂ and EC methods, there is a minor degree of circularity, as we use the EC GPP estimates to estimate ER_a. By applying it to another representative day, we prevent a full circularity (Table B1). In this campaign, we unfortunately could not determine the ER_a and ER_r signals independently from EC, which would be recommended for a full comparison. This would have been possible by using chamber measurements. We expect only minor changes in ER_a from branch/leaf chamber measurements compared to the values we derived, because our ER_a is in the range of expected values compared to previous studies (e.g. Jürgensen et al., 2021; Bloom, 2015; Fischer et al., 2015). The satisfactory comparison between the O₂ and the EC methods for the partitioning of the fluxes shows the potential of the O₂ method. The largest challenge for this method is to determine ER_forest with large enough accuracy, as this value is most variable and most difficult to determine based on the small O₂ gradients that we observed. Figure 8 shows that the ER_forest signal should be measured with an uncertainty of 0.05 or less to get results within the uncertainty range of the EC method. When such high accuracy is reached, the O₂ method has the potential to provide an alternative method for the separation of GPP and TER without relying on the regularly used temperature-based function as used for the EC method. Ishidoya et al. (2015) showed similar results, where the O₂ method also produced GPP and TER comparable to the EC method and the magnitude of the GPP and TER fluxes highly dependent on the derived ER_a and ER_r signals.

To allow for an independent comparison between the flux partitioning with the EC method and the O₂ method, such as was done by using δ¹³CO₂ by Wehr et al. (2016), we recommend measuring the ER_r and ER_a signals directly with chamber measurements (Seibt et al., 2004; Ishidoya et al., 2013). We also recommend to add at least one additional measurement height for the O₂ and CO₂ concentrations below the canopy to apply the storage correction for both the O₂ and the CO₂ fluxes (Aubinet et al., 2012) and to add a measurement in the free troposphere to better evaluate the effect of entrainment. Despite the high dependency on the accuracy of the ER, this study showed again, as did Ishidoya et al. (2015), that the O₂ method can be used to get a better understanding of the carbon cycle. To further develop this method, we need to expand the O₂ measurements for longer time series and more locations and to analyse how ER_forest varies over longer timescales, which can improve the global average value of ER (α_B) of 1.1 as used in global carbon budget studies such as Manning and Keeling (2006).

5 Conclusions

By continuously measuring atmospheric O₂ and CO₂ concentrations at two heights in Hyytiälä, Finland, we gained new insights into the diurnal variability of O₂ and CO₂ above a boreal forest, quantified by interpreting their exchange ratio (ER). We showed that the signal based on one measurement height of the O₂ and CO₂ concentrations (ER_atmos) is not representative for the exchange between the forest and the atmosphere only, but it instead includes other processes such as entrainment as well. To derive the ER of the forest specifically (ER_forest), we first determined the surface fluxes above the canopy of O₂ and CO₂ using the vertical gradient between the two measurement heights. We found that the most suitable method to calculate both the O₂ and CO₂ surface fluxes was to use the exchange coefficient calculated from the eddy covariance (EC) CO₂ flux and the vertical gradient of CO₂ measurements at three heights above the canopy. The ER_forest signals that resulted from the ratio of the mean O₂ and CO₂ fluxes varied between the daytime (0.92 ± 0.17 mol mol⁻¹) and nighttime (1.03 ± 0.05 mol mol⁻¹). The different ER_forest signals were composed of the ER of respiration (ER_r: 1.03 ± 0.05 mol mol⁻¹) and the ER of assimilation (ER_a: 0.96 ± 0.12 mol mol⁻¹). With these findings, we show improved methods to derive O₂ forest fluxes and to derive the variability in the different ER signals over a representative diurnal cycle. The ER_forest signal shows a clear diurnal cycle for this boreal forest, and the overall ratio is lower than 1.1 that is used in global carbon budget calculations. Finally, we show that these ER signals can be used to separate net ecosystem exchange (NEE) into gross primary production (GPP) and total ecosystem respiration (TER).

With only a few data sets of continuous measurements of both O₂ and CO₂ concentrations over forests, our data set is of high importance, specifically the availability of measurements at two heights that allow for calculation of O₂ and CO₂ fluxes. Our analyses can serve as a starting point for follow-up research using coupled land surface–atmosphere models to distinguish and quantify contributions of different processes to ER_atmos and ER_forest signals. Further understanding of these differences will help to fully exploit the advantages of atmospheric O₂ when unravelling the different components in the carbon cycle.

Appendix A: Equations to calculate the exchange coefficient, K

A1 Observation-based method

The gradient between three points is calculated with the following equation:

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{\mathit{\varphi }}\left(z\right)=a\cdot z^{\mathrm{2}}+b\cdot z+c,\text{(A2)}& {\displaystyle }& {\displaystyle }\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}\left(z\right)}{\partial z}}\right)=\mathrm{2}\cdot a\cdot z+b,\end{array}$$

where z [m] is the height above the displacement height (d [m]) (d is taken as $\mathrm{2}/\mathrm{3}\cdot$ canopy height); $\stackrel{\mathrm{‾}}{\mathit{\varphi }}$ is the average variable where the line is fitted through; and a, b, and c are the resulting fitted parameters. When only two vertical measurements are available, the gradient was determined using finite differences.

A2 Theoretical approach

For the MOST method, the following equations were used (Physick and Garratt, 1995):

$$\begin{array}{}\text{(A3)}& K={\displaystyle \frac{\mathit{\kappa }\cdot z\cdot u_{\ast }}{{\mathrm{\Phi }}_{\mathrm{H}}\left(\frac{z}{L}\right){\mathit{\varphi }}_{\text{rsl}}\left(\frac{z}{L}\right)}},\end{array}$$

where K is the exchange coefficient [m² s⁻¹], κ = 0.4 is the von Kármán constant, u_∗ [m s⁻¹] is the friction velocity, Φ_H [–] indicates the stability function, and ϕ_rsl [–] indicates the contribution of the roughness sublayer (RSL). The Φ_H was calculated with the following (Dyer, 1974):

$$\begin{array}{}\text{(A4)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{H}}\left({\displaystyle \frac{z}{L}}\right)={\left(\mathrm{1}-\mathrm{16}{\displaystyle \frac{z}{L}}\right)}^{-\mathrm{1}/\mathrm{2}}& {\displaystyle }\text{ when }z/L<\mathrm{0},\text{(A5)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{H}}\left({\displaystyle \frac{z}{L}}\right)=\mathrm{1}+\mathrm{5}{\displaystyle \frac{z}{L}}& {\displaystyle }\text{ when }z/L>\mathrm{0},\end{array}$$

where L [m] is the Obukov length, which was based on the following equation (Dyer, 1974):

$$\begin{array}{}\text{(A6)}& L={\displaystyle \frac{-u_{\ast }^{\mathrm{3}}}{\mathit{\kappa }\left(\frac{g}{{\mathit{\theta }}_{\mathrm{v}}}\right)\left(\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}\right)}},\end{array}$$

where θ_v [K] is the virtual potential temperature, $\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}$ [K m s⁻¹] is the virtual surface heat flux, and g [m s⁻²] is the acceleration due to gravity. Because the flux was measured close to the canopy, the roughness surface layer (RSL) could become important. The RSL needs an additional length scale (ϕ) and can be calculated with the following equation (de Ridder, 2010):

$$\begin{array}{}\text{(A7)}& {\mathit{\varphi }}_{\mathrm{H}\text{RSL}}\left({\displaystyle \frac{z}{z_{\ast }}}\right)=\mathrm{1}-e^{-\mathit{\mu }\frac{z}{z_{\ast }}}.\end{array}$$

Here z_∗ [m] indicates the height of the RSL above the displacement height and we take that as (2 ⋅ canopy height − d), and μ is a constant of 0.95 [–].

By integrating Eq. (7) with Eq. (A3) for K, we get the following equation that was used for the integrated method (Physick and Garratt, 1995):

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}\mathit{\varphi }\left(z_{\mathrm{2}}\right)-\mathit{\varphi }\left(z_{\mathrm{1}}\right)=& {\displaystyle \frac{\left(\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\varphi }}^{\prime }}\right)}{\mathit{\kappa }\cdot u_{\ast }}}\left[\ln \left({\displaystyle \frac{z_{\mathrm{2}}}{z_{\mathrm{1}}}}\right)-{\mathrm{\Psi }}_{\mathrm{H}}\left({\displaystyle \frac{z_{\mathrm{2}}}{L}}\right)\right.\\ & \left.+{\mathrm{\Psi }}_{\mathrm{H}}\left({\displaystyle \frac{z_{\mathrm{1}}}{L}}\right)+{\mathit{\psi }}_{\text{RSL}}\left({\displaystyle \frac{z}{L}},{\displaystyle \frac{z}{z_{\ast }}}\right)\right],\end{array}\end{array}$$

where Ψ_H [–] represents the integrated stability functions for heat, and ψ [–] is the integrated function to account for the roughness sublayer (RSL) effect. Ψ_H was calculated with (Paulson, 1970)

$$\begin{array}{}\text{(A9)}& {\displaystyle }& {\displaystyle }\left.\begin{array}{rl}{\displaystyle }{\mathrm{\Psi }}_{\mathrm{H}}\left({\displaystyle \frac{z}{L}}\right)& {\displaystyle }=\mathrm{2}ln\left({\displaystyle \frac{\mathrm{1}+x^{\mathrm{2}}}{\mathrm{2}}}\right)\\ {\displaystyle }x& {\displaystyle }=(\mathrm{1}-\mathrm{16}z/L{)}^{\mathrm{1}/\mathrm{4}}\end{array}\right\}& {\displaystyle }\text{ when }z/L<\mathrm{0},\text{(A10)}& {\displaystyle }& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{H}}\left({\displaystyle \frac{z}{L}}\right)=-\mathrm{5}{\displaystyle \frac{z}{L}}& {\displaystyle }\text{ when }z/L>\mathrm{0}.\end{array}$$

The function of the integrated RSL length scale (ψ_RSL) [–] was calculated with (de Ridder, 2010)

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}{\mathit{\psi }}_{\text{RSL}}\left({\displaystyle \frac{z}{L}},{\displaystyle \frac{z}{z_{\ast }}}\right)\approx & {\mathrm{\Phi }}_{\mathrm{H}}\left[\left(\mathrm{1}+{\displaystyle \frac{\mathit{\nu }}{\mathit{\mu }z/z_{\ast }}}\right){\displaystyle \frac{z}{L}}\right]\\ & \cdot {\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}\ln \left(\mathrm{1}+{\displaystyle \frac{\mathit{\lambda }}{\mathit{\mu }z/z_{\ast }}}\right)e^{-\mathit{\mu }z/z_{\ast }},\end{array}\end{array}$$

where ν and λ are both parameters, taken as 0.5 and 1.5, respectively.

Appendix B: Figures and tables

Figure B1 The corrected and uncorrected half-hourly average O₂ concentrations at 23 m for the PMKS deviations (a), together with the corrected and uncorrected gradient of O₂ (b). The time series for the 2019 measurement period of the PMKS (c) and the air temperature (d) are also shown. The shaded areas indicate the dates that were selected for the two aggregate representative days, i.e. 13 through 15 June (grey) and 7 through 12 July 2019 (yellow).

Table B1 Here are the data used to calculate different variables that we calculated in our study for the two aggregate days. Part (a) indicates the data that were used for the first representative aggregate day during 7 through 13 July 2019, and part (b) indicates the data that were used for the second aggregate day between 13 and 15 June 2019. The data used are mostly from the period of the respective aggregate day, except when indicated otherwise.

n/a = not applicable

Figure B2 The footprint of the O₂ and CO₂ surface fluxes at 53.6 m height (which is the geometric height), determined with the gradient method for the days during 7 through 12 July 2019 at Hyytiälä. The lines and contours indicate the contributions to the footprint from 10 % to 90 % in steps of 10 %. The plus sign (+) indicates the location of the tower. This figure was created with the method by Kljun et al. (2015).